Antialias=on

Antialias_Threshold=0.2
Antialias_Depth=3
Input_File_Name=rotation_curve_planet.pov
Output_File_Name=rotation_curve_planets

Initial_Frame=1
Final_Frame=320   
Initial_Clock=.0
Final_Clock=1

Cyclic_Animation=on
Pause_when_Done=off
#version 3.5;

#include "colors.inc"
#include "mrg_misc3.inc"
#include "planets.inc"
global_settings {
  assumed_gamma 1.0
  ambient_light color rgb 1 //<1,.8,.7>*.7
}
background{color rgb .0}

// ultra_wide_angle lens for wide, rectangular field of view


//The Physics

//acceleration of gravity as function of location... we'll keep this in x-y plane
/*
//this particular choice for g(x,y) should reproduce the galactic rotation curves
//#declare grav = function(xx,yy){1/sqrt(xx*xx+yy*yy)}

//this particular choice for g(x,y) is for point mass
#declare grav = function(xx,yy){1/(xx*xx+yy*yy)};

//this particular choice for g(x,y) is for inside uniform mass density
//#declare grav = function(xx,yy){sqrt(xx*xx+yy*yy)};

//speed for circular orbit at radius R
#declare Vorb=function(R){sqrt(R*grav(R,0))};


//derivative functions
#declare fx=function(xx,yy,vx,vy){vx};
#declare fy=function(xx,yy,vx,vy){vy};
#declare fvx=function(xx,yy,vx,vy){-xx*grav(xx,yy)/sqrt(xx*xx+yy*yy)};
#declare fvy=function(xx,yy,vx,vy){-yy*grav(xx,yy)/sqrt(xx*xx+yy*yy)};

//macro to evolve coordinates  //4th order Runge Kutta
#macro position_update(mx,my,mvx,mvy,tf)
     #local mdt=.01*tf;
          //initial conditions
     #local x1=mx;
     #local y1=my;
     #local vx1=mvx;
     #local vy1=mvy;
     //loop
     #local mt=0;
     #while(mt<tf)
          #local k1x=fx(x1,y1,vx1,vy1)*mdt;
          #local k1y=fy(x1,y1,vx1,vy1)*mdt;
          #local k1vx=fvx(x1,y1,vx1,vy1)*mdt;
          #local k1vy=fvy(x1,y1,vx1,vy1)*mdt;

          #local k2x=fx(x1+.5*k1x,y1+.5*k1y,vx1+.5*k1vx,vy1+.5*k1vy)*mdt;
          #local k2y=fy(x1+.5*k1x,y1+.5*k1y,vx1+.5*k1vx,vy1+.5*k1vy)*mdt;
          #local k2vx=fvx(x1+.5*k1x,y1+.5*k1y,vx1+.5*k1vx,vy1+.5*k1vy)*mdt;
          #local k2vy=fvy(x1+.5*k1x,y1+.5*k1y,vx1+.5*k1vx,vy1+.5*k1vy)*mdt;

          #local k3x=fx(x1+.5*k2x,y1+.5*k2y,vx1+.5*k2vx,vy1+.5*k2vy)*mdt;
          #local k3y=fy(x1+.5*k2x,y1+.5*k2y,vx1+.5*k2vx,vy1+.5*k2vy)*mdt;
          #local k3vx=fvx(x1+.5*k2x,y1+.5*k2y,vx1+.5*k2vx,vy1+.5*k2vy)*mdt;
          #local k3vy=fvy(x1+.5*k2x,y1+.5*k2y,vx1+.5*k2vx,vy1+.5*k2vy)*mdt;

          #local k4x=fx(x1+k3x,y1+k3y,vx1+k3vx,vy1+k3vy)*mdt;
          #local k4y=fy(x1+k3x,y1+k3y,vx1+k3vx,vy1+k3vy)*mdt;
          #local k4vx=fvx(x1+k3x,y1+k3y,vx1+k3vx,vy1+k3vy)*mdt;
          #local k4vy=fvy(x1+k3x,y1+k3y,vx1+k3vx,vy1+k3vy)*mdt;

          #local x2=x1+(k1x+2*k2x+2*k3x+k4x)/6;
          #local y2=y1+(k1y+2*k2y+2*k3y+k4y)/6;
          #local vx2=vx1+(k1vx+2*k2vx+2*k3vx+k4vx)/6;
          #local vy2=vy1+(k1vy+2*k2vy+2*k3vy+k4vy)/6;

          #local x1=x2;
          #local y1=y2;
          #local vx1=vx2;
          #local vy1=vy2;

          #local mt=mt+mdt;
     #end
     //use global rather than local variables to send out the update values

     #declare xout=x2;
     #declare yout=y2;
     #declare vxout=vx2;
     #declare vyout=vy2;
#end
*/

