#! /usr/bin/perl
#
##### Ammann Beenker tiling by the grid method #####
#
my($i,$j,$j1,$j2,$j3,$k1,$k2,$k3,$l);# integer variables
my(@t);# integer arrays
my($x,$y,$x0,$y0,$xx,$yy,$xx1,$yy1,$xx2,$yy2,
$x1,$x2,$x3,$x4,$y1,$y2,$y3,$y4,
$ex,$ey,$e11,$e12,$e21,$e22,$f11,$f12,$f21,$f22,
$xmin,$ymin,$xmax,$ymax,$eps,$p,$area);# real variables
my(@a0,@a1,@a2);# real arrays
my(%e); # real hash variables
my($s1,$s2,$str); # string variables
my(@s); # string arrays
#
#:: Basic constants
#
$e{'1,1'}=1.0;$e{'2,1'}=0.0;$e{'3,1'}=1.0;$e{'4,1'}=0.0;
#
$e{'1,2'}=sqrt(2.0)/2.0; $e{'2,2'}=sqrt(2.0)/2.0;
$e{'3,2'}=-sqrt(2.0)/2.0;$e{'4,2'}=sqrt(2.0)/2.0;
#
$e{'1,3'}=0.0;$e{'2,3'}=1.0;$e{'3,3'}=0.0;$e{'4,3'}=-1.0;
#
$e{'1,4'}=-sqrt(2.0)/2.0;$e{'2,4'}=sqrt(2.0)/2.0;
$e{'3,4'}=sqrt(2.0)/2.0;$e{'4,4'}=sqrt(2.0)/2.0;
#
$eps=0.001;
#
#:: Make the header of the povray files
#
$s[1] ="// perspective (default) camera";
$s[2] ="camera {";
$s[3] ="  orthographic";
$s[4] ="  location  <0, 0, 7>";
$s[5] ="  look_at   <0, 0, 0>";
$s[6] ="  right     18.2*image_width/image_height";
$s[7] ="  up 18.2";
$s[8] ="}";
$s[9] ="background { color rgb <1,1,1> }";
$s[10]="// create a regular point light source";
$s[11]="light_source {";
$s[12]="  0*x                  // light's position (translated below)";
$s[13]="  color rgb <1,1,1>    // light's color";
$s[14]="  translate <5, 3, 7>";
$s[15]="  // light_source { ...";
$s[16]="  // put this inside a light_source to make it parallel";
$s[17]="  parallel";
$s[18]="  point_at <5, 3, 0>";
$s[19]="}";
$s[20]="";
$s[106] ="  right     3.2*image_width/image_height";
$s[107] ="  up 3.2";
#
open(OUT1,">tiling.pov");
open(OUT2,">grids.pov");
open(OUT3,">window.pov");
#
#:: Make header
#
for($i=1;$i<=20;$i++){
    printf OUT1 "$s[$i]\n";
    printf OUT2 "$s[$i]\n";
    if($i==6||$i==7){
	printf OUT3 "$s[100+$i]\n";
    }else{
	printf OUT3 "$s[$i]\n";
    }
}
#
#:: Projection of the unit cell onto the internal space.
#
for($i=1;$i<=4;$i++){
#
    $j1=$i  % 4 + 1;
    $j2=$j1 % 4 + 1;
    $j3=$j2 % 4 + 1;
#
    for($k1=0;$k1<=1;$k1++){
	for($k2=0;$k2<=1;$k2++){
	    for($k3=0;$k3<=1;$k3++){
#
		$s1="3,$j1";$s2="3,$j2";$s3="3,$j3";
		$xx=$e{$s1}*$k1+$e{$s2}*$k2+$e{$s3}*$k3;
		$s1="4,$j1";$s2="4,$j2";$s3="4,$j3";
		$yy=$e{$s1}*$k1+$e{$s2}*$k2+$e{$s3}*$k3;
#
		$s1="3,$i";$xx1=$xx+$e{$s1};
		$s1="4,$i";$yy1=$yy+$e{$s1};
#
		if($i==1){
		    $str ="0 $k1 $k2 $k3";
		    $str1="1 $k1 $k2 $k3";
		}elsif($i==2){
		    $str ="$k3 0 $k1 $k2";
		    $str1="$k3 1 $k1 $k2";
		}elsif($i==3){
		    $str ="$k2 $k3 0 $k1";
		    $str1="$k2 $k3 1 $k1";
		}elsif($i==4){
		    $str ="$k1 $k2 $k3 0";
		    $str1="$k1 $k2 $k3 1";
		}
#
		$xx=-$xx;   $yy=-$yy;
		$xx1=-$xx1; $yy1=-$yy1;
		printf OUT3 "cylinder{ <%7.5f,%7.5f,0>,",$xx,$yy;
		printf OUT3 "<%7.5f,%7.5f,0>,0.01 pigment{rgb<1,.8,.8>}}\n",$xx1,$yy1;
#
	    }
	}
    }
#
}
#
#:: Border lines of the box
#
$xmin=-15.0;$ymin=-15.0;
$xmax=15.0;$ymax=15.0;
#$xmin=-6.0;$ymin=-4.0;
#$xmax=6.0;$ymax=4.0;
#printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xmin,$ymin;
printf OUT2 "mesh2{vertex_vectors{4, <%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0> } face_indices{2, <0,1,2>,<0,2,3> }  pigment {rgbf <0,0,0>}}\n",$xmin,$ymin,$xmax,$ymin,$xmax,$ymax,$xmin,$ymax;
#
#:: Draw four sets of grids (interval = 2.0)
#
$grinterval=2.0;
$shift[1]=0.0;
$shift[2]=0.5;
$shift[3]=-0.55;
$shift[4]=0.77;
for($j=1;$j<=4;$j++){ # No.j grids
#
    $s1="1,$j"; $s2="2,$j";
    $ex=$e{$s1};$ey=$e{$s2};
#
    for($i=-10;$i<=10;$i++){
#
	if(abs($ex)>$eps && abs($ey)>$eps){
	    $x1=(($i+$shift[$j])*$grinterval-$ey*$ymin)/$ex;$y1=$ymin;
	    $x2=(($i+$shift[$j])*$grinterval-$ey*$ymax)/$ex;$y2=$ymax;
	    $x3=$xmin;$y3=(($i+$shift[$j])*$grinterval-$ex*$xmin)/$ey;
	    $x4=$xmax;$y4=(($i+$shift[$j])*$grinterval-$ex*$xmax)/$ey;
	} elsif(abs($ex)<$eps && abs($ey)>$eps){
	    $x3=$xmin;$y3=(($i+$shift[$j])*$grinterval)/$ey;
	    $x4=$xmax;$y4=$y3;
	    $x1=-100.0; $y1=$y3;
	    $x2=100.0; $y2=$y3;
	} elsif(abs($ex)>$eps && abs($ey)<$eps){
	    $x1=(($i+$shift[$j])*$grinterval)/$ex;$y1=$ymin;
	    $x2=$x1;$y2=$ymax;
	    $x3=$x1;$y3=-100.0;
	    $x4=$x1;$y4=100.0;
	}
#
	if($x1>=$xmin-$eps && $x1<=$xmax+$eps && 
	   $x2>=$xmin-$eps && $x2<=$xmax+$eps ){
	    $xx1=$x1; $yy1=$ymin;
	    $xx2=$x2; $yy2=$ymax;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	} elsif($y3>=$ymin-$eps && $y3<=$ymax+$eps && 
		$y4>=$ymin-$eps && $y4<=$ymax+$eps ){
	    $xx1=$xmin; $yy1=$y3;
	    $xx2=$xmax; $yy2=$y4;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	} elsif($x1>=$xmin-$eps && $x1<=$xmax+$eps && 
		$y4>=$ymin-$eps && $y4<=$ymax+$eps ){
	    $xx1=$x1; $yy1=$ymin;
	    $xx2=$xmax; $yy2=$y4;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	} elsif($x1>=$xmin-$eps && $x1<=$xmax+$eps && 
		$y3>=$ymin-$eps && $y3<=$ymax+$eps ){
	    $xx1=$x1; $yy1=$ymin;
	    $xx2=$xmin; $yy2=$y3;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	} elsif($x2>=$xmin-$eps && $x2<=$xmax+$eps && 
		$y3>=$ymin-$eps && $y3<=$ymax+$eps ){
	    $xx1=$x2; $yy1=$ymax;
	    $xx2=$xmin; $yy2=$y3;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	} elsif($x2>=$xmin-$eps && $x2<=$xmax+$eps && 
		$y4>=$ymin-$eps && $y4<=$ymax+$eps ){
	    $xx1=$x2; $yy1=$ymax;
	    $xx2=$xmax; $yy2=$y4;
	    printf OUT2 "cylinder{ <%7.5f,%7.5f,0>,",$xx1,$yy1;
	    printf OUT2 "<%7.5f,%7.5f,0>,0.03 pigment{rgb<1,.8,.8>}}\n",$xx2,$yy2;
	}
    }
#
}
#
#:: List up all the intersections within the box
#
for($j1=1;$j1<=3;$j1++){
    $s1="1,$j1"; $s2="2,$j1";
    $e11=$e{$s1};$e12=$e{$s2};
    for($j2=$j1+1;$j2<=4;$j2++){
	$s1="1,$j2"; $s2="2,$j2";
	$e21=$e{$s1};$e22=$e{$s2};
#
	$det=$e11*$e22-$e12*$e21;
	$f11=$e22/$det;
	$f12=-$e12/$det;
	$f21=-$e21/$det;
	$f22=$e11/$det;
# inverse matrix : f
#
	for($i1=-10;$i1<=10;$i1++){
	    for($i2=-10;$i2<=10;$i2++){
		$xx=($i1+$shift[$j1])*$grinterval;
		$yy=($i2+$shift[$j2])*$grinterval;
#
		$xx1=$f11*$xx+$f12*$yy;
		$yy1=$f21*$xx+$f22*$yy;
#
		if($xx1>=$xmin-$eps && $xx1<=$xmax+$eps &&
		   $yy1>=$ymin-$eps && $yy1<=$ymax+$eps){
		    printf OUT2 "sphere {<%10.6f,%10.6f,0.0>,0.015 pigment {color rgb <0.7,0.7,0.3>}}\n", $xx1,$yy1;
#
#                   indices for the tile vertices
		    for($j=1;$j<=4;$j++){
			if($j==$j1){
			    $t[$j]=$i1;
			}elsif($j==$j2){
			    $t[$j]=$i2;
			}else{
			    $s1="1,$j"; $s2="2,$j";
			    $p=($xx1*$e{$s1}+$yy1*$e{$s2})/$grinterval;
			    $p=$p-$shift[$j];
			    if($p<-$eps){
				$t[$j]=int($p-1);
			    } else{
				$t[$j]=int($p);
			    }
			}
		    }
#
		    $x1=0.0; $y1=0.0;
		    for($j=1;$j<=4;$j++){
			$s1="1,$j"; $s2="2,$j";
			$x1=$x1+($t[$j]+$shift[$j])*$e{$s1};
			$y1=$y1+($t[$j]+$shift[$j])*$e{$s2};
		    }
#
		    $s1="1,$j1"; $s2="2,$j1";
		    $x2=$x1-$e{$s1};
		    $y2=$y1-$e{$s2};
#
		    $s1="1,$j2"; $s2="2,$j2";
		    $x3=$x2-$e{$s1};
		    $y3=$y2-$e{$s2};
#
		    $x4=$x1-$e{$s1};
		    $y4=$y1-$e{$s2};
#
		    $xx1=$x2-$x1; $yy1=$y2-$y1;
		    $xx2=$x4-$x1; $yy2=$y4-$y1;
		    $area=abs($xx1*$yy2-$yy1*$xx2);
#
		    if($area<0.9){
			printf OUT1 "mesh2{vertex_vectors{4, <%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0> } face_indices{2, <0,1,2>,<0,2,3> }  pigment {rgbf <0,0.8,0>}}\n",$x1,$y1,$x2,$y2,$x3,$y3,$x4,$y4;
		    }else{
			printf OUT1 "mesh2{vertex_vectors{4, <%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0>,<%7.5f,%7.5f,0> } face_indices{2, <0,1,2>,<0,2,3> }  pigment {rgbf <1,0.6,0>}}\n",$x1,$y1,$x2,$y2,$x3,$y3,$x4,$y4;
		    }
		    printf OUT1 "cylinder{ <%7.5f,%7.5f,0>,",$x1,$y1;
		    printf OUT1 "<%7.5f,%7.5f,0>,0.02 pigment{rgb<0,0,0>}}\n",$x2,$y2;
		    printf OUT1 "cylinder{ <%7.5f,%7.5f,0>,",$x2,$y2;
		    printf OUT1 "<%7.5f,%7.5f,0>,0.02 pigment{rgb<0,0,0>}}\n",$x3,$y3;
		    printf OUT1 "cylinder{ <%7.5f,%7.5f,0>,",$x3,$y3;
		    printf OUT1 "<%7.5f,%7.5f,0>,0.02 pigment{rgb<0,0,0>}}\n",$x4,$y4;
		    printf OUT1 "cylinder{ <%7.5f,%7.5f,0>,",$x4,$y4;
		    printf OUT1 "<%7.5f,%7.5f,0>,0.02 pigment{rgb<0,0,0>}}\n",$x1,$y1;
#
		    $x1=0.0; $y1=0.0;
		    for($j=1;$j<=4;$j++){
			$s1="3,$j"; $s2="4,$j";
			$x1=$x1+($t[$j]+$shift[$j])*$e{$s1};
			$y1=$y1+($t[$j]+$shift[$j])*$e{$s2};
		    }
#
		    $s1="3,$j1"; $s2="4,$j1";
		    $x2=$x1-$e{$s1};
		    $y2=$y1-$e{$s2};
#
		    $s1="3,$j2"; $s2="4,$j2";
		    $x3=$x2-$e{$s1};
		    $y3=$y2-$e{$s2};
#
		    $x4=$x1-$e{$s1};
		    $y4=$y1-$e{$s2};
#
		    printf OUT3 "sphere {<%10.6f,%10.6f,0.0>,0.02 pigment {color rgb <0.7,0.7,0.3>}}\n", $x1,$y1;
		    printf OUT3 "sphere {<%10.6f,%10.6f,0.0>,0.02 pigment {color rgb <0.7,0.7,0.3>}}\n", $x2,$y2;
		    printf OUT3 "sphere {<%10.6f,%10.6f,0.0>,0.02 pigment {color rgb <0.7,0.7,0.3>}}\n", $x3,$y3;
		    printf OUT3 "sphere {<%10.6f,%10.6f,0.0>,0.02 pigment {color rgb <0.7,0.7,0.3>}}\n", $x4,$y4;
#
		}
	    }
	}
#
    }
}
#
close(OUT1);
close(OUT2);
close(OUT3);
#
