#version unofficial Megapov 0.7

  #init_spline{"m_large",<0,300>,<13.999,300>,<14,300>,<18,0>,<18.001,0>,<40,0>}
  #init_spline{"m_small",<0,10>,<13.999,10>,<14,10>,<18,0>,<18.001,0>,<40,0>}
  #declare m_large=eval_spline("m_large",clock);
  #declare m_small=eval_spline("m_small",clock);
  #declare r_large=3;
  #declare r_small=1.5;
  #declare s_large=10;
  #declare s_small=2;
  #declare x_large=0-scene_clock;
  #declare z_large=0+scene_clock*20;

  #declare m_sun=8000;
  #declare r_sun=30;
  #declare s_sun=50;
  #declare x_sun=-200;
  #declare z_sun=0;

  #declare rad=12;
  #declare rstep=scene_clock*scene_clock*pi*2+pi;
  #declare x_small=rad*cos(rstep)+x_large;
  #declare z_small=rad*sin(rstep)+z_large;
  #declare a = sqrt (((x_small-x_large)^2)+((z_small-z_large)^2));
  #declare mx=0;  
  
  #declare sun_f = function { m_sun/(sqrt((x-x_sun)*(x-x_sun)+(z-z_sun)*(z-z_sun))+s_sun)  }
  #declare large = function { m_large/(sqrt((x-x_large)*(x-x_large)+(z-z_large)*(z-z_large))+s_large) }
  #declare small = function { m_small/(sqrt((x-x_small)*(x-x_small)+(z-z_small)*(z-z_small))+s_small) }

  #declare universe = function { y+sun_f+large+small }
  
  #macro eval_universe(vv)
    #declare alt=
    m_sun/(sqrt((vv.x-x_sun)*(vv.x-x_sun)+(vv.z-z_sun)*(vv.z-z_sun))+s_sun)
    +m_large/(sqrt((vv.x-x_large)*(vv.x-x_large)+(vv.z-z_large)*(vv.z-z_large))+s_large)
    +m_small/(sqrt((vv.x-x_small)*(vv.x-x_small)+(vv.z-z_small)*(vv.z-z_small))+s_small);
    alt
  #end
    
  #macro f_sun(X)
    -m_sun/(s_sun+X)
  #end
  
  #macro f_large(X)
    -m_large/(s_large+X)
  #end
  
  #macro f_small(X)
    f_large(X)-m_small/(s_small+abs(X-a))
  #end

  #macro f_env(X)
    f_large(X)-m_small/(s_small+X-a)
  #end

  #macro check_large(LY)
    #declare ly=LY;
    #declare psi=0;
    #declare check=0;
    #while(psi<pi/2)
      #declare xx=r_large*sin(psi);
      #declare yy=ly-r_large*cos(psi);
      #if (f_small(xx)>yy)
        #declare check=1;        
      #end 
      #declare psi=psi+pi/720;
      #end
    check
    #end
      
  #macro check_outside(MX,MY)
    #declare mx=MX;
    #declare my=MY;
    #declare psi=0;
    #declare check=0;
    #while(psi<pi/2)
      #declare xx=mx+r_small*sin(psi);
      #declare yy=my-r_small*cos(psi);
      #if (f_env(xx)>yy)
        #declare check=1;
        #end 
      #declare psi=psi+pi/360;
      #end
    check
    #end

  #macro check_inside(MX,MY) 
    #declare mx=MX;
    #declare my=MY;
    #declare psi=0;
    #declare check=0;
    #while(psi<pi/2)
      #declare xx=mx+r_small*sin(-psi);
      #declare yy=my-r_small*cos(-psi); 
      #if (f_small(xx)>yy)
        #declare check=1;
        #end
      #declare psi=psi+pi/360;
      #end
    check
    #end

  // ++++++++++++++++++ Large Mass ++++++++++++++++++++

  #declare xpos=0;
  #declare zpos=0;

  #declare rr=2;
  #declare mm=4;

  #declare ymax=r_large;
  #declare ymin=f_large(0);

  #declare yy = (ymax+ymin)/2;
  #declare step = 0;
 
  #while (step < 20) 
    #declare step = step+1;
    #if (check_large(yy))  
      #declare ymin = yy;
      #declare yy = (yy + ymax) / 2;
    #else  
      #declare ymax = yy;
      #declare yy = (yy + ymin) / 2;
    #end
  #debug concat("Earth y=",str(yy,1,3)," min=",str(ymin,1,3)," max=",str(ymax,1,3),"\n")
  #end
    
  #declare aa=mm;    
  #declare bb=0;
  #declare cc=-mm*(rr^2);

  #declare pp=(3*bb-(aa^2))/3; 
  #declare qq=cc+2*(aa^3)/27-(aa*bb/cc);

  #declare DD=(pp/3)^3+(qq/2)^2;
                               
                               

  #if (DD>=0)
    #declare uu=(-qq/2+sqrt(DD))^(1/3);
    #declare vv=(-qq/2-sqrt(DD))^(1/3);
    #declare xx=sqrt(-aa/3+uu+vv);
  #end

  #if (DD<0)                                   
    #declare cospsi=-qq/(2*sqrt((abs(pp)/3)^3));
    #declare psi=acos(cospsi);
    #declare xx=(-aa/3+2*sqrt(abs(pp)/3)*cos(psi/3));
  #end
                        
  #declare px=xx;
  #declare py=-mm/px; 


  #declare sy=sqrt(rr*rr-px*px);

  object {
    earth
    scale r_large 
    #declare xz = sqrt(x_large*x_large+z_large*z_large)   
    rotate <0,0,-90>
    rotate <0,180,0>
    rotate <degrees(xz/r_large),0,0>
    rotate <0,asin(-1/15),0>
    #declare dist=vlength(<x_large,0,z_large>-<x_sun,0,z_sun>);
    #declare yy=yy+f_sun(dist);
    translate <x_large,yy,z_large>
    #declare earthpos=<x_large,yy,z_large>;    
    }    

  // ++++++++++++++++++ Small Mass ++++++++++++++++++++

  #declare h0 = f_large(a);
  #declare h1 = f_small(a);

  //#render concat("h0 ",str(h0,10,10)," h1 ",str(h1,10,10)," a ",str(a,10,10),"\n")

  #declare mx = a;
  #declare my = h0;

  #while (check_outside(mx,my))
    #declare my = my + r_small;
  #end


  #declare m0 = my;
  #declare my = (m0 + h1) / 2;

  #declare mmax = m0;
  #declare mmin = h1;

  #declare inside = 0;

  #while (!inside)
    #declare step = 0;
    #declare mmax = m0;
    #declare mmin = h1;  
    #while (step < 15) 
      #declare step = step+1;
      #if (check_outside(mx,my))  
        #declare mmin = my;
        #declare my = (my + mmax) / 2;
      #else  
        #declare mmax = my;
        #declare my = (my + mmin) / 2;
      #end
      #debug concat("Moon x=",str(mx,1,3),"y=",str(my,1,3)," min=",str(mmin,1,3)," max=",str(mmax,1,3),"\n")
    #end
    #declare inside = check_inside(mx,my);
    #if (!check_inside(mx,my))
      #declare mx = mx - r_small/50;
    #end
  #end
       
  #declare move=mx/a;

  #declare dy=abs(yy-my)+r_small;

  #declare aa=sqrt(mx*mx+dy*dy);
  #declare phi = -acos(dy/aa);

  motion_blur {
    object {
      moon
      #declare scene_clock=clock/13;                               
      #declare rstep=scene_clock*scene_clock*pi*2+pi;
      scale r_small
      translate <0,aa,0>
      rotate <0,degrees(rstep*aa/2/pi)-90,0>
      rotate <0,0,degrees(phi)>
      rotate <0,-degrees(rstep),0> 
      #declare moonpos=<0,aa,0>;
      #declare moonpos=vrotate(moonpos,<0,degrees(rstep*aa/2/pi)-90,0>);
      #declare moonpos=vrotate(moonpos,<0,0,degrees(phi)>);
      #declare moonpos=vrotate(moonpos,<0,-degrees(rstep),0>);
      #declare moonpos=moonpos+<x_large,yy-r_small,z_large>;
      #declare dist=vlength((moonpos*<1,0,1>)-<x_sun,0,z_sun>);
      #declare yyy=yy+f_sun(dist);
      translate <x_large,yyy-r_small,z_large>
      }
    }