// HOSE / Mirror '00.03
// By Tsutomu HIGO  E-mail: nj2t-hg@asahi-net.or.jp

#version 3.1;
global_settings {
  assumed_gamma 2.2
  }
  
light_source {<2, 5, -5> color <0.5,0.5,0.4>}
light_source {<-5, 5, 0> color <0.5,0.5,0.5>}
 
camera {
  location <-5.0,-1.0,-12.0>
  angle 50
  look_at <1.2, 0, 0>
  }
   
//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 Woman(Atwistn,Atwistw, Aallx,Aally,Aallz, Pallx,Pally,Pallz, Colr,Colg,Colb)
#fopen MyFile "lphsmirr.prn" read
#read (MyFile, Hmax,Zmax,Rdg)
//initialize the coordinates for mesh Cp(x,y,z,r)
#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

#declare Z=0; #while (Z<Zmax+2)
      #declare Prx=0.03;
    #declare Cc[Z][0]=Rdg*8*(Z-1)/Zmax-4*Rdg;
    #declare Cc[Z][1]=Rdg*0;
    #declare Cc[Z][2]=Rdg*0;
    #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)
    #if (Z<6) #declare Arot=0; 
    #else
      #if (Z<8) #declare Arot=(Z-5)/2*Atwistn;
      #else
        #if (Z<12) #declare Arot=Atwistn;
        #else #declare Arot=Atwistn+(Z-11)/(Zmax-11)*Atwistw;
        #end
      #end 
    #end//twist
    #declare Hhmax=H/Hmax-Arot/360;
    #read (MyFile,D,D,Cp[H][Z][3])
    #declare Cp[H][Z][0]=Cc[Z][0]+Cp[H][Z][3]*(Vxx*cos(2*pi*Hhmax)+Vyx*sin(2*pi*Hhmax));
    #declare Cp[H][Z][1]=Cc[Z][1]+Cp[H][Z][3]*(Vxy*cos(2*pi*Hhmax)+Vyy*sin(2*pi*Hhmax));
    #declare Cp[H][Z][2]=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
#fclose MyFile
 
mesh {
    #declare Z=2; #while (Z<Zmax+2)
    #declare H=0; #while (H<Hmax)
      Trisd (H,Z)
    #declare H=H+1; #end
    #declare Z=Z+1; #end
  texture {
    pigment { color rgb <Colr,Colg,Colb>}
    finish { reflection 0 ambient 0.4 phong 0.2 brilliance 0.5 diffuse 0.5 }
    }
  rotate z*Aallz rotate y*Aally rotate x*Aallx
  translate <Pallx,Pally,Pallz>
  no_shadow
  }
    
#fclose MyFile

#end

union {
//Woman(Atwistn,Atwistw, Aallx,Aally,Aallz, Pallx,Pally,Pallz, Colr,Colg,Colb)
  Woman (-20,2,0,0,-90, -20,0,0, 1,0.99,0.97)
  scale 0.1
  }

#declare Mirror=
  merge {
    box { <0, -4, -2>  < 0.001,  2.5,  2> }
    cylinder { <0, 2.5, 0>  < 0.001,  2.5,  0> 2}
    texture {
      pigment { color rgb <0.8,0.9,1>}
      finish { reflection 0.9 ambient -0.1 brilliance 0.5}
      }
    }

union {
  merge { Mirror }
  intersection {
    merge { Mirror translate z*-2}
    box { <-1, -5, 3>  < 1,  5,  -2> }    
    rotate y*19.5
    translate z*-2
    }      
  intersection {
    merge { Mirror translate z*2}
    box { <-1, -5, -3>  < 1,  5,  2> }    
    rotate y*-25
    translate z*2
    }  
  rotate y*-30 translate <3,0,0>
  }