#version unofficial MegaPov 0.5;

#include "colors.inc"

global_settings { assumed_gamma 1.0 }

camera
 {
  location  <0.0, 0.0, -2.5>
  direction z
  right     x
  up        (image_height/image_width)*y
  look_at   <0.0, 0.0,  0.0>
 }

sky_sphere
 {
  pigment
   {
    gradient y
    color_map { [0.0 color blue 0.6] [1.0 color rgb 1] }
   }
 }

light_source
 {
  0*x
  color red 1.0  green 1.0  blue 1.0
  translate <-30, 30, -30>
 }

// macro to aproximate elipse

#declare elipse_a = 1;
#declare elipse_b = elipse_a * 0.6 ;

#macro put_ob(obj,dx,dy)
  object { obj translate <dx,dy,0> }
  object { obj translate <-dx,dy,0> }
  object { obj translate <-dx,-dy,0> }
  object { obj translate <dx,-dy,0> }
#end

// along x or y axis - small spheres

#declare step=0.005;
#declare ob=sphere { 0.0, (step/2) texture { pigment { color Gray } } }

#declare dx = 0;
#while ( dx <= 1 )
  #if ( elipse_a >= dx )
    put_ob( ob , dx , ( elipse_b / elipse_a ) * sqrt( ( elipse_a - dx ) * ( elipse_a + dx ) ) )
  #end
  #declare dx=dx+step;
#end

#declare dy = 0;
#while ( dy <= 1 )
  #if ( elipse_b >= dy )
    put_ob( ob , ( elipse_a / elipse_b ) * sqrt( ( elipse_b - dy ) * ( elipse_b + dy ) ) , dy )
  #end
  #declare dy=dy+step;
#end

// selfapproximation with distribution

#declare nseg=4;
#declare step_r = elipse_a / nseg ;
#declare seg=sphere { 0.0, (step_r/3) texture { pigment { color White } finish { reflection .5 } } }
#declare max_error = 1;
#declare table=array[nseg+1][4]

  // 0 - angle
  // 1 - x
  // 2 - y
  // 3 - length

#declare table[0][0] = 0;
#declare table[0][1] = elipse_a;
#declare table[0][2] = 0;
#declare table[0][3] = 0;

#declare dx=1;
#while ( dx<=nseg )
  #declare table[dx][0] = tan(radians(dx*90/nseg));
  #declare dx=dx+1;
#end

#macro calc()
  #declare length=0;
  #declare dx=1;
  #while ( dx<=nseg )
    #declare table[dx][1] = ( 1 / ( elipse_a^2 ) ) + ( ( table[dx][0] / elipse_b )^2 );
    #declare table[dx][1] = sqrt( 1 / table[dx][1] );
    #declare table[dx][2] = table[dx][0] * table[dx][1] ;
    #declare table[dx][3] = ( table[dx][2] - table[dx-1][2] )^2;
    #declare table[dx][3] = table[dx][3] + ( ( table[dx][1] - table[dx-1][1] )^2 );
    #declare table[dx][3] = sqrt( table[dx][3] ) ;
    #declare length = length + table[dx][3] ;
    #declare dx=dx+1;
  #end
#end

#while (max_error>0.000005)
  #declare max_error=0;
  #declare dy=1;
  #while ( dy<nseg )
    calc()
    #declare cur_error=table[dy][3]/(length/nseg);
    #declare table[dy][0] = (table[dy][0]-table[dy-1][0])/cur_error + table[dy-1][0];
    #declare dy=dy+1;
  #end
  calc()
  #warning concat("optimal - ",str(length/nseg,0,30),"\n")
  #declare dy=1;
  #while ( dy<nseg )
    #declare cur_error=abs((table[dy][3]/(length/nseg))-1);
    #if (max_error < cur_error)
      #declare max_error=cur_error;
    #end
    #warning concat(str(dy,7,0)," - ",str(table[dy][3],0,30),"\n")
    #declare dy=dy+1;
  #end
  #warning concat("  error - ",str(max_error,0,30),"\n\n")
#end

#declare dx=0;
#while ( dx<=nseg )
  put_ob( seg , table[dx][1] , table[dx][2] )
  #declare dx=dx+1;
#end