#include "colors.inc"
                           
#macro MatRho4() #end

#declare SCAL = 25;    // 80: a few; 10: many;           
#declare EndScale = 3;

#macro RhombusRho(P1,P2,P3,P4)
  mesh {
    triangle {P1, P2, P3}
    triangle {P2, P4, P3}
  }
#end


#declare c36 = cos(radians(36));
#declare s36 = sin(radians(36));
#declare c72 = cos(radians(72));
#declare s72 = sin(radians(72));

#declare height = 1;

#declare P11 = <     0,  0*height,    0>;
#declare P12 = <   c72, .5*height, -s72>;
#declare P13 = <   c72, .5*height,  s72>;
#declare P14 = < 2*c72,  1*height,    0>;

#macro ThinUpRho()  
  object { RhombusRho(P11,P12,P13,P14) }
#end

#declare P21 = <     0, 1.5*height,    0>;
#declare P22 = <   c72,   1*height, -s72>;
#declare P23 = <   c72,   1*height,  s72>;
#declare P24 = < 2*c72,  .5*height,    0>;

#macro ThinDownRho() 
  object { RhombusRho(P21,P22,P23,P24) }   
#end

#declare P31 = <     0,  0*height,    0>;
#declare P32 = <   c36, .5*height, -s36>;
#declare P33 = <   c36, .5*height,  s36>;
#declare P34 = < 2*c36,  1*height,    0>;

#macro ThickUpRho() 
  object { RhombusRho(P31,P32,P33,P34) }  
#end

#declare P41 = <     0, 1.5*height,    0>;
#declare P42 = <   c36,   1*height, -s36>;
#declare P43 = <   c36,   1*height,  s36>;
#declare P44 = < 2*c36,  .5*height,    0>;

#macro ThickDownRho() 
  object { RhombusRho(P41,P42,P43,P44) }  
#end

#declare PAGE_WIDTH  = 400;
#declare PAGE_HEIGHT = 300;  

#declare Debug = 1;
#declare re = 0;
#declare im = 1;

#macro Roof(R,RT)
  #declare RT = int(R);
  #if (RT < R) #declare RT = RT+1; #end
  //#if (Debug > 3) #write (DEBUG , concat("        // Roof(",str(R,0,5),") = ",str(RT,5,0),"\n")) #end
#end

#macro Floor(R,RT) 
  #declare RT = int(R);
  #if (RT > R) #declare RT = RT-1; #end
  //#if (Debug > 3) #write (DEBUG , concat("        // Floor(",str(R,0,5),") = ",str(RT,5,0),"\n")) #end
#end

#macro Round(R,RT)
  Floor(R+.5,RT)
  //#if (Debug > 3) #write (DEBUG , concat("        // Round(",str(R,0,5),") = ",str(RT,5,0),"\n")) #end
#end

