#!/usr/bin/perl -w
#cxe- spirit - 99.10.28
#cxe- This is a special purpose program designed for a very special thing-
#cxe- to tile a sphere with trapezoidal shapes. You really have to see the 
#cxe- final product or be me to make sense of it, BUT there are lots of
#cxe- good things to look at here. I have lots of general geometry
#cxe- manipulation functions that I'm sure I'll be using again.
#cxe- This program should load data from a specially formatted vector file.
#cxe- The format can best be parsed from AutoCAD's .dxf format and
#cxe- specifically focuses on pline entities. This is the format of the
#cxe- data file:
#cxe- POLYLINE
#cxe- VERTEX -0.12345678 -0.12345678
#cxe- VERTEX 0.12345678 0.12345678
#cxe-  .. .. ..
#cxe- SEQEND CLOSED
#cxe- Note that the "closed" is optional and will be handy when
#cxe- reconstructing the entities.
#cxe- Actually, I just realized that it doesn't matter what you do.
#cxe- the non point stuff is irrelavent and gets totally ignored. Do
#cxe- whatever you need to do with that stuff. But if it's a point
#cxe- and needs to be dealt with, make sure it starts with "VERTEX".
#cxe- Beware, however, since the "VERTEX" gets stripped out in the end.

#cxe-       =============== SETUP ===============
#cxe- Set up the number pi...
$pi = atan2(0,-1);    #cxe-  Why isn't this better than Larry's way?
$DIGITS = 8;	      #cxe-  How many decimals to output - set it here.
$VERBOSE = "0"; #cxe-  Set this to something for a printout of data

#cxe- User inputs filename here.
print STDOUT "Enter the filename to process: " ;
chop($ORIGFILE=<STDIN>) ;

#cxe- Here's a tricky bit. This copies the specified file to a temporary 
#cxe- file specific to this program using $$, the PID. The 2>&1 part takes
#cxe- any STDERR and pipes it to STOUT. The theory here is that if the copy
#cxe- is happy, no ERR will show up and we can use that as a test to see if
#cxe- there were any problems. If not (the else) open the tempfile for use.
if ( `cp $ORIGFILE spirit.$$.tmp 2>&1`)
{ print ("Sorry, file problems.\n") and die;}
else 
{
$DAFILE="spirit.$$.tmp";
#TESTER- print ("Check out: ", $DAFILE, "\n");
$DAFILE= "+<".$DAFILE; #cxe-  Hmm just in case I need to do so I and O.
}

#cxe-       =============== TESTING CALLS ===============
if ( $VERBOSE ) {  #cxe- Prints out the raw data first if we want..
 GETDATA();
 $COUNTER=0;
 foreach $VERTEX (@VERTICES)
 { print ( $COUNTER++ , "==>" , $VERTEX ,"<==\n" ) ; }
}

#cxe- Some subroutine testing...
#TESTER- SCALE (@VERTICES, .5);
#TESTER- ROTATE (@VERTICES, .785398163397);
#TESTER- TRANSLATE (@VERTICES, -5,5,3);
#TESTER- TAPER (@VERTICES, 1, 1, .5); 
#TESTER- print ("Here's CIRCIRINT->", CIRCIRINT (2.2, -5.3, 1.5), "<-\n" );
#TESTER- ROUND (@VERTICES, 3);

#TESTER- foreach $VERTEX (@VERTICES)
#TESTER- { print ("New GLOBAL vector! -->", $VERTEX, "\n"); }

#TESTER- die "\n\n...Ok, show\'s over. Move along...\n";


#cxe-       =============== MAIN PROCESSING ===============
$FACETS = 20; #cxe- Number of shpere tiles to create
$ALPHA = 0; #cxe- This is the main angle from the XY plane to the insert point
$BASELENGTH = 1; #cxe- This starts at 1 and defines the size of the sphere which
	       #cxe- must be perfectly circumscribed around a 24-sided polygon
	       #cxe- with faces that equal this. It can only get smaller and
	       #cxe- will be used as a scale factor since shapes should be 1.