//test
/*
#declare xold=1;
#declare yold=0;
#declare vx0=0;
#declare vy0=Vorb(1)*1.2;
#declare dt=.3;
#declare tst=dt;
#declare drw=1;
#while(tst<30)

     position_update(xold,yold,vx0,vy0,dt);
     cylinder{<xold,yold,0>,<xout,yout,0>,.0125
          pigment{color rgb <1,0,1>}
     }
     #declare xold=xout;
     #declare yold=yout;
     #declare vx0=vxout;
     #declare vy0=vyout;

     #declare tst=tst+dt;
#end

object{sun scale .2}
*/
//endtest
#declare orbit_scale=4;
#declare radiusT_scale=20;
#declare radiusJ_scale=8;

#declare planet_object = array[10];
#declare planet_position = array[10];
#declare planet_omega = array[10];
#declare planet_phase = array[10];

#declare planet_object[0]=object {sun scale sun_display*.7 no_shadow};
#declare planet_position[0]=<0,0,0>;
#declare planet_omega[0]=0;

#declare planet_object[1]=object {mercury scale mercury_display*radiusT_scale no_shadow};
#declare planet_position[1]=orbit_scale*mercury_orbit*x;
#declare planet_omega[1]=1/mercury_period;

#declare planet_object[2]=object {venus scale venus_display*radiusT_scale no_shadow};
#declare planet_position[2]=orbit_scale*venus_orbit*x;
#declare planet_omega[2]=1/venus_period;

#declare planet_object[3]=object {earth scale earth_display*radiusT_scale no_shadow};
#declare planet_position[3]=orbit_scale*earth_orbit*x;
#declare planet_omega[3]=1/earth_period;

#declare planet_object[4]=object {mars scale mars_display*radiusT_scale no_shadow};
#declare planet_position[4]=orbit_scale*mars_orbit*x;
#declare planet_omega[4]=1/mars_period;

#declare planet_object[5]=object {jupiter scale jupiter_display*radiusJ_scale no_shadow};
#declare planet_position[5]=orbit_scale*jupiter_orbit*x;
#declare planet_omega[5]=1/jupiter_period;

#declare planet_object[6]=object {saturn scale saturn_display*radiusJ_scale no_shadow};
#declare planet_position[6]=orbit_scale*saturn_orbit*x;
#declare planet_omega[6]=1/saturn_period;

#declare planet_object[7]=object {uranus scale uranus_display*radiusJ_scale no_shadow};
#declare planet_position[7]=orbit_scale*uranus_orbit*x;
#declare planet_omega[7]=1/uranus_period;

#declare planet_object[8]=object {neptune scale neptune_display*radiusJ_scale no_shadow};
#declare planet_position[8]=orbit_scale*neptune_orbit*x;
#declare planet_omega[8]=1/neptune_period;

#declare planet_object[9]=object {pluto scale pluto_display*radiusT_scale*8 no_shadow};
#declare planet_position[9]=orbit_scale*pluto_orbit*x;
#declare planet_omega[9]=1/pluto_period;

#declare R1 = seed(39);      
#declare pcount=0;
#while (pcount<10)
     #declare planet_phase[pcount]=360*rand(R1);
     #declare pcount=pcount+1;
#end
#declare planet_phase[7]=60;
#declare planet_phase[8]=90;
#declare planet_phase[9]=100;

