#version unofficial Megapov 0.7

//#render "glow {}\n"
//#render "test\n"

//#render concat(" ",str(step/maxstep*100,1,0),"% of spiral ",str(num+1,1,0)," ",chr(13))   

/*
#declare new_per_frame=80;
#declare max_age=20;
#declare max_particles=max_age*new_per_frame;
#declare particle_position=array[max_particles]
#declare particle_speed=array[max_particles]
#declare particle_age=array[max_particles]
#declare pointer=0;
#declare frame=frame_number;
#declare step=1/200; 
#declare tmax=max_age*step;
#declare initial_time=comet_clock-tmax-step;
#declare time=initial_time; 
#declare maxv=2;
#declare gas=maxv/2;
#declare wind=5;
#declare factor=0.9;
*/

#declare new_per_frame=80;
#declare new_per_frame2=80;
#declare max_age=12;
#declare max_age2=12;
#declare max_particles=max_age*new_per_frame;
#declare max_particles2=max_age2*new_per_frame2;
#declare particle_position=array[max_particles]
#declare particle_speed=array[max_particles]
#declare particle_age=array[max_particles]
#declare particle_position2=array[max_particles2]
#declare particle_speed2=array[max_particles2]
#declare particle_age2=array[max_particles2]
#declare pointer=0;
#declare pointer2=0;
#declare frame=frame_number;
#declare step=1/25; 
#declare tmax=max_age*step;
#declare initial_time=comet_clock-tmax;
#declare time=initial_time; 
#declare maxv=0.5;
#declare maxv2=1.5;
#declare gas=maxv/2;
#declare wind=3;
#declare wind2=1.5;
#declare factor=0.9;
#declare factor2=0.85;

#declare rnd=seed(0);

#declare sun_position=<-200,0,0>;

/*
#macro comet_position(time)
  //#declare pos=<time*2,0,0>;
  
  //#declare pos=vrotate(<-1,0,0>,<0,360*time,0>);
  //#declare pos=pos-<0.5,0,0>;
  //#declare pos=pos*<2,1,1>;
  //#declare pos=pos*75;
  //pos
  
  #declare st=0;
  #declare e=0.8;
  #declare a=1; 
  #declare GM=2;
  #declare theta=-0.8*pi;
  #declare xpos=0;
  #declare ypos=0;
  #while (st<=time)
    #declare rr=a*(1-e*e)/(1+e*cos(theta));
    #declare dtheta=(1-e*e)*step/rr/rr*GM*a;
    #declare theta=theta+dtheta;
    #declare rr=a*(1-e*e)/(1+e*cos(theta));
    #declare xpos=rr*cos(theta);
    #declare ypos=rr*sin(theta);
    #declare st=st+step; 
    //sphere { <xpos,ypos,0>,0.02 finish {ambient 1} pigment {color rgb <1,1,1>} }
    #debug concat("st=",str(st,1,2),"xpos=",str(xpos,1,2)," ypos=",str(ypos,1,2),"\n")
  #end
  //camera { location <0,0,-10> look_at <0,0,0> angle 30}
  //plane { z,0 finish { ambient 1} pigment {checker }}
  #declare pos=<-xpos-0.5,0,-ypos>; 
  #declare pos=pos*400;
  pos
#end 
*/

/*                    
#macro comet_position(time)
  #declare e=0.9;
  #declare a=1; 
  #declare GM=2;
  #declare theta=(((1+cos(time*2*pi))/2)^10)*2*pi;
  #declare cos_theta = cos(theta);
  #declare sin_theta = sin(theta);
  #declare rr=a*(1-e*e)/(1+e*cos_theta);
  #declare xpos=rr*cos_theta;
  #declare ypos=rr*sin_theta;
  #declare pos=<xpos+2*e-0.5,0,-ypos>; 
  #declare pos=pos*400;
  pos
#end 
*/