$BHFACTOR = .5; #cxe- The ratio of baselength to base height
$BASEHEIGHT = $BHFACTOR * $BASELENGTH;
              #cxe- Baseheight starts off at .5 since 1st len is one. 

$ZOIDRATIO = ( ($BASELENGTH - (2*$BASEHEIGHT*tan($pi/24))) / $BASELENGTH ); 
             #cxe-  This is the
             #cxe- constant ratio between the top of a trapezoid with base of
             #cxe- BASELENGTH, a vanishingpoint angle of 15 degrees (pi/24),
             #cxe- and a height of BASEHEIGHT.
$SPHERERAD = (($BASELENGTH / 2) / tan($pi/24)) ; #cxe-  This is kind of an
	  #cxe- approximation since this is the middle of the shape (lowest
          #cxe- point). But the other points should work out too. It is
          #cxe- my main control arc's radius.

for ($I = 1; $I <= $FACETS; $I++) {

#cxe- Get refreshed stock of data
GETDATA();

#cxe- Calculate base point (X is of course 0)
$BASEPOINTY = -$SPHERERAD * cos($ALPHA);
$BASEPOINTZ = $SPHERERAD * sin($ALPHA);


#cxe- Scale the information so that the main baseline is the right fit.
#cxe- A bit redundant on the first one, but it makes mini-versions as
#cxe- it moves along.
#TESTER- print ("Test just before Scaling. BL=", $BASELENGTH, " #1=");
#TESTER- print ($VERTICES[0],"\n");
SCALE (@VERTICES, $BASELENGTH);
#TESTER- print ("Test just after Scaling. BL=", $BASELENGTH, " #1=");
#TESTER- print ($VERTICES[0],"\n");

#cxe- Call CIRCIRINT to find the angle that this block lays on
#cxe- Its parameters read (main circle X, main circle Y, small rad)
$LEANANGLE =  CIRCIRINT( (-$BASEPOINTY), (-$BASEPOINTZ) , $BASEHEIGHT) ;

#cxe- Find offsets for the next basepoint which is this top midpoint.
$DELTAY = $BASEHEIGHT * cos($LEANANGLE);
$DELTAZ = $BASEHEIGHT * sin($LEANANGLE);

#cxe- And the point itself...
$NEXTBASEPTY = ($BASEPOINTY + $DELTAY);
$NEXTBASEPTZ = ($BASEPOINTZ + $DELTAZ);
$NEXTALPHA = atan2($NEXTBASEPTZ,abs($NEXTBASEPTY));

#cxe- Find next baselength which is this top taper length.
$NEXTBASELEN = ( $NEXTBASEPTY * tan($pi/24) * -2 );
#cxe- This is what the trapeziod starts at it's default shape that
#cxe- needs to be expanded a bit. The difference between NEXTBASELEN
#cxe- and TOPISNOW is the whole reason for the TAPER function.
$TOPISNOW = ( $BASELENGTH * $ZOIDRATIO);

#cxe- Do the taper with ratio of "next base length" to "top as it is now" to
#cxe- fill out the top so it becomes the same size as the bottom of the next. 
#cxe- Base scale is 1 and the height is baseheight. 
#cxe- TAPER is called with (oh, and VERT,grrr) height, bottom, top
#TESTER-print ("Top factor:", $NEXTBASELEN/$TOPISNOW , "\n");
#TESTER-print ("BaseHeight:",$BASEHEIGHT,"\n");
TAPER (@VERTICES, $BASEHEIGHT, $NEXTBASELEN, $TOPISNOW );

#cxe-  Rotate the points to their new angle
ROTATE( @VERTICES, ($LEANANGLE) );

#cxe-  Move the points to their new home
TRANSLATE (@VERTICES, 0, $BASEPOINTY, $BASEPOINTZ);

#cxe- Clean up the numbers to be a nice size
ROUND(@VERTICES, $DIGITS);

#TESTER- foreach $VERTEX (@VERTICES)
#TESTER- { print ("New GLOBAL vector! -->", $VERTEX, "\n"); }

#cxe- It's a wrap...
#cxe- I added this "feature" to skip writing data for every other panel.
#cxe- This is so I can run it 4x and get scripts with each image in each
#cxe- position. Change "if" to "unless" for the even ones.
#cxe- And, to get every facet, just comment out this control.
if ( $I%2 )
 {
 PUTDATA();
 }

#cxe- ok, we're just about done with this facet. Lets pass on the King's
#cxe- crown, scpetre and throne to the Prince.
$ALPHA = $NEXTALPHA;
$BASELENGTH = $NEXTBASELEN;
$BASEHEIGHT = $BASELENGTH * $BHFACTOR;

} #cxe- end of main for loop - another facet is done

