#include "mrg_misc3.inc"
global_settings
{
  assumed_gamma 1.0
  ambient_light .2
}                                          
background{color rgb 0}

#declare lclock=clock+.25;

#declare cscale = 100;
#declare cpos =<0.0, 2, -5.0>*12*cscale+20*x; 
camera {
  location  cpos
  look_at   <25, -1.2,  0.0>
  right     x*image_width/image_height
  angle 60/cscale
}
//stellar color


#declare star_color = spline{//adapted from main sequence data table on wikipedia
        linear_spline
//   M/M  pc          class	R/R       L/L 		T
     158,<.239,.420,1>//	O2  	16 		   ?     	54,000
     58,<.243,.424,1>//	O5	14		   ?      46,000
     16,<.267,.447,1>//	B0	5.7		16,000	29,000
     5.4,<.341,.522,1>//	B5	3.7		750	     15,200
     2.6,<.486,.647,1>//	A0	2.3		63	     9,600
     1.9,<.612,.741,1>//	A5	1.8		24	     8,700
     1.6,<.694,.8,1	>//  F0	1.6	     9.0	     7,200
     1.35,<.831,.894,1>//F5	1.2		4.0	     6,400
     1.08,<.93,.957,1>//	G0	1.05		1.45	     6,000
     1.0,<.992,.996,1>//	G2	1.0		1.0	     5,900
     0.95,<1,.965,.914>//G5	0.98		0.70	     5,500
     0.83,<1,.914,.796>//K0	0.89		0.36	     5,150
     0.62,<1,.796,.569>//K5	0.75		0.18 	4,450
     0.47,<1,.682,.384>//M0	0.64		0.075	3,850
     0.25,<1,.541,.22>//	M5	0.36		?	     3,200
}

#declare rpm=1;



#declare tts=0;
#declare dtt=.005;


//binary star physical parameters
/*
a	     b	     c
0.1925	1.527	-0.9486
polynomial fit for m=a*r^2+b*r+c
*/
#declare r_of_m = function(mm) { (-1.527+sqrt(1.527*1.527-4*0.1925*(-0.9486-mm)))/(2*0.1925)  };
#declare L_of_m = function(mm) { pow(mm,3.5)  };

#declare m1=.8; //star masses, always have m1<=m2!
#declare m2=1.5;
#declare r1=r_of_m(m1);  //stellar radii
#declare r2=r_of_m(m2);
#declare L1=L_of_m(m1);  //stellar luminosities
#declare L2=L_of_m(m2);
#declare d=5*(r1+r2);      //separation
#declare d1=m2*d/(m1+m2);  //about center of mass
#declare d2=m1*d/(m1+m2);
#declare x1_t=function(tt){ -d1*cos(2*pi*rpm*tt)}
#declare z1_t=function(tt){ -d1*sin(2*pi*rpm*tt)}
#declare v1_t=function(tt){2*pi*rpm*d1*sin(2*pi*rpm*tt)}
#declare vmax=2*pi*rpm*d1;
#declare x2_t=function(tt){ d2*cos(2*pi*rpm*tt)}
#declare z2_t=function(tt){ d2*sin(2*pi*rpm*tt)}
#declare v2_t=function(tt){-2*pi*rpm*d2*sin(2*pi*rpm*tt)}
#declare p1=<x1_t(lclock),0,z1_t(lclock)>;
#declare p2=<x2_t(lclock),0,z2_t(lclock)>;
#declare c = 2.0*2*pi*d1; //speed of waves, set to 3* sources max speed






union{
     sphere{p1,r1
          pigment{color rgb star_color(m1)}
          finish{ambient 4*pow(m1,.25) diffuse 0}
          }
     sphere{p2,r2
          pigment{color rgb star_color(m2)}
          finish{ambient 4*pow(m2,.25) diffuse 0}
          }
     torus{vlength(p1),.05
          pigment{color rgbt star_color(m1)+<0,0,0,.6>}
          finish{ambient 4*pow(m1,.25) diffuse 0}
          }
     torus{vlength(p2),.05
          pigment{color rgbt star_color(m2)+<0,0,0,.6>}
          finish{ambient 4*pow(m2,.25) diffuse 0}
          }

}
#declare dx=.05;

