// Persistence of Vision Ray Tracer Scene Description File
// File: ?.pov
// Vers: 3.5
// Desc: Basic Scene Example
// Date: mm/dd/yy
// Auth: ?
//

#version 3.5;

#include "colors.inc"
#include "stars.inc"
#include "shapes.inc"
#include "mrg_misc3.inc"

global_settings {
  assumed_gamma 1.0
 max_trace_level 12

}

/*   
sphere {0,200
 hollow
  texture {
    Starfield1
  }
}
*/
#declare cel_sphr_rad=40;
#declare stars=union{
     #declare numb=0;
     #declare dum=0;
     #declare ra=0;
     #declare dec=0;
     #declare mag=0;
     #fopen MyFile "stars_by_mag.csv" read
     #read (MyFile,num,ra,dec,mag)
     #while (mag<=6)
          sphere { <0, 0, -cel_sphr_rad>, 0.0010*pow(1-mag/7,.33)*cel_sphr_rad
               texture{
                    pigment{color rgb 1*(1-mag/7)}
                    finish{ambient 2 diffuse 0}
               }
               rotate x*dec //rotate -z towards y for "altitude" like declination
               rotate -y*ra //right ascension in a left handed coord system...
                
          }   
          #read (MyFile,num,ra,dec,mag)
     
     #end
     #fclose MyFile
};
object{stars}
camera {
  location  <0.0,+ 1, -3.0>*5.0
  direction 1.5*z
  right     x*image_width/image_height
  look_at   <0.0, 0.0,  0.0>
  angle 12
}

light_source {
  <0, 0, 0>            // light's position (translated below)
  color rgb <1, 1, 1>  // light's color
  translate <-30, 30, -30>
}


#declare lclock=clock+12;
#declare q=lclock*360*0;
#declare T1=ltrans((lclock-1.0)/1);  //fade in wells
#declare T2=ltrans((lclock-2.0)/1);  //fade in contours
#declare T3=ltrans((lclock-4.0)/1);  //fade out wells/ contours  and fade in lobe markers
#declare T4=ltrans((lclock-5.0)/1);  //fade in lobes
#declare T5=trans((lclock-7.0)/3);  //expand star2
#declare T6=ltrans((lclock-10.0)/3);  //mass exchange init 
#declare DT6=3;
#declare T7=ltrans((lclock-13.0)/10);  //mass exchange  


  //old script
#declare t1=trans((lclock-0)/3);//expand to red giant    
#declare t1p5=ltrans((lclock-3.0)/3);  //mass exhange flow
#declare t2=ltrans((lclock-6)/10); 
  //old script parameters revised
#declare t1=T5;
#declare t1p5=T6;
#declare t2=T7;

   
#declare m1_0=1;
#declare m2_0=4;
#declare mT=m1_0+m2_0;
#declare r1_0=1;
#declare r2_0=m1_0/m2_0*r1_0;
#declare l_0=r1_0+r2_0;
#declare w_0=1;

#declare m1=m1_0*(1-t2)+m2_0*t2;
#declare m2=m2_0*(1-t2)+m1_0*t2;
#declare l=l_0*(m1_0*m1_0*m2_0*m2_0)/(m1*m1*m2*m2);
#declare w=w_0*pow(l_0/l,1.5);
#declare r1=m2/mT*l;
#declare r2=m1/mT*l;


#declare mf1=sqrt(m1);
#declare mf2=sqrt(m2);
#declare xb=(mf1*r2-mf2*r1)/(mf1+mf2);
#declare uinv=function{sqrt((x+r1)*(x+r1)+y*y+z*z)*sqrt((x-r2)*(x-r2)+y*y+z*z)/
     (m2*sqrt((x+r1)*(x+r1)+y*y+z*z)+m1*sqrt((x-r2)*(x-r2)+y*y+z*z))};
#declare p_thresh3=uinv(xb,0,0);
#declare try3=function { uinv(x,y,z)-p_thresh3};
     

//"flow"
#if(t1p5>0)
     cylinder{0,t1p5*(-r1-xb)*x,.010 translate xb*x
      texture {
            pigment {          // (---surface color---)
              marble           // some pattern
              color_map {      // color map
                [0.1 color rgbt <1,0,0,.7> ]
                [0.75 color rgbt <1,0,0,1>]
              }
              turbulence 1.5   // some turbulence
              scale <3,1,1>    // transformations
            }
            normal {           // (---surface bumpiness---)
              marble 0.3       // some pattern (and intensity)
              turbulence 1.5   // some turbulence
              scale <3,1,1>    // transformations
            }
            finish {           // (---surface finish---)
              ambient 0.8
              specular 0.0     // not shiny
            }
            scale .10
            translate -(r1+xb)/DT6*lclock //flow velocity
     }

           
           no_shadow no_reflection
     }
#end


