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

#version 3.1;
global_settings {
  max_trace_level 9
  assumed_gamma 2.2
  }

light_source {<500, 1000, -2000> color <0.5,0.5,1>}
light_source {<-6000, -3000, -4000> color <0.54,1,0.5>}

camera {
  location <0,0,-10000>
  angle 58
  look_at <-400, 600, 0>
  }

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

//initialize the coordinates for mesh Cp(x,y,z,r)
#declare Vmax=144;
#declare Hmax=2*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

//radius caliculation and the translation from polar to xyz coordinsates
#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 Hm=H/Hmax;
    #declare Vm=Vvmax+pi/20*sin(2*pi*Vvmax);
    #declare Cp[H][V][3]=Rdg*(1+6*(1+0.6*cos(2*pi*Vm+pi))*pow(1+0.3*cos(16*2*pi*Hm+pi),6)*pow(sin(pi*Vm),16))*(1+0.05*sin(8*pi*Hm));
    #declare Cp[H][V][0]=Cp[H][V][3]*sin(pi*Vm+0.08*Cp[H][V][3]/Rdg)*cos(2*pi*Hm+0.1*Cp[H][V][3]/Rdg);//+0.1*sin(27*Hm+23*Vm));
    #declare Cp[H][V][1]=Cp[H][V][3]*sin(pi*Vm+0.08*Cp[H][V][3]/Rdg)*sin(2*pi*Hm+0.1*Cp[H][V][3]/Rdg);//+0.1*sin(37*Hm+13*Vm));
    #declare Cp[H][V][2]=-Cp[H][V][3]*cos(pi*Vm+0.08*Cp[H][V][3]/Rdg);//+Cp[H][V][3]/Rdg)-0.004*pow(Rdg-Cp[H][V][3],2);
  #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 <1,1,1>}
      finish { reflection 1 ambient -0.20 phong 1 brilliance 1 diffuse 1 }
      }
  rotate <70,40,0>
  no_shadow
  }