#cxe- Remove the temp file or wonder about it.
$DAFILE =~ s/\+<(.*)/$1/;
if (`rm $DAFILE`){
print ("Hmm...had a spot of trouble getting rid of the temp file ");
print ("called:\n", $DAFILE, " Know what I mean?\n");
}

die "\n\n...Ok, show\'s over. Move along...\n";
#cxe- End of MAIN

#cxe-       =============== SUBROUTINES ===============
#cxe- GETDATA subroutine- This brings in the data from the user specified
#cxe- file. I made it a function so that I could call it repeatedly to 
#cxe- periodically refresh my data for each panel.
#cxe- Zoom through file and put all vertices in an array. Also strip off
#cxe- text and crap and just leave space separated number data. I've also
#cxe- gone ahead and added a Z value of 0.0 since this data is assumed to
#cxe- to be coming in without it.
sub GETDATA {
@VERTICES = 0;
open DAFILE or die ;
while ( <DAFILE> ) {
  if ( /VERTEX/ ) {
   chop;
   s/VERTEX (.*)/$1 0\.0/ ;
   push @VERTICES, $_ ;
   }
}
close DAFILE;
return
} #cxe-  end of sub GETDATA

#cxe- PUTDATA subroutine- This goes through the original file and
#cxe- whenever it sees the word VERTEX, it replaces the whole line with
#cxe- the next item in the processed @VERTICES array. Theoretically, there
#cxe- should be one VERTEX entry for each $VERTICES[n] entry and they
#cxe- should be related.
sub PUTDATA {
$NEWFILE = ">>".$ORIGFILE.".out";
open DAFILE or die "Hmm...couldn't open the datafile during final output";
open NEWFILE or die "Hmm...couldn't open the out put file...";
print ("Writing to file --", $NEWFILE,"\n");
print ("Outputting another facet...\n");
my $I = 1;
while ( <DAFILE> ) {
 if ( /VERTEX/ ) { 
 print ("VERTICES[", $I,"] = ", $VERTICES[$I], "\n");
 print NEWFILE ($VERTICES[$I++], "\n");
 }
 else  {
 print NEWFILE ($_ );
 }
}
close DAFILE;
close NEWFILE;
print ("\n");
return;
} #cxe-  end of sub PUTDATA

#cxe- CIRCIRINT subroutine- This function finds a particular point
#cxe- where two special case circles intersect. It's very focused to what I
#cxe- need (a generalized one would be nice). Basically this function
#cxe- compares a quarter arc that is in quadrant 2 and a small quarter arc
#cxe- of quadrant 1 whose center point is on the first arc. Basically if I
#cxe- lay a penny down on a clock edge at around the 10:00 point, I'm
#cxe- looking for the point where the penny intersects the clockedge toward
#cxe- midnight - not towards nine. The function doesn't return a point
#cxe- (that's not what I need). Instead, it returns the angle from the
#cxe- center of the penny to the upper right intersection point. The
#cxe- parameters are the X and Y points of the "clock" center and the
#cxe- diameter of the penny. The centerpoint of the penny on the edge of the
#cxe- clock will be assumed to be 0,0. The strategy here is to just bisect
#cxe- the maximum and minimum possibilities until we get to a point where
#cxe- the distance from the proposed point is close enough to the diameter
#cxe- of the clock to assure us that we are virtually on its edge.

