// HOSE / Dolphin '00.03
// 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.8,0.8,1>}

camera {
  location <0,0,-100>
  angle 56
  look_at <0, 15, 0>
  }

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


//initialize the coordinates for mesh Cp(x,y,z,r)
#declare Zmax=60;
#declare Hmax=36;
#declare Cp=array[Hmax][Zmax+2][4]
#declare Cc=array[Zmax+2][3]
#declare Z=0; #while (Z<Zmax+1)
  #declare I=0; #while (I<3)
    #declare Cc[Z][I]=0;
  #declare I=I+1; #end
#declare H=0; #while (H<Hmax)
  #declare I=0; #while (I<4)
    #declare Cp[H][Z][I]=0;
  #declare I=I+1; #end
#declare H=H+1; #end
#declare Z=Z+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

#macro Vrect (P)
  #if (Cc[P][2]-Cc[P-1][2]=0)
    #declare Cc[P-1][2]=Cc[P-1][2]-0.000001;
  #end
  #declare Vvx=Cc[P][0]-Cc[P-1][0];
  #declare Vvy=Cc[P][1]-Cc[P-1][1];
  #declare Vvz=Cc[P][2]-Cc[P-1][2];
  #declare Vvl=sqrt(Vvx*Vvx+Vvy*Vvy+Vvz*Vvz)+0.000001;
  #declare Vvx=Vvx/Vvl;
  #declare Vvy=Vvy/Vvl;
  #declare Vvz=Vvz/Vvl;
  #declare Vxz=(Cc[P-1][0]-Cc[P][0])/(Cc[P][2]-Cc[P-1][2]);
  #declare Vxl=sqrt(1+Vxz*Vxz);
  #declare Vxx=1/Vxl;
  #declare Vxy=0;
  #declare Vxz=Vxz/Vxl;
  #declare Vyz=-(Cc[P][1]-Cc[P-1][1])*(Cc[P][2]-Cc[P-1][2])/(pow(Cc[P][0]-Cc[P-1][0],2)+pow(Cc[P][2]-Cc[P-1][2],2));
  #declare Vyx=Vyz*(Cc[P][0]-Cc[P-1][0])/(Cc[P][2]-Cc[P-1][2]);
  #declare Vyl=sqrt(Vyx*Vyx+1+Vyz*Vyz);
  #declare Vyx=Vyx/Vyl;
  #declare Vyy=1/Vyl;
  #declare Vyz=Vyz/Vyl;
  #if (Cc[P][2]-Cc[P-1][2]<0)
    #declare Vxx=-Vxx;
    #declare Vxz=-Vxz;
  #end      
  #if (Cc[P][2]-Cc[P-1][2]=0)
    #declare Cc[P-1][2]=Cc[P-1][2]+0.000001;
  #end
#end

