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

#version 3.1;
global_settings {
  max_trace_level 6
  assumed_gamma 2.2
  }

camera {
  location <0,0,-5000>
  angle 0.6
  look_at <0, 0, 0>
  }

//***************** Mesh treatment *****************
//coloring
#macro Col (Ph,Pv)
   texture {
    #switch ( Ph )
    #case( 0 ) pigment { color rgbt <1-Z/Zmax, 0, Z/Zmax, 0> } #break
    #case( 2 ) pigment { color rgbt <0, Z/Zmax, 1-Z/Zmax, 0> } #break
    #case( 1 ) pigment { color rgbt <Z/Zmax, 1-Z/Zmax, 0, 0> } #break
    #end
    finish { reflection 0 ambient 0.2 diffuse 1 brilliance 1 }
    }
#end

//frame of the surface
#macro Trifd (Ph,Pv,R)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
  #declare Pv1=mod(Pv-1+Zmax,Zmax);
  #if (Pv1=0) 
    #declare Pv1=Zmax;
    #declare Phs=mod(Ph+1+Hmax,Hmax);
    #declare Ph1s=mod(Ph1+1+Hmax,Hmax);
  #else
    #declare Phs=Ph;
    #declare Ph1s=Ph1;
  #end
    cylinder {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      <Cp[Ph1s][Pv1][0],Cp[Ph1s][Pv1][1],Cp[Ph1s][Pv1][2]> R
      Col (Ph,Pv)
      }
    cylinder {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]> R
      Col (Ph,Pv)
      }
    cylinder {
      <Cp[Phs][Pv1][0],Cp[Phs][Pv1][1],Cp[Phs][Pv1][2]>
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> R
   texture {
    pigment { color rgb<1, 1, 1> }
    finish { reflection 0 ambient 0.8 diffuse 0.5 brilliance 1 }
    }
      }
#end

//triangle of the surface
#macro Trid (Ph,Pv)
  #declare Green = Col (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
  #declare Pv1=mod(Pv-1+Zmax,Zmax);
  #if (Pv1=0)
    #declare Pv1=Zmax;
    #declare Phs=mod(Ph+1+Hmax,Hmax);
    #declare Ph1s=mod(Ph1+1+Hmax,Hmax);
  #else
    #declare Phs=Ph;
    #declare Ph1s=Ph1;
  #end   
    triangle {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      <Cp[Ph1s][Pv1][0],Cp[Ph1s][Pv1][1],Cp[Ph1s][Pv1][2]>
      <Cp[Phs][Pv1][0],Cp[Phs][Pv1][1],Cp[Phs][Pv1][2]>
      texture { Green}
      }
    triangle {
      <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]>
      <Cp[Ph1s][Pv1][0],Cp[Ph1s][Pv1][1],Cp[Ph1s][Pv1][2]>
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]>
      texture { Green}
      }
#end

//normal vector of the point
#macro Vd (Ph,Pv)
  #declare Phs=mod(Ph+Hmax,Hmax);
  #declare Pvs=mod(Pv+Zmax,Zmax);
  #declare Ph1s=mod(Ph-1+Hmax,Hmax);
  #declare Pv1s=mod(Pv-1+Zmax,Zmax);
  #if (Pvs=0)
    #declare Pvs=Zmax;
  #end
  #if (Pv1s=0)
    #declare Pv1s=Zmax;
  #end
  <(Cp[Ph1s][Pvs][1]-Cp[Phs][Pvs][1])*(Cp[Phs][Pv1s][2]-Cp[Phs][Pvs][2])-(Cp[Phs][Pv1s][1]-Cp[Phs][Pvs][1])*(Cp[Ph1s][Pvs][2]-Cp[Phs][Pvs][2]),
  (Cp[Phs][Pv1s][0]-Cp[Phs][Pvs][0])*(Cp[Ph1s][Pvs][2]-Cp[Phs][Pvs][2])-(Cp[Ph1s][Pvs][0]-Cp[Phs][Pvs][0])*(Cp[Phs][Pv1s][2]-Cp[Phs][Pvs][2]),
  (Cp[Ph1s][Pvs][0]-Cp[Phs][Pvs][0])*(Cp[Phs][Pv1s][1]-Cp[Phs][Pvs][1])-(Cp[Phs][Pv1s][0]-Cp[Phs][Pvs][0])*(Cp[Ph1s][Pvs][1]-Cp[Phs][Pvs][1])>