#macro comet_position(time)
  #declare e=0.9;       // eccentricity
  #declare a=400;       // semi-major axis
  #declare T=60;        // period
  #declare M0=-4;       // offset
  
  #declare sun_pos=<-200,0,0>;

  #declare M=-time*2*pi/T+M0;    // mean anomaly

  #declare iteration=0;
  #declare EE = M;       // eccentric anomaly
  #while (iteration < 100)
    #declare EE = M + e * sin(EE);
    #declare iteration = iteration + 1;
  #end
  #declare rr = a * (1 - e * cos(EE));
  
  #declare cos_f = (a * (1 - e*e) - rr)/rr/e;
  #declare f = acos(cos_f);
  #declare sin_f = sin(f);
  #if (mod(-M,2*pi) < pi)
    #declare sin_f = - sin_f;
  #end
  
 //   #render concat (str(cos_f,1,3),"\n")
  
  #declare xpos = -rr*cos_f;
  #declare ypos = rr*sin_f;
  
  #declare pos=<xpos,0,ypos>+sun_pos;
  pos
#end
    
#macro c_speed()       
  #declare speed=old_comet_pos-comet_pos;
  speed
#end

#macro move_comet()  
  #debug "in move_comet()\n"
  #declare time=time+step;
  #declare frame=frame+1;
  #declare old_comet_pos=comet_pos;
  #declare comet_pos=comet_position(time);
  #declare comet_speed=c_speed();
#end

#macro move_particle(num)
  #debug "in move_particle()\n"
  #declare vdelta=particle_position[num]-sun_position;  
  #declare d=vlength(vdelta);
  #declare delta=vnormalize(vdelta);
  //if (delta <> 0)
    //#declare a=1/d;
    //#declare va=delta*a*10;
    //  #declare va=<0,0,0>;                                          
    #declare sun_speed=delta*wind;
    #declare particle_speed[num]=particle_speed[num]*factor+sun_speed*(1-factor);
    #declare particle_position[num]=particle_position[num]+particle_speed[num]+<0,1,0>; 
    #declare particle_age[num]=particle_age[num]+1;
  //#end
#end

#macro move_particle2(num)
  #debug "in move_particle()\n"
  #declare vdelta=particle_position2[num]-sun_position;  
  #declare d=vlength(vdelta);
  #declare delta=vnormalize(vdelta);
  //if (delta <> 0)
    //#declare a=1/d;
    //#declare va=delta*a*10;
    //  #declare va=<0,0,0>;                                          
    #declare sun_speed=delta*wind2;
    #declare particle_speed2[num]=particle_speed2[num]+sun_speed;
    #declare particle_position2[num]=particle_position2[num]+particle_speed2[num]; 
    #declare particle_age2[num]=particle_age2[num]+1;
  //#end
#end

#macro move_all()
  #declare n=0;
  #while (n<max_particles)
    move_particle(n)
    #declare n=n+1;
  #end
  #declare n=0;
  #while (n<max_particles2)
    move_particle2(n)
    #declare n=n+1;
  #end
#end

#macro init()
  #debug "in init()\n"
  #declare n=0;
  #while (n<max_particles)    
    //#debug concat(str(max_particles,1,0),"\n")
    #declare particle_position[n]=<0,0,0>;
    #declare particle_speed[n]=<0,0,0>;
    #declare particle_age[n]=0;
    #declare n=n+1;
  #end
  #declare n=0;
  #while (n<max_particles2)    
    //#debug concat(str(max_particles,1,0),"\n")
    #declare particle_position2[n]=<0,0,0>;
    #declare particle_speed2[n]=<0,0,0>;
    #declare particle_age2[n]=0;
    #declare n=n+1;
  #end  
#end
    
#macro create_particle(pos,speed,num,sv)
  #declare x_speed=-speed.x+(rand(rnd)-0.5)*maxv*2;
  #declare y_speed=-speed.y+(rand(rnd)-0.5)*maxv*2;
  #declare z_speed=-speed.z+(rand(rnd)-0.5)*maxv*2;
  #declare particle_position[num]=pos;
  #declare particle_speed[num]=<x_speed,y_speed,z_speed>+sv*gas*0;
  #declare particle_age[num]=0;
#end

