#include "math.inc"
#include "colors.inc"    // The include files contain
#include "stones.inc"    // pre-defined scene elements
#include "finish.inc"
#include "functions.inc"
#include "geodesic_math.inc"
#declare lclock=clock+0.;
#declare lc4= lclock*lclock*lclock*lclock;
global_settings {
  assumed_gamma 1.0
  ambient_light rgb <1, 1, 1>
}

background{Gray}
// create a regular point light source
light_source {
  0*x                  // light's position (translated below)
  color rgb <1,1,1>    // light's color
  translate <20, 40, 0>
}


#declare surf_size=3;
isosurface {
  function { fns(x, y, z) }        // alternative declared function
  contained_by { box { -surf_size, surf_size } }  // container shape
  threshold 0.0                     // optional threshold value for isosurface [0.0]
  accuracy 0.001                      // accuracy of calculation [0.001]
  max_gradient 2.556                      // maximum gradient the function can have [1.1]
  //evaluate 5, 1.2, 0.95             // evaluate the maximum gradient
  //max_trace 1                       // maybe increase for use in CSG [1]
  //all_intersections                 // alternative to 'max_trace'
  open                              // remove visible container surface
  texture{pigment{color rgb .6*<.2,.2,.8>} }
}
#declare Sl=.2+1.2*lclock; //total length of geodesic
#declare geowidth=.01;
#declare ng = 20; //number of geodesic segments
#declare ds = Sl/ng; //cylinder length for drawing geodesic
#declare nsub =3;  // number of subdividions to help integration accuraccy
#declare dds=ds/nsub;

#macro drawgeodesic(q1,q1v,q2,q2v)
        #declare iseg=0;
        #while(iseg<ng)
                #declare q1n=q1;
                #declare q2n=q2;
                #declare isubseg=0;
                #while (isubseg<nsub)
                        //acceleration from geodesic equation
                        #declare q1a=-(C111(q1n,q2n)*q1v*q1v+C112(q1n,q2n)*q1v*q2v+C121(q1n,q2n)*q2v*q1v+C122(q1n,q2n)*q2v*q2v);
                        #declare q2a=-(C211(q1n,q2n)*q1v*q1v+C212(q1n,q2n)*q1v*q2v+C221(q1n,q2n)*q2v*q1v+C222(q1n,q2n)*q2v*q2v);
                        //increment coordinate from approximate cons acceleration
                        #declare q1n=q1n+q1v*dds+.5*q1a*dds*dds;        
                        #declare q2n=q2n+q2v*dds+.5*q2a*dds*dds;
                        //increment velocity
                        #declare q1v=q1v+q1a*dds;
                        #declare q2v=q2v+q2a*dds;
                        #declare isubseg=isubseg+1;        
                #end
                cylinder { <xs(q1,q2), ys(q1,q2), zs(q1,q2)>, <xs(q1n,q2n), ys(q1n,q2n), zs(q1n,q2n)>, geowidth
                                  texture{pigment{color rgb .9*<.2,.8,.8> } finish{ambient 1}}
                                  no_shadow
                                  }
                #declare q1=q1n;
                #declare q2=q2n;
                #declare iseg=iseg+1;               
        #end
#end

#declare seperation=.1;

#declare q1=-seperation;    //initial data
#declare q1v=0;
#declare q2=.06;
#declare q2v=1;
drawgeodesic(q1,q1v,q2,q2v)


#declare q1=seperation;    //initial data
#declare q1v=0;
#declare q2=.06;
#declare q2v=1;
drawgeodesic(q1,q1v,q2,q2v)

#declare q1=0;             //initial data
#declare q1v=0;
#declare q2=.06;
#declare q2v=1;
drawgeodesic(q1,q1v,q2,q2v)

//where to point the camera:
#declare xc= xs(q1n,q2n);
#declare yc= ys(q1n,q2n);
#declare zc= zs(q1n,q2n);