/*countours   */
#if(T2-T3>0)  
     intersection{
          isosurface {
               function { uinv(x,y,z)-uinv(xb,0,1.0)}        accuracy 0.01
                       max_gradient 2
                       contained_by{sphere{0,3}}
                       //finish {phong 0.5 phong_size 100}
          }
          box{3*(x+z),-3*(x+z)-.001*y}
          pigment { function {20*try3(x,y,z) }
              color_map { [0.0   rgbt <.3,.3,.3,1-(T2-T3)>]
                          [0.08   rgbt 1]
                          [0.92   rgbt 1]
                          [1.0   rgbt <.3,.3,.3,1-(T2-T3)>]
              }
          }
          finish{ambient 1 diffuse .5}
          no_shadow no_reflection
          rotate y*q
     }
#end

#if(T4>0)
     /* left lobe */
       
     isosurface {
          function { try3(x,y,z)-.0001}        accuracy 0.0051
             max_gradient 1
             contained_by{sphere{-r1*x,abs(xb+r1)}}
             hollow
          pigment{color rgbt <.7,.7,1,.9+.1*(1-T4)>} no_shadow no_reflection
              scale <1.00,1.04,1.04>
     
         rotate y*q
     } 
     /**/
      //right lobe        
     isosurface {
          function { try3(x,y,z)-.0001}        accuracy 0.0051
             max_gradient .5
             contained_by{sphere{r2*x,abs(r2-xb)}}
             hollow
          pigment{color rgbt <1,.7,.7,.9+.1*(1-T4)>}  no_shadow no_reflection
         rotate y*q
         scale <1.00,1.04,1.04>
     }
#end

/**/ 
//left star-simple sphere
sphere{-r1*x,.05*pow(m1_0*m1/m1_0,.333)
     pigment{color rgb <1,1,1.2*t2>}
     finish{ambient .5 diffuse .2}
      no_shadow no_reflection
    rotate y*q
}

/*
#declare s1_thresh=uinv(-r1,.05*pow(m1*m1/m1_0,.33),0); 
isosurface {
     function { uinv(x,y,z)-s1_thresh}        accuracy 0.00025
        max_gradient .9
        contained_by{sphere{-r1*x,abs(xb+r1)}}
     pigment{color rgb Yellow}
     finish{ambient .5 diffuse .2}
      no_shadow no_reflection
    rotate y*q
}
*/ 
//right star-simple sphere
sphere{r2*x,.05*pow(m2_0*m2/m2_0,.333)
     pigment{color rgb <1,1,1.2>}
     finish{ambient .5 diffuse .2}
      no_shadow no_reflection
    rotate y*q
}

//right star
/*#declare s2_thresh=uinv(r2,.05*pow(m2*m2/m2_0,.33),0); 
isosurface {
     function { uinv(x,y,z)-s2_thresh}        accuracy 0.00025
        max_gradient .89
        contained_by{sphere{r2*x,abs(r2-xb)}}
        hollow
     pigment{color rgb White}
          finish{ambient .5 diffuse .2}
      no_shadow no_reflection
    rotate y*q
}
*/
#declare nf=50;
#declare tl=.4;
#declare s3_thresh=uinv(r2,.05*m2/m2_0,0)*(1-t1)+t1*uinv(xb,0,0); 


/*

intersection{
     isosurface {
          function { uinv(x,y,z)-s3_thresh +.005*f_noise3d(x*nf,y*nf,z*nf)}        
             accuracy 0.00025
             max_gradient .4
             contained_by{sphere{r2*x,abs(r2-xb)}}
     }     
     box{3*(-x-y-z),3 }             
     texture { 
               pigment {function {(uinv(x,y,z)-uinv(r2,0,0))/(uinv(xb,1.4,0)-uinv(r2,0,0))}
                   color_map { 
                              [0 rgbt <.6,.6,.6,tl>]
                              [0.2 rgbt <.7,.1,.1,tl>]
                              [.5 rgbt <.5,.08,.05,tl>]
                              [.7 rgbt <.5,.05,.0,tl>]
                              [1 rgbt <.3,0,0,tl>]
     
                   }
               }
               finish{ambient 1.0 diffuse .0}
      
     }
     no_reflection no_shadow
}
*/


 //gas inflation 
#if(T5>0)
     isosurface {
          function { uinv(x,y,z)-s3_thresh+.0000}        
          accuracy 0.00025
             max_gradient .52
             contained_by{sphere{r2*x,abs(r2-xb)}}
          hollow
          interior{ 
               media{ 
                    emission 1
                    density{ 
                         function { uinv(x,y,z)-s3_thresh}
                          density_map{ 
                              [0 rgb <.7,.7,1>]
                              [0.8 rgb <.7,.1,0>]
                              [.999 rgb <.5,.0,0>]
                              [.9999 rgb <.3,.0,0>]
                              [1 rgb <0,0,0>]
                         }
                    }
               }
               media{scattering{1,<.3,0,0>}}
          }
     
             texture { 
                     pigment { color rgbt <1,1,1,1> } 
                     finish{ambient .0 diffuse .0}
      
            }
            no_reflection no_shadow
     }