#end

//smooth_triangle of the surface
#macro Trisd (Ph,Pv)
  #declare Green = Col (Ph,Pv)
  #declare Ph1=mod(Ph-1+Hmax,Hmax);
  #declare Pv1=mod(Pv-1+Zmax,Zmax);
  #if (Pv1=0)
    #declare Pv1=Zmax;
    #declare Pht=mod(Ph+1+Hmax,Hmax);
    #declare Ph1t=mod(Ph1+1+Hmax,Hmax);
  #else
    #declare Pht=Ph;
    #declare Ph1t=Ph1;
  #end   
    smooth_triangle {
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd(Ph,Pv)
      <Cp[Ph1t][Pv1][0],Cp[Ph1t][Pv1][1],Cp[Ph1t][Pv1][2]> Vd(Ph1t,Pv1)
      <Cp[Pht][Pv1][0],Cp[Pht][Pv1][1],Cp[Pht][Pv1][2]> Vd(Pht,Pv1)
      texture { Green}
      }
    smooth_triangle {
      <Cp[Ph1][Pv][0],Cp[Ph1][Pv][1],Cp[Ph1][Pv][2]> Vd(Ph1,Pv)
      <Cp[Ph1t][Pv1][0],Cp[Ph1t][Pv1][1],Cp[Ph1t][Pv1][2]> Vd(Ph1t,Pv1)
      <Cp[Ph][Pv][0],Cp[Ph][Pv][1],Cp[Ph][Pv][2]> Vd (Ph,Pv)
      texture { Green}
      }
#end


//***************** Hose coordinates *****************
#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 Ring (Atwist, Aallx,Aally,Aallz, Pallx,Pally,Pallz, Vallys,Vallyc,Vallyd,Vallz)
#declare Z=0; #while (Z<Zmax+2)
    #declare Cc[Z][2]=Rdg*cos(2*pi*(Z-1)/Zmax);
    #declare Cc[Z][1]=0;
    #declare Cc[Z][0]=Rdg*sin(2*pi*(Z-1)/Zmax);
    #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*abs(Z-2)/Zmax+0.000001;
      Vrect (Z)
    #end
  #declare H=0; #while (H<Hmax)
    #declare Arot=Z/Zmax*Atwist; //twist
    #declare Hhmax=H/Hmax-Arot/360;
    #declare Cp[H][Z][3]=Rdgg;
    #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

#switch ( Id )
  #case( 1 )
  union {
    #declare Z=1; #while (Z<Zmax+1)
    #declare H=0; #while (H<Hmax)
      Trifd (H,Z,0.1) 
    #declare H=H+1; #end
    #declare Z=Z+1; #end
    }
  #break
  #case( 2 ) 
  mesh {
    #declare Z=1; #while (Z<Zmax+1)
    #declare H=0; #while (H<Hmax)
      Trid (H,Z) 
    #declare H=H+1; #end
    #declare Z=Z+1; #end
    }
  #break
  #case( 3 )
  mesh {
    #declare Z=1; #while (Z<Zmax+1)
    #declare H=0; #while (H<Hmax)
      Trisd (H,Z) #break
    #declare H=H+1; #end
    #declare Z=Z+1; #end
    }
  #break
#end

#end

#declare Id=1; #while (Id<4)
#declare Rdg=10;  //radius of circle
#declare Rdgg=3;  //radius of ring

//initialize the coordinates for mesh Cp(x,y,z,r)
#declare Zmax=24;
#declare Hmax=3;
#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

union {
   Ring (480,0,0,0, 0,0,0, 4,1.5,pi/2,1.5)  
   rotate <-45,15,0>
   translate <13*(Id-2)-3+3*abs(Id-2),-9*pow(-1,Id),0>
   no_shadow
   }
#declare Id=Id+1; #end
   
light_source {<0, 0, -5000> color <1,1,1>}

background { color rgb <0,0,0.1> }