sub CIRCIRINT
{
#cxe- So parameters read (main circle X, main circle Y, small radius)
my $ITERATIONS = 5000;
my $OKERROR = .000001; #cxe- CHANGE THIS TO SUIT REQUIREMENTS
my $SMCIR_RAD = pop; 
my $LGCIR_Y = pop; #cxe-  should be negative
my $LGCIR_X = pop; #cxe-  should be positive
my $LOTRY = 0; 		#cxe-  starts out horizontal pointing into the clock
my $HITRY = $pi / 2;	#cxe-  vertical outside of the clock
#cxe- Although this is a good general guess, it might be smarter to start
#cxe- guessing at a low angle or so since the values are asymtotic to 0.
my $AVGTRY = $pi / 4;	#cxe-  a good starting guess (45deg)

#TESTER- print ("AVGTRY=", $AVGTRY, " : ");
#TESTER- print ("SMCIR_RAD=", $SMCIR_RAD, " : ");
#TESTER- print ("LGCIR_Y=", $LGCIR_Y, " : ");
#TESTER- print ("LGCIR_X=", $LGCIR_X, "\n");

my ( $SM_2_LG, $TEMPX, $TEMPY, $ERROR, $I );

#cxe- distance formula from 2 circles.
$SM_2_LG = dist(0, 0, $LGCIR_X, $LGCIR_Y);

for ( $I=1; $I <= $ITERATIONS; $I++ )
{
  $TEMPX = cos($AVGTRY) * $SMCIR_RAD;
  $TEMPY = sin($AVGTRY) * $SMCIR_RAD;
  $ERROR = ($SM_2_LG - dist($TEMPX,$TEMPY,$LGCIR_X,$LGCIR_Y));
#TESTER- print ("I=", $I, " : ");
#TESTER- print ("AVGTRY=", $AVGTRY, " : ");
#TESTER- print ("SMCIR_RAD=", $SMCIR_RAD, " : ");
#TESTER- print ("TEMPX=", $TEMPX, " : ");
#TESTER- print ("TEMPY=", $TEMPY, " : ");
#TESTER- print ("ERROR=", $ERROR, "\n");

  if ( abs($ERROR) < $OKERROR ) 
  { return $AVGTRY; }
  elsif ( $ERROR <  0 )
  {
   $HITRY = $AVGTRY;
   $AVGTRY = (($AVGTRY + $LOTRY) / 2); 
  }
  else				 #cxe- it must be less than 
  { $LOTRY = $AVGTRY;
    $AVGTRY = (($AVGTRY + $HITRY) / 2);
  }
}
print ("Gosh, I ran ",$ITERATIONS," times and still didn't");
print (" solve within target.\n");
print ("The error was down to: ", $ERROR,"\n");
print ("The SM_2_LG variable (big radius) is: ", $SM_2_LG,"\n");
print ("The closest to that I got was: ");
print (dist($TEMPX,$TEMPY,$LGCIR_X,$LGCIR_Y),"\n");
return $AVGTRY;

} #cxe- end of sub CIRCIRINT