#macro Dolphin(Atwist, Aallx,Aally,Aallz, Pallx,Pally,Pallz, Vallys,Vallyc,Vallyd,Vallz)
#declare Z=0; #while (Z<Zmax+2)
      #declare Prx=0.03;   //
    #declare Cc[Z][0]=Rdg*10*(Z-1)/Zmax-5*Rdg;
    #declare Cc[Z][1]=Rdg*Vallys*sin(Z/Zmax*Vallyc*pi+Vallyd)-Rdg*Vallys*sin(Vallyd);
    #declare Cc[Z][2]=Rdg*Vallz*pow(Z/Zmax,4);
    #if (Z=0)
      #declare Zzmax=0;
      #declare Vvx=1;
      #declare Vvy=0;
      #declare Vvz=0;
      #declare Vxx=1;
      #declare Vxy=0;
      #declare Vxz=0;
      #declare Vyx=0;
      #declare Vyy=1;
      #declare Vyz=0;
   #else
      #declare Zzmax=0.999998*(Z-1)/Zmax+0.000001;
      Vrect (Z)
    #end
  #declare H=0; #while (H<Hmax)
    #declare Arot=Z/Zmax*Atwist; //twist
    #declare Hhmax=H/Hmax-Arot/360;
    //body
    #declare Pa=sin(0.95*pi*(Z-1)/Zmax+pi*0.1);
      #declare Pb=0.9;
      #declare Pc=1-Z/Zmax;
      #declare Pd=0.6;
    #if (Pa<0.16) #declare Pa=0.0; #end
    #if (Z/Zmax<0.01) #declare Pd=0.0; #end
    //dorsal
    #if (abs(Z/Zmax-0.47)<0.10)
      #if (abs(H/Hmax-1/4)<1/3/Hmax+0.00001) #declare Pse=0.5+0.5*sin(2*pi*H/Hmax); #else #declare Pse=0; #end
      #declare Ped=300000000000*Pse*pow(Z/Zmax-0.37,8)*pow(Z/Zmax-0.57,4);
    #else  
      #declare Ped=0;
    #end
    #if (abs(Z/Zmax-0.45)<0.15)
      #if (abs(H/Hmax-1/4)<1/Hmax+0.00001) #declare Pse=0.7+0.7*sin(2*pi*H/Hmax); #else #declare Pse=0; #end
      #declare Pe=30000000*Pse*pow(Z/Zmax-0.30,6)*pow(Z/Zmax-0.60,4);
    #else  
      #declare Pe=0;
    #end
      #declare Pe=Pe+Ped;
    //head
    #if (Z/Zmax<0.15)
        #declare Phs=pow(Z/Zmax-0.1,0.25)/4;
    #else
        #declare Phs=pow(0.05,0.25)/4;
    #end
    #if (Z/Zmax>0.1) 
      #declare Ahd=pow((Z-Zmax*0.1)/(Z+Zmax-Zmax*0.1),0.05)*pi/2;
      #declare Ah=pi/2-(Z/Zmax-0.1)*pi;
      #if (abs(Ah-2*pi*H/Hmax)<Ahd)
        #declare Phl=100*Phs*(1-Z/Zmax)/Zmax*pow(cos(pi/2/Ahd*(2*pi*H/Hmax-Ah)),6);
      #else
        #if (abs(Ah+2*pi-2*pi*H/Hmax)<Ahd)
          #declare Phl=100*Phs*(1-Z/Zmax)/Zmax*pow(cos(pi/2/Ahd*(2*pi*H/Hmax-Ah)),6);
        #else
          #declare Phl=0;
        #end
      #end
      #declare Ah=pi/2+(Z/Zmax-0.1)*pi;
      #if (abs(Ah-2*pi*H/Hmax)<Ahd)
        #declare Phr=100*Phs*(1-Z/Zmax)/Zmax*pow(cos(pi/2/Ahd*(2*pi*H/Hmax-Ah)),6);
      #else
        #declare Phr=0;
      #end
    #else
      #declare Phl=0; #declare Phr=0;
    #end
      #declare Ph=Phr+Phl;
    //tail
    #if (abs(Z/Zmax-0.9)<0.1)
      #declare Pso=pow((0.5+0.5*cos(4*pi*H/Hmax)),8)*Z/(1.2-cos(4*pi*H/Hmax));
    #else
      #declare Pso=0;
    #end
        #declare Po=1200*Pso*pow(Z/Zmax-0.8,5)*(Z/Zmax-0.78)*(1-Z/Zmax)*(5-Z/Zmax);  
    //pectoral
    #if (abs(Z/Zmax-0.30)<0.04)
      #if (abs(sin(2*pi*H/Hmax)+0.6)<0.15) #declare Psm=0.5+0.5*cos(2*pi*H/Hmax*12);
      #else
        #declare Psm=0;
      #end
    #else
      #declare Psm=0;
    #end
        #declare Pm=400000*Psm*pow(Z/Zmax-0.26,3)*(0.34-Z/Zmax)*(1-Z/Zmax);
    //
    #if (Z<2)
      #declare Cp[H][Z][3]=0.01*Rdg;  
    #else
      #declare Cp[H][Z][3]=Rdg*(Pe+Po+Pm+Ph)+Pd*Rdg*(Pb+0.16*(
                  0.8+Z/Zmax*cos(4*pi*H/Hmax-pi))) //ratio X-Y
                  *(0.15*(Zmax-Z+0.2*Zmax)/Zmax
                  +pow(Pa,1.2*Zmax/(0.3*Zmax+(Z-1))) //ratio along Z
                  *(1.1+Pc*sin(2*pi*H/Hmax) //biase
                  ));
    #end
    #declare Cp[H][Z][0]=Prx*Vvx*Cp[H][Z][3]*Cp[H][Z][3]+Cc[Z][0]+Cp[H][Z][3]*(Vxx*cos(2*pi*Hhmax)+Vyx*sin(2*pi*Hhmax));
    #declare Cp[H][Z][1]=Prx*Vvy*Cp[H][Z][3]*Cp[H][Z][3]+Cc[Z][1]+Cp[H][Z][3]*(Vxy*cos(2*pi*Hhmax)+Vyy*sin(2*pi*Hhmax));
    #declare Cp[H][Z][2]=Prx*Vvz*Cp[H][Z][3]*Cp[H][Z][3]+Cc[Z][2]+Cp[H][Z][3]*(Vxz*cos(2*pi*Hhmax)+Vyz*sin(2*pi*Hhmax));
  #declare H=H+1; #end
#declare Z=Z+1; #end

mesh {
    #declare Z=2; #while (Z<Zmax+1)
    #declare H=0; #while (H<Hmax)
      Trisd (H,Z)
    #declare H=H+1; #end
    #declare Z=Z+1; #end
  texture {
    pigment { color rgb <0.6,0.6,0.6>}
    finish { reflection 1 ambient 0.3 phong 1 brilliance 1 diffuse 1 }
    normal { bumps 0.5 scale 20 }
    }
  rotate <Aallx,Aally,Aallz>
  translate <Pallx,Pally,Pallz>
  }
#end

//Dolphin(Atwist, Aallx,Aally,Aallz, Pallx,Pally,Pallz, Vallys,Vallyc,Vallyd,Vallz)
Dolphin (-20,-30,150,0, 5,0,0, 0.5,1.2,pi/2,1)
Dolphin (-40,15,230,24, 33,28,50, 5,0.5,0,-1)
Dolphin (20,-25,180,5, 10,23,100, 1.8,0,-pi*0.3,0)

//background {color rgbt <0.0,0.8,0.5>}

fog {
  fog_type 1
  distance 200
  color rgb <0,0.5,0.8>
  }