#end
/* */

/* coordinate axis
cylinder{-x,x,.01 pigment{color rgb <.2,.2,.6>}  no_shadow no_reflection} 
cylinder{-y,y,.01 pigment{color rgb <.2,.2,.6>}  no_shadow no_reflection} 
cylinder{-z,z,.01 pigment{color rgb <.2,.2,.6>}  no_shadow no_reflection} 

 */
 
 
//lobe demarcation
#if(T3>0)
     #declare u_f=function(x,z){uinv(x,0,z)-uinv(xb,0,0)} ;     //2-D potential
     #declare uz_f=function(s){u_f(s,0)} ;      //1-D potential
     #declare uz_fp=function(s){(uz_f(s+1e-4)-uz_f(s))*1e4} ;  //potential derivative
     #declare u_fx=function(xx,zz) {(u_f(xx+1e-6,zz)-u_f(xx,zz))*1e6}   //f_x
     #declare u_fz=function(xx,zz) {(u_f(xx,zz+1e-6)-u_f(xx,zz))*1e6}   //f_z
     #declare LLcolor=<0,0,1,1-T3>; 
     #declare RLcolor=<1,1,0,1-T3>; 
      /* */
     #declare pr=.005; //lobe perimeter radius 
     #declare ds=10;
     #declare sold=-(r1+abs(r1-abs(xb)));
     #declare nn=1;
     #while(ds>.001&nn<10)
          #declare snew=sold-uz_f(sold)/uz_fp(sold);
          #declare ds=abs(snew-sold);
          #declare nn=nn+1; 
          #declare sold=snew;
     #end
     #declare p_dl=.005;
     #declare pold=<snew,0,0>;
     #declare fold=<u_fx(pold.x,pold.z),0,u_fz(pold.x,pold.z)>;
     #declare pnew=pold+p_dl*vnormalize(vcross(y,fold));
     //#declare ds_c=vnormalize(<u_fx(pold.x,pold.z),0,u_fz(pold.x,pold.z)>);//direction of "correction"
     //#declare pnew=pnew+u_f(pold.x,pold.z)/(u_fx(pold.x,pold.z)*ds_c.x-u_fz(pold.x,pold.z)*ds_c.z)*ds_c;
     cylinder{pold, pnew,pr  pigment{color rgbt LLcolor} no_shadow no_reflection}
     cylinder{pold, pnew,pr  pigment{color rgbt LLcolor} scale <1,1,-1>no_shadow no_reflection}
     #while(pnew.x<=xb)
          #declare pold=pnew;
          #declare fnow=<u_fx(pold.x,pold.z),0,u_fz(pold.x,pold.z)>;
          #declare pnew=pold+p_dl*vnormalize(vcross(y,fnow));
          cylinder{pold, pnew,pr  pigment{color rgbt LLcolor} no_shadow no_reflection}
          cylinder{pold, pnew,pr  pigment{color rgbt LLcolor} scale <1,1,-1>no_shadow no_reflection}
     #end
     
     #declare ds=10;
     #declare sold=+(r2+abs(r2-abs(xb)));
     #declare nn=1;
     #declare nn=1;
     #while(ds>.001&nn<10)
          #declare snew=sold-uz_f(sold)/uz_fp(sold);
          #declare ds=abs(snew-sold);
          #declare nn=nn+1; 
          #declare sold=snew;
     #end
     
     #declare pold=<snew,0,0>;
     #declare fnow=<u_fx(pold.x,pold.z),0,u_fz(pold.x,pold.z)>;
     #declare pnew=pold+p_dl*vnormalize(vcross(y,fnow));
     cylinder{pold, pnew,pr  pigment{color rgbt RLcolor} no_shadow no_reflection}
     cylinder{pold, pnew,pr  pigment{color rgbt RLcolor} scale <1,1,-1>no_shadow no_reflection}
     #while(pnew.x>=xb)
          #declare pold=pnew;
          #declare fnow=<u_fx(pold.x,pold.z),0,u_fz(pold.x,pold.z)>;
          #declare pnew=pold+p_dl*vnormalize(vcross(y,fnow));
          cylinder{pold, pnew,pr  pigment{color rgbt RLcolor} no_shadow no_reflection}
          cylinder{pold, pnew,pr  pigment{color rgbt RLcolor} scale <1,1,-1>no_shadow no_reflection}
     #end
#end 

#if(T1-T3>0)
     //contour
        
     isosurface {
          function {y+log(1/uinv(x,0,z))}        accuracy 0.01
             max_gradient 54 
             open
             contained_by{sphere{0,2}}
          pigment { rgbt <.7,.7,.7,.8+.2*(1-T1+T3)> }
          no_shadow no_reflection
         rotate y*q
         translate .3*y
     }    
#end