#declare f=20;
#declare sig=function(tt){1*cos(f*2*pi*tt)}
#declare xmax=40;
//sets sig_out to be sig at point xx at time tt, emitted at tp  1-D wave please
#macro sig1_rec(xx,tt,tp1,sig1_out)
     #local told=tt;  //told =first guess when g=xx-x_t(tp)-c*(tt-tp)=0
     #local gold=xx-x1_t(told)-c*(tt-told);
     #local tnew=tt-.1;  //tnew =next guess when g=xx-x_t(tp)-c*(tt-tp)=0
     #local gnew=xx-x1_t(tnew)-c*(tt-tnew);
     #declare tp1=tnew;
     #while (abs(tnew-told)>.0001)
          #declare tp1=tnew-gnew*(tnew-told)/(gnew-gold);
          #local told=tnew;
          #local gold=gnew;
          #local tnew=tp1;
          #local gnew=xx-x1_t(tnew)-c*(tt-tnew);        
     #end
     #declare sig1_out=sig(tp1);
#end

#declare xplot=x1_t(lclock);
#declare tsig=0;
#declare te1=0;
sig1_rec(xplot,lclock,te1,tsig);
#declare first1=<xplot,tsig,0>;
#while (xplot<xmax)
     #declare xplot=xplot+dx;
     #declare tsig=0;
     #declare te1=0;
     sig1_rec(xplot,lclock,te1,tsig);
     #declare next1=<xplot,tsig,0>;
     cylinder{first1, next1,.08 
          pigment{color rgb <.5-.5*v1_t(te1)/vmax,.5+.5*v1_t(te1)/vmax*(1-v1_t(te1)/vmax),.5+.5*v1_t(te1)/vmax>} 
          finish{ambient 3}
          translate z1_t(lclock)*z
     }
     #declare first1=next1;     
#end

//sets sig_out to be sig at point xx at time tt, emitted at tp  1-D wave please
#macro sig2_rec(xx,tt,tp,sig2_out)
     #local told=tt;  //told =first guess when g=xx-x_t(tp)-c*(tt-tp)=0
     #local gold=xx-x2_t(told)-c*(tt-told);
     #local tnew=tt-.1;  //tnew =next guess when g=xx-x_t(tp)-c*(tt-tp)=0
     #local gnew=xx-x2_t(tnew)-c*(tt-tnew);
     #declare tp=tnew;
     #while (abs(tnew-told)>.0001)
          #declare tp=tnew-gnew*(tnew-told)/(gnew-gold);
          #local told=tnew;
          #local gold=gnew;
          #local tnew=tp;
          #local gnew=xx-x2_t(tnew)-c*(tt-tnew);        
     #end
     #declare sig2_out=sig(tp);
#end
#declare xplot=x2_t(lclock);
#declare tsig=0;
#declare te2=0;
sig2_rec(xplot,lclock,te2,tsig);
#declare first=<xplot,tsig,0>;
#while (xplot<xmax)
     #declare xplot=xplot+dx;
     #declare tsig=0;
     #declare te2=0;
     sig2_rec(xplot,lclock,te2,tsig);
     #declare next=<xplot,tsig,0>;
     cylinder{first, next,.08 
          pigment{color rgb <.5-.5*v2_t(te2)/vmax,.5+.5*v2_t(te2)/vmax*(1-v2_t(te2)/vmax),.5+.5*v2_t(te2)/vmax> } 
          finish{ambient 3}
          translate z2_t(lclock)*z
     }
     #declare first=next;     
#end

cylinder{(xmax-0)*x,60*x, d
          pigment{color rgb 1 } 
          finish{ambient 3}
     scale <1,.25,1>
}

// create a TrueType text shape
text {
  ttf             // font type (only TrueType format for now)
  "crystal.ttf",  // Microsoft Windows-format TrueType font file name
  "spectrometer",      // the string to create
  .2,              // the extrusion depth
  0               // inter-character spacing
          pigment{color rgb <0,.5,0> } 
          finish{ambient 3}
     scale 2.5
     translate -d*z+xmax*x+2*x+7*y
}
union{
     box{-6*x,6*x+6*y+.01*z}
     cylinder{.2*y,5.8*y,.1 
          pigment{color rgb <.5-.5*v1_t(te1)/vmax,.5+.5*v1_t(te1)/vmax*(1-v1_t(te1)/vmax),.5+.5*v1_t(te1)/vmax> } 
          finish{ambient 3}
          translate v1_t(te1)/vmax*5*x-.1*z
     
     }
     cylinder{.2*y,5.8*y,.1 
          pigment{color rgb <.5-.5*v2_t(te2)/vmax,.5+.5*v2_t(te2)/vmax*(1-v2_t(te2)/vmax),.5+.5*v2_t(te2)/vmax> } 
          finish{ambient 3}
          translate v2_t(te2)/vmax*5*x-.1*z
     
     }
     translate -d*z+xmax*x+9*x
}     
     