#macro create_particle2(pos,speed,num,sv)
  #declare x_speed=-speed.x+(rand(rnd)-0.5)*maxv2*2;
  #declare y_speed=-speed.y+(rand(rnd)-0.5)*maxv2*2;
  #declare z_speed=-speed.z+(rand(rnd)-0.5)*maxv2*2;
  #declare particle_position2[num]=pos;
  #declare particle_speed2[num]=<x_speed,y_speed,z_speed>+sv*gas*0;
  #declare particle_age2[num]=0;
#end

#macro create_block()
  #debug "in create_block()\n"
  #declare rnd=seed(frame);
  #declare limit=pointer+new_per_frame;
  #declare limit2=pointer2+new_per_frame2;
  #declare sv=vnormalize(sun_position-comet_pos);
  #while (pointer<limit)
    create_particle(comet_pos,comet_speed,pointer,sv)
    #declare pointer=pointer+1;
  #end
  #if (pointer >= max_particles)
    #declare pointer=0;
  #end
  #while (pointer2<limit2)
    create_particle2(comet_pos,comet_speed,pointer2,sv)
    #declare pointer2=pointer2+1;
  #end
  #if (pointer2 >= max_particles2)
    #declare pointer2=0;
  #end
#end 

#macro display_comet()
  #render concat ("comet_position x=",str(comet_pos.x,1,2)," y=",str(comet_pos.y,1,2), " z=",str(comet_pos.z,1,2),"\n") 
  #declare n=0;
  #while (n<max_particles)    
    //#debug concat(str(max_particles,1,0),"\n")
    #declare col=((max_age-particle_age[n])/max_age+0.75);
    #declare dist=vlength(particle_position[n]*<1,0,1>-sun_position);
    #declare pos=particle_position[n]+<0,f_sun(dist),0>+<0,1,0>;
    glow {
      type 0
      location <0,0,0>
      color rgb <.02,.02,.1>
      size 0.05*col*0.5
      radius 3
      fade_power 1 translate pos
      }

    #declare n=n+1;
  #end
  #declare n=0;
  #while (n<max_particles2)    
    //#debug concat(str(max_particles,1,0),"\n")
    #declare col=((max_age2-particle_age2[n])/max_age2+0.75);
    #declare dist=vlength(particle_position2[n]*<1,0,1>-sun_position);
    #declare pos=particle_position2[n]+<0,f_sun(dist),0>+<0,1.5,0>;
    glow {
      type 0
      location <0,0,0>
      color rgb <.2,.2,.2>
      size 0.05*col*0.5
      radius 3
      fade_power 1 translate pos
      }

    #declare n=n+1;
  #end
#end


#if (!dont_display_comet)  
#render "Particle system initializing...\n\n"
#declare comet_pos=comet_position(time);
#declare old_comet_pos=comet_pos;  
init()
#while (time<(initial_time+tmax)) 
  #debug "in loop\n"
  #render concat(str(initial_time+tmax-time,1,5),"   \r")
  #declare rnd=seed(time);
  move_comet()
  create_block()
  move_all()
#end 
#end

//#declare comet_clock=(clock-13)/24;

#macro display_comet_path(sphere_size,init_clock,fin_clock,clock_step,highlight)
  #declare splineclock = init_clock;
  #while (splineclock <= fin_clock) 
    #render concat(str(splineclock,1,3),"\n")
    #declare progress = (splineclock-init_clock)/(fin_clock-init_clock);
    #declare pos = comet_position((splineclock));
    sphere {
      pos,sphere_size
      pigment { color rgb <1,progress,progress> }
      finish { ambient 1 diffuse 0 }
      }
    #declare splineclock = splineclock + clock_step;
    sphere { 
      #declare pos = comet_position((highlight-0)/1);
      pos,sphere_size*3
      pigment { color rgb <0,0,1> }
      finish { ambient 1 diffuse 0 }
      }
    #end
    sphere {
      sun_pos,sphere_size*5
      pigment { color rgb <1,1,0> }
      finish { ambient 1 diffuse 0 }
      }     
  #end

#macro print_cometpos(time)
  #declare pos = comet_position(time);
  #render concat ("comet's position at t=",str(time,1,2)," is <",str(pos.x,1,2),"/",str(pos.y,1,2),"/",str(pos.z,1,2),">\n")
#end