#cxe- TAPER subroutine- WOW! HERE IT IS! THE function that I actually
#cxe- need that AutoCAD utterly lacks. Here's what this does; it basically
#cxe- scales points in one axis an amount which is dependant on another
#cxe- axis. I've kept it to the simple case that I need which is scaling the
#cxe- X values based on their Y position. In simpler words, imagine a
#cxe- collection of points that form an image that's roughly square. With
#cxe- this function, you can "squeeze" the sides of the square as if you
#cxe- were pinching it toward the bottom or top. The effect on an actual
#cxe- square would be a perfect trapeziod. The inputs here are the data
#cxe- points in an array of "X Y Z" strings, followed by the X scale factor
#cxe- at Y=1. The X scale factor at Y=0 is assumed to be 1. So let's say you
#cxe- had a square that was 1x1 that you needed to turn into a trapezoid
#cxe- with a bottom of 1 and a top of .2; you'd use a scale factor of .2. If
#cxe- the square was 5x5 and you needed the top to be 3.5, you could scale
#cxe- the top in using .14 which is (3.5top/5base)/5height.
sub TAPER
{
#cxe- in factor formula topedge then bottomedge then height
#cxe- so in call on "stack" it's backwards: height, bottom, top
my $TIN = pop(); #cxe- $TopIsNow (what is the top like to start with)
my $NBL  = pop(); #cxe- $NextBaseLen (what it should be)
my $BH = pop(); #cxe-How tall is this trapeziod, probably smaller.

foreach $VERTEX (@_)
 {
 my $X = $VERTEX;	#cxe-  Extract the X component
 $X =~ s/^(.*) .* .*$/$1/;
 my $Y = $VERTEX;	#cxe-  Extract the Y component (shouldn't change)
 $Y =~ s/^.* (.*) .*$/$1/;
 my $Z = $VERTEX; #cxe-  Extract the Z component for use in reconstruction
 $Z =~ s/^.* .* (.*)$/$1/;

#TESTER- print ("Y:",$Y," BH:",$BH," NBL:",$NBL," TIN:",$TIN,"\n");
my $FACTOR = (( ($Y/$BH) * ( ($NBL/$TIN) -1 ) ) +1) ; 

 $X *= $FACTOR; #cxe- That should do it. 

 $VERTEX = $X." ".$Y." ".$Z; 
 }
} #cxe- end of sub TAPER



#cxe- TRANSLATE subroutine- This function translates coordinates from one
#cxe- place to another. It is fed an array of coordinate strings ("-0.12 -0.12
#cxe- -0.12") followed by three offsets, one for X, Y, and Z. This function also
#cxe- could have some of the axes removed for optimization, but I went ahead and
#cxe- set it up for the general case.
sub TRANSLATE
{
my $ZOFFSET = pop;
my $YOFFSET = pop;
my $XOFFSET = pop;

foreach $VERTEX (@_)
 {
 my $X = $VERTEX;	#cxe-  Extract the X component
 $X =~ s/^(.*) .* .*$/$1/;
 $X += $XOFFSET;
 my $Y = $VERTEX;	#cxe-  Extract the Y component
 $Y =~ s/^.* (.*) .*$/$1/;
 $Y += $YOFFSET;
 my $Z = $VERTEX;	#cxe-  Extract the Z component
 $Z =~ s/^.* .* (.*)$/$1/;
 $Z += $ZOFFSET;
 $VERTEX = $X." ".$Y." ".$Z; 
 }
} #cxe-  end of sub TRANSLATE

#cxe- ROUND subroutine- This is a simple subroutine that does ordinary
#cxe- rounding. I didn't see it listed as a built-in function, but
#cxe- I wanted clean numbers. This function takes an array of point
#cxe- strings ("-0.1234 -0.1234 -0.1234") followed by the number of
#cxe- places you want them rounded to.
#cxe- It appears to work...A nice addition might be padding too.
sub ROUND
{
my $PLC = pop;
$PLC = 10 ** $PLC; #cxe- 3places=1000, 5=100000 etc.
foreach $VERTEX (@_)
 {
 my $X = $VERTEX;	#cxe-  Extract the X component
 $X =~ s/^(.*) .* .*$/$1/;

 ($X == 0) || ($X =  ( int($X*$PLC+    (.5*($X/abs($X)))    )/$PLC ));

 my $Y = $VERTEX;	#cxe-  Extract the Y component
 $Y =~ s/^.* (.*) .*$/$1/;
 ($Y == 0) || ($Y =  ( int($Y*$PLC+    (.5*($Y/abs($Y)))    )/$PLC ));

 my $Z = $VERTEX;	#cxe-  Extract the Z component
 $Z =~ s/^.* .* (.*)$/$1/;
 ($Z == 0) || ($Z =  ( int($Z*$PLC+    (.5*($Z/abs($Z)))    )/$PLC ));

 $VERTEX = $X." ".$Y." ".$Z; 
 }
} #cxe-  end of sub ROUND


