#version unofficial Megapov 0.7

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

  #declare sun_f = function { m_sun/(sqrt((x-x_sun)*(x-x_sun)+(z-z_sun)*(z-z_sun))+s_sun)  }

  #macro f_sun(X)
    -m_sun/(s_sun+X)
  #end
  
  #declare universe = function { y+sun_f }
  
  #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);
    alt
  #end
  
  #macro check_sun(LY)
    #declare ly=LY;
    #declare psi=0;
    #declare check=0;
    #while(psi<pi/2)
      #declare xx=r_sun*sin(psi);
      #declare yy=ly-r_sun*cos(psi);
      #if (f_sun(xx)>yy)
        #declare check=1;        
      #end 
      #declare psi=psi+pi/720;
      #end
    check
    #end
      

  // ++++++++++++++++++ Sun ++++++++++++++++++++

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

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

  #declare ymax=r_sun;
  #declare ymin=f_sun(0);

  #declare yy = (ymax+ymin)/2;
  #declare step = 0;
 
  #while (step < 20) 
    #declare step = step+1;
    #if (check_sun(yy))  
      #declare ymin = yy;
      #declare yy = (yy + ymax) / 2;
    #else  
      #declare ymax = yy;
      #declare yy = (yy + ymin) / 2;
    #end
    #debug concat("Sun 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 {
    sun
    //scale r_sun 
    #declare xz = sqrt(x_sun*x_sun+z_sun*z_sun)   
    rotate <0,0,-90>
    rotate <0,180,0>
    rotate <degrees(xz/r_sun),0,0>
    rotate <0,asin(-1/15),0>
    translate <x_sun,yy,z_sun>  
    }    

