/* Program file=raysph.c                        */
/*  Raytracing program for spheres  
                         cited from "minimal ray tracer"
*/

#include <stdio.h>
#include <math.h>
#include "raysph_gl.h"
typedef struct { double x,y,z; } vec;

/* set parameters */
#define TOL 1e-7
#define SIZE 250        /* resolution of picture in x and y */
#define AOV 35          /* total angle of view in degrees */
#define NSPHERE 2       /* number of spheres */

struct sphere {vec cen, color; double rad,kd,ks,deg;} \
               *s, *best, sph[] \
     = {0., 0., 0.,  1., 0.1, .1,   1.3,   0.6, 0.5, 30.0,  
      -2., -1., -.9, 0., 1.0, 0.2,  1.4,   .3,  0.6, 90.,
 }; /* sphere: (x y z) (r g b)  rad  kd ks  deg */
vec U, light,origin, back = {.1, .1, .8};   /* back-groud color */
double the_light=-50., pi_light=30.;       /* light direction */
double ambient = 0.15;                     /* ambient light */
double xv=0., yv=0., zv=0.;                /* look at */ 
double r_view=6.0, the_view=0., pi_view=0.; /* look from */

double u, b, tmin;
double ct, st, cf, sf;

/* ************* vector normarizing ***************************: */
vec vunit(vec A)
{
 double d;
     d = sqrt(A.x*A.x + A.y*A.y + A.z*A.z);
     A.x = A.x/d; A.y = A.y/d; A.z = A.z/d;
     return A;
}
/* ************* intersection test for ray/spheres *************: */
struct sphere *intersect(vec P, vec D)
{
    best = 0;
    tmin = 1e30;
    s = sph+NSPHERE;
    while (s-- >sph)
       {U.x=s->cen.x-P.x;U.y =s->cen.y-P.y;U.z=s->cen.z- P.z;
        b = D.x*U.x + D.y*U.y + D.z*U.z;
        u = b*b-(U.x*U.x + U.y*U.y+U.z*U.z)+s->rad*s->rad;
        u = u>0 ? sqrt(u) : 1e31;
        u = b-u>TOL ? b-u : b+u;
        tmin = u>=TOL && u<tmin ?
            best = s, u : tmin;
        }
    return best;
}
/* ****** ray tracing ****************************************: */
vec trace( vec P, vec D)
{
    vec N,  H, color;
    struct sphere *s;
    double  gamma,cos_alfa,specular, diffuse;

    s = intersect(P, D);
    if (s==0) return back;
    /* P(x,y,z): Intersection point between ray and sphere */
    P.x=tmin*D.x+P.x; P.y =tmin*D.y+P.y;P.z=tmin*D.z + P.z;
    N.x=-P.x+s->cen.x;N.y =-P.y+s->cen.y;N.z=-P.z+ s->cen.z;
    N = vunit(N); /* N(x,y,z): normal vetcor at P */
    cos_alfa = light.x*N.x + light.y*N.y + light.z*N.z;
    specular = 0.;
    if(cos_alfa>0.) 
      { diffuse = cos_alfa*s->kd + ambient;
        H.x=light.x + D.x; H.y =light.y+D.y; H.z=light.z + D.z;
       gamma = (H.x*N.x + H.y*N.y + H.z*N.z)/sqrt(H.x*H.x + H.y*H.y + H.z*H.z);
        if(gamma>0.) specular = s->ks*pow(gamma,s->deg);
    }
    else diffuse = ambient;

    color.x = s->color.x*diffuse + specular; 
    color.y = s->color.y*diffuse + specular; 
    color.z = s->color.z*diffuse + specular;
    return color;
}
/* ***** calculate coefficents for peresepective transformation ****: */
void preper(double the,double fai)
{
        double fef,tef;

        fef=3.1415/180.*fai; tef=3.1415/180.*the;
        ct=cos(tef);     st=sin(tef);
        cf=cos(fef);     sf=sin(fef);
}
/* ************ perspective transformation ********************: */
void trans(double xt,double yt,double zt,double *px,double *py,double *pz)
{
      double xy;

      xy=ct*xt+st*yt;         *pz=(r_view-cf*xy-sf*zt);
      *px=(ct*yt-st*xt);      *py=(cf*zt-sf*xy);
}
/* ******* transformation from world-cood. to eye-cood. *******: */
void tranfrm()
{
      preper(the_view,pi_view);
      s = sph+NSPHERE;
      while (s-->sph){
        trans(s->cen.x-xv,s->cen.y-yv,s->cen.z-zv, \
        &s->cen.x, &s->cen.y, &s->cen.z);
      }
      trans(light.x,light.y,light.z,&light.x,&light.y,&light.z);
      light.z =  light.z - r_view;
}
/* ************** light vector ******************************: */
vec lightvec(double the_light, double pi_light)
{
     vec L;
     L.x = -cos(pi_light*3.14/180.)*cos(the_light*3.14/180.);
     L.y = -cos(pi_light*3.14/180.)*sin(the_light*3.14/180.);
     L.z = -sin(pi_light*3.14/180.);
     return L;
}
/* ********** main function *********************************: */
void main()
{
   int x,y;
   vec COL,RAY;
   double scale;

   light = lightvec(the_light,pi_light);
   tranfrm();
   scale = SIZE/2/tan(AOV/114.5915590261)/r_view;  

   SETSCREEN(SIZE,SIZE)

   /* scan start */
                             RAY.z = r_view;
   for(y=0; y<=SIZE; y++)  { RAY.y = (SIZE/2 - y)/scale;      
      for(x=0; x<SIZE; x++){ RAY.x = (x - SIZE/2)/scale; 
        COL = trace(origin, vunit(RAY)); /* ray tracing */
        PLOT(x,SIZE-y,COL.x,COL.y,COL.z)
      }   /* loop for pixel */
   }      /* loop for scan-line */
   /* scan end */
   GEND
}