#cxe- SCALE subroutine- This function executes a linear scale based at 0,0,0.
#cxe- This is really quite simple as it just multiplies every value by a
#cxe- factor. This function is called with an array of space separated
#cxe- coordinates "-0.123 -0.123 -0.123" and a scale factor. It basically goes
#cxe- through and picks out each X,Y, and Z value and multiplies the factor to
#cxe- it. Then it reconstitutes it back to the global that called it. Some
#cxe- interesting enhancements could easily include separate scale factors for
#cxe- each axis. Another thing is an optimization - the Z scale is a waste, but
#cxe- I kept it in here for generality in case I use this again later. The
#cxe- various axes can be separated out for specific optimized scalings. They
#cxe- could all three be separated and called by an intermediate subroutine
#cxe- which would have arguments directing the effective axes.
sub SCALE
{
my $FACTOR = pop;
($FACTOR <= 0) && die "Whoops - you can't have a negative scale factor.";
#TESTER- print ("Start sub SCALE.\n");
#cxe- Sets up loop variable $VERTEX from list passed to sub (@_).
foreach $VERTEX (@_)
{
 my $X = $VERTEX;	#cxe-  Extract the X component
 $X =~ s/^(.*) .* .*$/$1/;
 $X *= $FACTOR;
 my $Y = $VERTEX;	#cxe-  Extract the Y component
 $Y =~ s/^.* (.*) .*$/$1/;
 $Y *= $FACTOR;
 my $Z = $VERTEX;	#cxe-  Extract the Z component
 $Z =~ s/^.* .* (.*)$/$1/;
 $Z *= $FACTOR;
$VERTEX = $X." ".$Y." ".$Z; 

#TESTER- print ("New vector -->", $VERTEX, "\n");
}
#TESTER- print ("Done with sub SCALE\n");
} #cxe-  end of sub SCALE

#cxe- ROTATE subroutine - Note that this is just implemented to rotate about
#cxe- the X axis (which is kind of odd); that can be changed pretty easily.
#cxe- It's also not very robust since I just need it to rotate within 1
#cxe- quadrant right now. That can be fixed up later.
#cxe- It takes a list of point strings followed by a radian angle.
sub ROTATE
{
my $ANGLE = pop;
foreach $VERTEX (@_)
 {
 my $X = $VERTEX;	#cxe-  Extract the X component for future reconstruction
 $X =~ s/^(.*) .* .*$/$1/;
#TESTER- print ("ANGLE: X=", $X);
 my $Y = $VERTEX;	#cxe-  Extract the Y component
 $Y =~ s/^.* (.*) .*$/$1/;
#TESTER- print (" Y=", $Y);
 my $Z = $VERTEX;	#cxe-  Extract the Z component
 $Z =~ s/^.* .* (.*)$/$1/;
#TESTER- print (" Z=", $Z);
 my $DIST = dist( $Y, $Z, 0, 0 );
#TESTER- print (" DIST=", $DIST);
 my $ANGLE1 = atan2($Z,$Y);
#TESTER- print (" ANGLE1=", $ANGLE1);
 $ANGLE1 += $ANGLE;
#TESTER- print (" New ANGLE1=", $ANGLE1);
 $Y = cos($ANGLE1) * $DIST;
#TESTER- print (" New Y=", $Y);
 $Z = sin($ANGLE1) * $DIST;
#TESTER- print (" New Z=", $Z);
 $VERTEX= $X." ".$Y." ".$Z;
#print (" New VERTEX=",$VERTEX,"\n");
 }
} #cxe-  end of sub ROTATE

#cxe- The distance formula for 2d points.
sub dist { sqrt( ( $_[0] - $_[2] )**2 + ( $_[1] - $_[3] )**2); }

#NOTE- atan2(Y,X) = builtin
#NOTE- sin(rad) = builtin
#NOTE- cos(rad) = builtin
#NOTE- sqrt(X) = builtin
#NOTE- $pi = atan2(1,1) * 4; #cxe- That's the way it is in the Camel book...
sub tan { sin($_[0]) / cos($_[0]) }
sub acos { atan2( sqrt( 1 - $_[0] * $_[0]), $_[0] ) }
sub asin { atan2( $_[0], sqrt( 1 - $_[0] * $_[0])) }