#include "mrg_misc3.inc"
global_settings
{
  assumed_gamma 1.0
  ambient_light .2
}                                          
background{color rgb 0}

#declare lclock=clock+.0;

#declare cscale = 5;
#declare cpos =<0.0, 0, -5.0>*1.4*cscale; 
camera {
  location  cpos
  look_at   <0.0, -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
}
//macro to calculate the area of overlap between two circles of radii R1, R2>=R1 centers seperated by delter
#macro set_eclipsed_area(Area,R1,R2,delter)
     #local dd=abs(delter);
     #local range1=(R2-R1);
     #local range2=sqrt(R2*R2-R1*R1);
     #local range3=(R1+R2);
     #declare Area=-1;
     #switch (dd)
          #range (0,range1)
               #declare Area=pi*R1*R1;
          #break
          #range (0,range2)  //this is for xm <=0
               #local xm=(R1*R1-R2*R2+dd*dd)/(2*dd);
               #local aa=pi*R1*R1;
               #local aa=aa-(R1*R1*acos(-xm/R1)+xm*sqrt(R1*R1-xm*xm));
               #local aa=aa+R2*R2*acos((dd-xm)/R2)-(dd-xm)*sqrt(R2*R2-(dd-xm)*(dd-xm));
               #declare Area=aa;
          #break
          #range (0,range3)  //this is for xm >= 0
               #local xm=(R1*R1-R2*R2+dd*dd)/(2*dd);
               #local aa=R1*R1*acos(xm/R1)-xm*sqrt(R1*R1-xm*xm);
               #local aa=aa+R2*R2*acos((dd-xm)/R2)-(dd-xm)*sqrt(R2*R2-(dd-xm)*(dd-xm));
               #declare Area=aa;
          #break
          #else
          #declare Area=0;

     #end
#end
#declare rpm=5;
//orbital plane axis, tilting
//#declare tilta=function(tt){sin(2*pi*tt)*12}
#declare tilta=function(tt){max(int(rpm*tt),0)*12/rpm}



#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 p1_0=-d1*x;
#declare p2_0=d2*x;
#declare p1=vrotate (p1_0,rpm*360*y*lclock);
#declare p1=vrotate (p1, tilta(lclock)*x);
#declare p2=vrotate (p2_0,rpm*360*y*lclock);
#declare p2=vrotate (p2, tilta(lclock)*x);
//Area Spline: used because
#declare tdd=0;
#declare dtd=.005;
#declare Areas = function{
     spline {
          quadratic_spline 
          #while(tdd<1.4*(r1+r2))
               #declare Area=-1;
               set_eclipsed_area(Area,r1,r2,tdd)
               tdd,<Area,0,0>
               #declare tdd=tdd+dtd;
          #end     
          
     }
}

#macro L_tot(tt)
    // #local tp1=vrotate (p1_0,rpm*360*orb_axis(tt)*tt);
     //#local tp2=vrotate (p2_0,rpm*360*orb_axis(tt)*tt);

     #local tp1=vrotate (p1_0,rpm*360*y*tt);
     #local tp1=vrotate (tp1, tilta(tt)*x);
     #local tp2=vrotate (p2_0,rpm*360*y*tt);
     #local tp2=vrotate (tp2, tilta(tt)*x);
     #local deltv=(tp1-tp1.z*z)-(tp2-tp2.z*z);
     #local delt=vlength(deltv);
     #local dz=tp1.z-tp2.z;
     #declare tArea = 0;
     #local Lm=0;
     #if(dz>0)
       #local Lm=L1*Areas(delt).x/(pi*r1*r1);
     #end
     #if(dz<0)
       #local Lm=L2*Areas(delt).x/(pi*r2*r2);
     #end
     #declare Ltot=(L1+L2-Lm)/(L1+L2);
#end






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}
          rotate x*tilta(lclock)
          }
     torus{vlength(p2),.05
          pigment{color rgbt star_color(m2)+<0,0,0,.6>}
          finish{ambient 4*pow(m2,.25) diffuse 0}
          rotate x*tilta(lclock)
          }

     scale  .2
}



//plot the macro


union{
     // plot
     #declare lth=.015;
     #declare pt=0;
     #declare dpt=.0005;
     #declare vscale = 1;
     #declare hscale=2.5;
     #declare tscale=2.5*rpm;
     #declare Area = -1;
     #declare apt=pt+lclock;
     #declare Ltot=0;
     L_tot(apt)
     #declare this=<pt*tscale,Ltot,0>;
     //marker for current brightness
     cylinder{0,-.1*z,3*lth
          pigment{color rgb <.8,1,.3>} finish{ambient 4 diffuse 0}  
          translate this-3*lth*x     
          }
     union{
          #while(pt>-1.0/rpm)
               #declare pt=pt-dpt;
               #declare apt=(pt+lclock);
               #declare Ltot=0;
               L_tot(apt)
               #declare next=<pt*tscale,Ltot,0>;
               cylinder{this, next,lth pigment{color rgb <1,0,0>} finish{ambient 3}}
               #declare this=next;
          #end 
     }
     cylinder{+.0*x,-(hscale+.0)*x,lth pigment{color rgb 1}finish{ambient 5 diffuse 0}}     
     text { ttf "arial.ttf","net Luminosity",.02, 0
          pigment{color rgb <1,.7,.5>} finish{ambient 5 diffuse 0}  
          scale .18  rotate z*90
          translate +.17*x+(vscale/2-.45)*y -.01*z
     }
         
     box{-.0*x-(.25)*y+.1*z,(hscale+.0)*x+(vscale+.25)*y+.05*z 
          pigment {
               brick color rgb <0,1,0>, color rgb 0
               brick_size 1/4 mortar .025
               rotate x*90 translate .5*y
               translate -x*tscale*lclock
          }
          finish{ambient 2}
          translate -hscale*x
     }
     box{-.125*x-(.35)*y+.15*z,(hscale+.25)*x+(vscale+.35)*y+.35*z 
          pigment {color rgb <0,0,1> }
          finish{ambient 2}
          translate -hscale*x
     }
     scale 1.25
     translate -3.5*y+.3*x
}