#macro Delta(I, J)
  (#if (I = J)  1  #else  0  #end)
#end

#macro ZetaProd(A, R)
  #declare R[re] = A[0]*Zeta[0][re] + A[1]*Zeta[1][re] + A[2]*Zeta[2][re] + A[3]*Zeta[3][re] + A[4]*Zeta[4][re];
  #declare R[im] = A[0]*Zeta[0][im] + A[1]*Zeta[1][im] + A[2]*Zeta[2][im] + A[3]*Zeta[3][im] + A[4]*Zeta[4][im];
  //#if (Debug > 2) #write (DEBUG , concat("        // ZetaProd({",str(A[0],5,0),",",str(A[1],5,0),",",str(A[2],5,0),",",str(A[3],5,0),",",str(A[4],5,0),"}) = (",str(R[re],0,5),",",str(R[im],0,5),")\n")) #end
  
#end

#macro F(Z, R)         
  #local K = array[5] {0,0,0,0,0}         
  #local I = 0;
  #while (I < 5)
    Roof(Z[re] * Zeta[5-I][re] - Z[im] * Zeta[5-I][im] + Gamma[I], K[I])
    #local I = I+1;
  #end
  ZetaProd(K,R)
  //#if (Debug > 2) #write (DEBUG , concat("        // F((",str(Z[re],0,5),",",str(Z[im],0,5),")) = (",str(R[re],0,5),",",str(R[im],0,5),")\n")) #end
#end

#macro Inside(X,Y)
  (X <= Width+2 & Y <= Height+2 & X >= -Width-2 & Y >= -Height-2)
#end

#macro AddRhombus(Shape, ShapeStr, GR, RD)
  #local Dx = RD[re] - GR[re];
  #local Dy = RD[im] - GR[im];
  #local Angle = 0;
  #if (Dx = 0)
    #local Angle = 90;
  #else         
    #if (Dx > 0)
      #if (Dy >= 0)
        Round(atan2(Dy, Dx) * 180 / pi, Angle)
      #else                                              
        Round(- atan2(-Dy, Dx) * 180 / pi, Angle)
      #end
    #else
      #if (Dy >= 0)
        Round(180 - atan2(Dy, -Dx) * 180 / pi, Angle)
      #else
        Round(180 + atan2(-Dy, -Dx) * 180 / pi, Angle)
      #end
    #end
  #end

  #if (Angle <= 0) #local Angle = Angle + 360; #end

  #if (Inside(GR[re], GR[im]))
    object{ Shape  rotate y*(360-Angle)   translate <GR[re],0,GR[im]> }
    //#if (Debug > 0) #write (DEBUG , concat("        // object{ ", ShapeStr, "() rotate y*",str(360-Angle,5,0),"   translate <",str(GR[re],0,5),",0,",str(GR[im],0,5),"> }\n")) #end
  #end
#end

#macro KRExtremes(R, KRmin,KRmax)
  #local b1 = GridMin[re] * Zeta[5-R][re] - GridMin[im] * Zeta[5-R][im] + Gamma[R];
  #local b2 = GridMax[re] * Zeta[5-R][re] - GridMin[im] * Zeta[5-R][im] + Gamma[R];
  #local b3 = GridMin[re] * Zeta[5-R][re] - GridMax[im] * Zeta[5-R][im] + Gamma[R];
  #local b4 = GridMax[re] * Zeta[5-R][re] - GridMax[im] * Zeta[5-R][im] + Gamma[R];
  #declare KRmin = 0; Roof(min(min(b1,b2),min(b3,b4)), KRmin)
  #declare KRmax = 0; Floor(max(max(b1,b2),max(b3,b4)), KRmax)
  //#if (Debug > 1) #write (DEBUG , concat("      krmin=",str(KRmin,5,0)," krmax=",str(KRmax,5,0),"\n")) #end
#end

#macro KSExtremes(R,KR,S,N, KSmin,KSmax)
  #if (R = 0)
      #local b1 = -100000000000000.0;
      #local b2 =  100000000000000.0;
      #local b3 = ((KR-Gamma[R])*Zeta[5-S][re]-N*GridMin[im]) + Gamma[S];
      #local b4 = ((KR-Gamma[R])*Zeta[5-S][re]-N*GridMax[im]) + Gamma[S];
  #else
      #local b1 = ((KR-Gamma[R])*Zeta[5-S][im]-N*GridMin[re])/Zeta[5-R][im] + Gamma[S];
      #local b2 = ((KR-Gamma[R])*Zeta[5-S][im]-N*GridMax[re])/Zeta[5-R][im] + Gamma[S];
      #local b3 = ((KR-Gamma[R])*Zeta[5-S][re]-N*GridMin[im])/Zeta[5-R][re] + Gamma[S];
      #local b4 = ((KR-Gamma[R])*Zeta[5-S][re]-N*GridMax[im])/Zeta[5-R][re] + Gamma[S];
  #end

  #if (b1 > b2)
      #local b  = b1;
      #local b1 = b2;
      #local b2 = b;
  #end
  #if (b3 > b4)
      #local b  = b3;
      #local b3 = b4;
      #local b4 = b;
  #end

  #declare KSmin = 0; Floor(max(b1,b3),KSmin)
  #declare KSmax = 0; Roof(min(b2,b4), KSmax)
  //#if (Debug > 1) #write (DEBUG , concat("      ksmin=",str(KSmin,5,0)," ksmax=",str(KSmax,5,0),"\n")) #end
#end


#local H = 2 * pi / 5;

#declare Zeta = array[6][2]

#declare Zeta[0][re] = 1;
#declare Zeta[0][im] = 0;

#declare Zeta[1][re] = cos(H);
#declare Zeta[1][im] = sin(H);

#declare Zeta[2][re] = cos(2*H);
#declare Zeta[2][im] = sin(2*H);

#declare Zeta[3][re] = cos(2*H);
#declare Zeta[3][im] = -sin(2*H);

#declare Zeta[4][re] = cos(H);
#declare Zeta[4][im] = -sin(H);

#declare Zeta[5][re] = 1;
#declare Zeta[5][im] = 0;

#declare Gamma = array[5]
#declare Gamma[0] = .1;
#declare Gamma[1] = .1;
#declare Gamma[2] = .1;
#declare Gamma[3] = .1;
#declare Gamma[4] = - Gamma[0] - Gamma[1] - Gamma[2] - Gamma[3];

#declare Eta = array[2]
#declare Eta[re] = Gamma[0]*Zeta[0][re] + Gamma[1]*Zeta[1][re] + Gamma[2]*Zeta[2][re] + Gamma[3]*Zeta[3][re] + Gamma[4]*Zeta[4][re];
#declare Eta[im] = Gamma[0]*Zeta[0][im] + Gamma[1]*Zeta[1][im] + Gamma[2]*Zeta[2][im] + Gamma[3]*Zeta[3][im] + Gamma[4]*Zeta[4][im];
      
      
      

#declare Width  = PAGE_WIDTH / (2 * SCAL);
#declare Height = PAGE_HEIGHT / (2 * SCAL);

#local C = 1 + 2 * cos(pi/5) + cos(2*pi/5);

#declare GridMin = array[2]
#declare GridMax = array[2]
#declare GridMin[re] = (-Width  - Eta[re] - C)*2/5;
#declare GridMax[re] = ( Width  - Eta[re] + C)*2/5;
#declare GridMin[im] = (-Height - Eta[im] - C)*2/5;
#declare GridMax[im] = ( Height - Eta[im] + C)*2/5;

#declare MIN = array[2] {0,0}
#declare MAX = array[2] {0,0}
F(GridMin,MIN)
F(GridMax,MAX)



#macro Penta()
    #declare R = 0; 
    #while (R<5)  /* grid1 */
      //#if (Debug > 1) #write (DEBUG , concat("r = ",str(R,5,0), "\n")) #end

      #local KRmin = 0; #local KRmax = 0; KRExtremes(R, KRmin,KRmax)
      #local KR=KRmin;  
      #while (KR<=KRmax) // line of grid1
        //#if (Debug > 1) #write (DEBUG , concat("  kr = ",str(KR,5,0), "\n")) #end
	#local S = R+1;
	#while (S<5) // grid2
          //#if (Debug > 1) #write (DEBUG , concat("    s = ",str(S,5,0), "\n")) #end

	  #local N = Zeta[5-R][re] * Zeta[5-S][im] - Zeta[5-S][re] * Zeta[5-R][im];
	  #local KSmin = 0; #local KSmax = 0; KSExtremes(R,KR,S,N, KSmin,KSmax)
	  #local KS=KSmin;
	  #while (KS<=KSmax) // line of grid2
            //#if (Debug > 1) #write (DEBUG , concat("      ks = ",str(KS,5,0), "\n")) #end
            #local Z = array[2]
	    #local Z[re] = ( (KR - Gamma[R]) * Zeta[5-S][im] - (KS - Gamma[S]) * Zeta[5-R][im] ) / N;
	    #local Z[im] = ( (KR - Gamma[R]) * Zeta[5-S][re] - (KS - Gamma[S]) * Zeta[5-R][re] ) / N;

            #local KL = array[5] {0,0,0,0,0}
            
            #local I = 0;
            #while (I < 5)
              #if (R = I)
                #local KL[I] = KR;
	      #else
	        #if (S = I)
	          #local KL[I] = KS;
	        #else
	          Roof(Z[re] * Zeta[5-I][re] - Z[im] * Zeta[5-I][im] + Gamma[I], KL[I])
                #end //#if (S = I)
	      #end //#if (R = I)
              //#if (Debug > 2) #write (DEBUG , concat("          KLI = ",str(KL[I],5,0), "\n")) #end
              #local I = I+1;
            #end

	    #local FL = array[2] {0,0}
	    ZetaProd(KL, FL)
	    #local IndL = (KL[0] + KL[1] + KL[2] + KL[3] + KL[4]);
            //#if (Debug > 1) #write (DEBUG , concat("      indl = ",str(IndL,5,0), "\n")) #end

            #local KH = array[5] {0,0,0,0,0}
	    #local I = 0;
	    #while (I < 5)
 	      #local KH[I] = KL[I] + Delta(I,R) + Delta(I,S);
	      #local I = I+1;                        
	    #end

	    #local FH = array[2] {0,0}
	    ZetaProd(KH, FH)
	    #local IndH = (KH[0] + KH[1] + KH[2] + KH[3] + KH[4]);
            
            #local GR = array[2] {0,0}
            #local RD = array[2] {0,0}
	    #if (IndL = 1)
	      #local GR[re] = FL[re]; #local GR[im] = FL[im];
	      #local RD[re] = FH[re]; #local RD[im] = FH[im];
	      #local Up = 1;
	    #else // IndL = 2 
	      #if (IndL != 2) #error concat("// ERROR<<<<<<<<<<<<<<<", str(IndL,5,0), "\n")  #end //if
	      #local GR[re] = FH[re]; #local GR[im] = FH[im];
	      #local RD[re] = FL[re]; #local RD[im] = FL[im];
	      #local Up = 0;
	    #end //#if (IndL = 1)
	    #if ((S - R) = 2 | (S - R) = 3)
	      #if (Up)
		AddRhombus(ThinUp(), "ThinUp", GR, RD)
	      #else
 		AddRhombus(ThinDown(), "ThinDown", GR, RD)
	      #end //#if (Up)
	    #else
	      #if (Up)
		AddRhombus(ThickUp(), "ThickUp", GR, RD)
	      #else
		AddRhombus(ThickDown(), "ThickDown", GR, RD)
	      #end //#if (Up)
	    #end //#if ((S - R) = 2 | (S - R) = 3)
	    #local KS = KS+1;
	  #end //#while (KS<=KSmax) // line of GR2
	  #local S = S+1;
	#end //#while (S<5) // grid2
	#local KR = KR+1;
      #end //#while (KR<=KRmax) // line of grid1
      #local tmp = R;  #local R =  tmp  + 1; // a bug????
    #end // #while (R<5)  /* grid1 */

#end

#macro PentaPlane() 
    union {
      #macro ThinUp()    object { ThinUpRho()   } #end
      #macro ThinDown()  object { ThinDownRho() } #end
      #macro ThickUp()   object { ThickUpRho()  } #end
      #macro ThickDown() object { ThickDownRho()} #end
      union { Penta() }
    }                
#end

// for SCAL = 90: factor 0.5;
// for SCAL = 30: factor 0.5;
#declare EndScale = 3;
 
    
#declare Col = <0,0,0>;                            
            
object { PentaPlane()
    scale EndScale 
    texture { 
      pigment { color White}   
    }
}

//////////////////////////////////////////////////////////////////////////
//                                                                      //
//               CAMERA, LIGHT, BACKGROUND                              //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

camera        
  {                              
    location <100,100,50>
    direction <0,0,1> *3.5
    look_at  <0,-1.5,0>
  }

light_source { <30,200,-60> 
  color White * 1.5
}