#declare lclock=clock+.0;
#declare pcount=0;
#while (pcount<10)
     #declare pangle=-(planet_phase[pcount]+360*12*lclock*planet_omega[pcount])*y;

     object{planet_object[pcount] 
          translate planet_position[pcount]
          rotate pangle
     }
     #declare pcount=pcount+1;
#end
#declare cscale=2.3;
#declare ds=jupiter_orbit*orbit_scale;
#declare cpos=<0, .2*ds,-ds>*cscale;

// perspective (default) camera
camera {
  location  cpos 
  look_at   <0.0, 0.0,  0.0>
  right     x*image_width/image_height
  angle 70/cscale                          // field
}
/*
camera {
  location cpos               // position
  look_at  <0, 0, 0>                // view
  right x*image_width/image_height  // aspect
  angle 70/cscale                          // field
}*/
#declare bpos=cpos*.9;
union{
     box{-.2*x-.2*y,x+y+.01*z
          translate .01*z
          pigment{color rgb 1}
          finish{ambient 1.0 diffuse 0}
     }
     text { ttf "timesbd.ttf", 
          "speed", .01,  0
          rotate z*90
          scale .18
          translate -z*.01-.05*x+.25*y 
          pigment{color rgb <1,0,0>}
          finish{ambient 1 diffuse 0}
     }
     text { ttf "timesbd.ttf", 
          "distance", .01,  0
          scale .18
          translate -z*.01+.15*x-.15*y 
          pigment{color rgb <1,0,0>}
          finish{ambient 1 diffuse 0}
     }
     #declare vscale=planet_omega[1]*vlength(planet_position[1]);//max speed for merc
     #declare rscale=vlength(planet_position[9]);//max dist, go with nept
     #declare np=1;
     #while (np<=9)
          #declare xnew=vlength(planet_position[np])/rscale;
          #declare ynew=planet_omega[np]*vlength(planet_position[np])/vscale;
          //data points
          cylinder{-.01*z,.01*z,.025
               translate <xnew,ynew,0>
               pigment{color rgb z} finish{ambient .7 diffuse 0}
          }
          //connect the dots
          #if(np >1)
               cylinder{<xold,yold,0>,<xnew,ynew,0>,.01
                    pigment{color rgb z} finish{ambient .9 diffuse 0}
               }
          #end
          #declare xold=xnew;
          #declare yold=ynew;
          #declare np=np+1;
     #end
     cylinder{-.15*x,x,.01    pigment{color rgb 0}}
     cylinder{-.15*y,y,.01    pigment{color rgb 0}}
          rotate +x*12
          scale .6   translate -.90*y-.25*x-.95*x
          translate bpos
}
/**/
light_source {
  <0, 0, 0>            
  color rgb  <1,.8,.7>*2 //brighten up those planets 
}

// ----------------------------------------

#declare T2=1;
#declare cel_sphr_rad=50;
#declare stars=union{
     #declare ra=0;
     #declare dec=0;
     #declare mag=0;
     #declare con="aaa";
     #fopen MyFile "decimal_yale.csv" read
     #read (MyFile,ra,dec,mag,con)
     #declare star_origin=cel_sphr_rad*(z);
     #declare altitude_axis=vcross(star_origin,y);
     #declare altitude_axis=altitude_axis/vlength(altitude_axis);
     #while (mag<=6)
          sphere { star_origin, 0.0025*(1-mag/7)*cel_sphr_rad
               texture{
                    pigment{color rgb 1*(1-mag/10)}
                    finish{ambient 20 diffuse 0}
               }
               rotate altitude_axis*dec //rotate towards y for "altitude" like declination
               rotate -y*ra //right ascension in a left handed coord system...
                
          }   
          #read (MyFile,ra,dec,mag,con)
     
     #end
     #fclose MyFile
          #declare cline_rad=0.0055*.5*cel_sphr_rad;

};


//places winter solstice at -z, summer solstice at +z, autumn equinox at -x, spring equinox at +x
object { stars  rotate x*23}