// POLAR / Rotation '00.02
// By Tsutomu HIGO  E-mail: nj2t-hg@asahi-net.or.jp

#version 3.1;
global_settings {
  max_trace_level 3
  assumed_gamma 2.2
  }

light_source {<20000, 100000, -50000> color <0.5,0.5,1>}
light_source {<-60000, -30000, -40000> color <0.5,1,0.5>}

camera {
  location <0,0,-10000>
  angle 3.2
  look_at <-5, 5, 0>
  }

#declare Rdg=100;  //radius of arranging sphere

//initialize the coordinates for mesh Cp(x,y,z,r)
#declare Vmax=72;
#declare Hmax=4*Vmax;
#declare Cp=array[Hmax][Vmax+2][4]
#declare V=0; #while (V<Vmax+1)
#declare H=0; #while (H<Hmax)
#declare I=0; #while (I<4)
  #declare Cp[H][V][I]=0;
#declare I=I+1; #end
#declare H=H+1; #end
#declare V=V+1; #end

//triangle of the surface
#macro Trid (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
    triangle {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]>
      <Cp[Ph][Pv-1][0],Cp[Ph][Pv-1][1],Cp[Ph][Pv-1][2]>
      }
    triangle {
      <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]>
      <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]>
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      }
#end

//normal vector of the plar coordinates point
#macro Vd (Ph,Pv)
  <(Cp[mod(Ph-1+Hmax,Hmax)][Pv][1]-Cp[Ph][Pv][1])*(Cp[Ph][Pv-1][2]-Cp[Ph][Pv][2])-(Cp[Ph][Pv-1][1]-Cp[Ph][Pv][1])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][2]-Cp[Ph][Pv][2]),
  (Cp[Ph][Pv-1][0]-Cp[Ph][Pv][0])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][2]-Cp[Ph][Pv][2])-(Cp[mod(Ph-1+Hmax,Hmax)][Pv][0]-Cp[Ph][Pv][0])*(Cp[Ph][Pv-1][2]-Cp[Ph][Pv][2]),
  (Cp[mod(Ph-1+Hmax,Hmax)][Pv][0]-Cp[Ph][Pv][0])*(Cp[Ph][Pv-1][1]-Cp[Ph][Pv][1])-(Cp[Ph][Pv-1][0]-Cp[Ph][Pv][0])*(Cp[mod(Ph-1+Hmax,Hmax)][Pv][1]-Cp[Ph][Pv][1])>
#end

//smooth_triangle of the surface
#macro Trisd (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
    smooth_triangle {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd(Ph,Pv)
      <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]> Vd(Ph1,Pv-1)
      <Cp[Ph][Pv-1][0],Cp[Ph][Pv-1][1],Cp[Ph][Pv-1][2]> Vd(Ph,Pv-1)
      }
    smooth_triangle {
      <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]> Vd(Ph1,Pv)
      <Cp[Ph1][Pv-1][0],Cp[Ph1][Pv-1][1],Cp[Ph1][Pv-1][2]> Vd(Ph1,Pv-1)
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd (Ph,Pv)
      }
#end

#declare V=0; #while (V<Vmax+2)
  #declare H=0; #while (H<Hmax)
    #if (V=0) #declare Vvmax=0; #else #declare Vvmax=0.999998*(V-1)/Vmax+0.000001; #end
    #declare Cp[H][V][3]=2.8*Rdg*(0.01+pow(0.1+pow(sin(7*pi*Vvmax)*sin(24*pi*H/Hmax),6)*pow(sin(pi*Vvmax),50),0.2));
    #declare Cp[H][V][0]=Cp[H][V][3]*sin(pi*Vvmax)*cos(2*pi*H/Hmax+pow(Cp[H][V][3]*sin(pi*Vvmax)*0.004,2));//+Cp[H][V][3]*sin(pi*Vvmax)*0.002);
    #declare Cp[H][V][1]=Cp[H][V][3]*sin(pi*Vvmax)*sin(2*pi*H/Hmax+pow(Cp[H][V][3]*sin(pi*Vvmax)*0.004,2));//+Cp[H][V][3]*sin(pi*Vvmax)*0.002);
    #declare Cp[H][V][2]=Cp[H][V][3]*cos(pi*Vvmax);
  #declare H=H+1; #end
#declare V=V+1; #end

//main
mesh {
  #declare V=2; #while (V<Vmax+2)
    #declare H=0; #while (H<Hmax)
      Trisd (H,V)
    #declare H=H+1; #end
  #declare V=V+1; #end
  texture {
    pigment { color rgb <0,0.5,1>}
    finish { reflection 2 ambient -0.2 phong 1 brilliance 1 diffuse 1 }
    normal { bumps 1 scale 20 }
    }
  scale z*2/3
  rotate <60,320,0>
  no_shadow
  }

background {color rgbt <0,0,0.2>}