/*for pos curvature */
camera {
    location 2*<4, 2*lclock+2*lc4, 0>
    look_at  <xc,yc*(1-lc4)/2,zc*(1-.25*lc4*lc4)/2>
    angle 12+18*lclock
}
/*  for neg curvature
camera {
    location (1.8+.0*lc4)*<1.6*lc4*lc4, 2+3*lclock+1.2*lc4, -.05-1.5*lc4>
    look_at  <xc,yc,zc*(1-.25*lc4*lc4)/2>
    angle 14+34*lclock+0*lc4
}
*/

//cross pieces\
#declare cwidth=geowidth;   //width of markers
#declare cds=.2;                //spacing of cross markers
#declare cl= 6*seperation;      //length of cross markers
#declare ncds = 5;              //number of integration steps per tic
#declare csnow=0;               //initialize cross marker location
#declare q1=0;             //initial data
#declare q1v=0;
#declare q2=.06;
#declare q2v=1;
#declare dds=cds/ncds;
#declare q1n=q1;
#declare q2n=q2;
#while(csnow<Sl)
        #declare isubseg=0;
        #declare velocity=<x_1(q1n,q2n)*q1v+x_2(q1n,q2n)*q2v,
                                y_1(q1n,q2n)*q1v+y_2(q1n,q2n)*q2v,
                                z_1(q1n,q2n)*q1v+z_2(q1n,q2n)*q2v>;   //CHAIN RULE
        #declare grad=<fns_x(xs(q1n,q2n),ys(q1n,q2n),zs(q1n,q2n)),
                        fns_y(xs(q1n,q2n),ys(q1n,q2n),zs(q1n,q2n)),
                        fns_z(xs(q1n,q2n),ys(q1n,q2n),zs(q1n,q2n))>;  //gradient is parallel to normal
                        //problems with grad or velocity near "steep" points of curve
        #declare d_cend=vnormalize(vcross(velocity,grad))*cl/2;
        #declare here=<xs(q1n,q2n),ys(q1n,q2n),zs(q1n,q2n)>;
        #declare surface_offset=vnormalize(grad)*cwidth*(-0);   //"lift" cross markers out of surface a bit for visibility
        #declare l_end=here-d_cend+surface_offset;
        #declare r_end=here+d_cend+surface_offset;
        union{                                                //alternating colors on segments for cross markers
                cylinder{(6*l_end+0*r_end)/6,(5*l_end+1*r_end)/6,cwidth}                    
                cylinder{(4*l_end+2*r_end)/6,(3*l_end+3*r_end)/6,cwidth}                    
                cylinder{(2*l_end+4*r_end)/6,(1*l_end+5*r_end)/6,cwidth}                    
                texture{pigment{color rgb <1,1,1> } finish{ambient 1}}
                no_shadow
             }
        union{
                cylinder{(5*l_end+1*r_end)/6,(4*l_end+2*r_end)/6,cwidth}                    
                cylinder{(3*l_end+3*r_end)/6,(2*l_end+4*r_end)/6,cwidth}                    
                cylinder{(1*l_end+5*r_end)/6,(0*l_end+6*r_end)/6,cwidth}                    
                texture{pigment{color rgb <0,0,0> } finish{ambient 1}}
                no_shadow
             }
        #while (isubseg<ncds)
                //acceleration from geodesic equation
                #declare q1a=-(C111(q1n,q2n)*q1v*q1v+C112(q1n,q2n)*q1v*q2v+C121(q1n,q2n)*q2v*q1v+C122(q1n,q2n)*q2v*q2v);
                #declare q2a=-(C211(q1n,q2n)*q1v*q1v+C212(q1n,q2n)*q1v*q2v+C221(q1n,q2n)*q2v*q1v+C222(q1n,q2n)*q2v*q2v);
                //increment coordinate from approximate cons acceleration
                #declare q1n=q1n+q1v*dds+.5*q1a*dds*dds;        
                #declare q2n=q2n+q2v*dds+.5*q2a*dds*dds;
                //increment velocity
                #declare q1v=q1v+q1a*dds;
                #declare q2v=q2v+q2a*dds;
                #declare isubseg=isubseg+1;        
        #end
        #declare csnow=csnow+cds;
#end