PageRenderTime 190ms CodeModel.GetById 3ms app.highlight 176ms RepoModel.GetById 1ms app.codeStats 0ms

/mj/HalfOctant4.pl

https://bitbucket.org/maphew/cahill-keyes
Perl | 1202 lines | 858 code | 91 blank | 253 comment | 179 complexity | 2d06a06c705e6d799fdf81a217d8321a MD5 | raw file
   1#! /usr/bin/perl -w
   2#	Program HalfOctant4. (made from HalfOctant3.odt on 2010-04-03)
   3#
   4#	In this version, all parallels are equidistant from equator to pole, along meridians.
   5# There are no polar circular arcs. In other HalfOctant#.odt versions, dP was a hash, %dP,
   6# with parallel spacings for polar, supple, temperate and torrid zones. In this version, it is
   7# only an array, @dP, with a single parallel spacing value for each meridian.
   8
   9use Math::Trig;
  10
  11# Some constants for use by Subroutine MJtoG, which does coordinate axis transformation.
  12# Angle of rotation is 60�.	Point G is (10000,0) in MJ, and (5000,0) in G
  13$cos60 = cos (deg2rad(60));
  14$sin60 = sin (deg2rad(60));
  15$yTranslate = 10000 * $sin60;
  16
  17# Multidimensional arrays are done with hashes, with the keys used for i,j indeces.
  18
  19# Given input
  20$xM = 0;		$yM = 0;	# Point M is the origin of the axes
  21$MA = 940;				# Indent of octant's apex, with respect to triangle's
  22$xG = 10000;	$yG = 0;	# Point G, at center of base of octant
  23$DeltaFrigid = 104;			# Distance for each degree of latitude close to the pole
  24
  25# Other points and lengths of special interest
  26$xA = $MA;		$yA = 0;						# Point A, at the indent
  27$xN = $xG;		$yN = $xG * tan (deg2rad(30));			# Point N
  28($flag, $xB, $yB) = LineIntersection ($xM,$yM,30,$xA,$yA,45);	# Point B
  29if ($flag ne "Intersect") { die "$flag while calculating point B.\n"; }
  30$AG = $xG - $xA;
  31$AB = Length ($xA, $yA, $xB, $yB);
  32$EF = $AB;						# Sometimes I prefer AB, others EF
  33$MB = Length ($xM, $yM, $xB, $yB);
  34$MN = Length ($xM, $yM, $xN, $yN);
  35$xD = Interpolate ($MB, $MN, $xN, $yM);	# D is away from N as B is away from M
  36$yD = Interpolate ($MB, $MN, $yN, $yM);
  37$xF = $xG;
  38$yF = $yN - $MB;
  39$xE = $xN - $MA * sin (deg2rad(30));
  40$yE = $yN - $MA * cos (deg2rad(30));
  41$FG = $yF;
  42$BD = Length ($xB, $yB, $xD, $yD);
  43$EFG = $EF + $FG;
  44
  45# Calculate polar start, two joints, and equatorial end, for each meridian; hashes %xJ, %yJ.
  46# Also calculate a few other values for other hashes, on the way.
  47
  48$DeltaMEq = $EFG / 45;	# 45 meridian spacings along equator for half an octant
  49
  50foreach $m (0 .. 45) {	# Loop for each meridian.
  51	$M = $m . ',' ;	# part of hash key (array row number)
  52
  53	# Start point of minor meridians will be modified to be at parallel 85� after that
  54	# parallel has been calculated. For now, all meridians are set to start at point A.
  55	$xJ {$M . 0} = $xA;		$yJ {$M . 0} = $yA;
  56
  57	# End point of meridians, that is, equidistant intersections with equator;
  58	# These are also the meridian crossings for the zero-th parallel
  59	$distance = $m * $DeltaMEq;
  60	if ($distance <= $yF) {
  61		$xJ {$M . 3} = $xP {$M . 0} = $xG;
  62		$yJ {$M . 3} = $yP {$M . 0} = $distance;
  63	} else {	# Past joint at point F
  64		$extra = $distance - $yF;
  65		$xJ {$M . 3} = $xP {$M . 0} = Interpolate ($extra, $EF, $xF, $xE);
  66		$yJ {$M . 3} = $yP {$M . 0} = Interpolate ($extra, $EF, $yF, $yE);
  67	}
  68
  69	# Calculate both joints of meridians
  70
  71	if ($m == 0) {	# Straight, centre of octant meridian
  72		# This is a good place to calculate a bunch of things for centre line, AG.
  73
  74		# No joints; give coordinates of point G for the joints
  75		@xJ {$M . 1, $M . 2} = ($xG, $xG);
  76		@yJ {$M . 1, $M . 2} = ($yG, $yG);
  77		# Lenght of meridian segments; this whole meridian is a single segment
  78		@L {$M . 0, $M . 1, $M . 2, $M . 3, $M . 4} = ($AG, 0, 0, $AG, $AG);
  79		# May as well calculate parallel intersections along this straight line
  80		$dP [$m] = $L{$M . 4} / 90;
  81		foreach $p (1 .. 89) {
  82			$xP {$M . $p} = $xG - $p * $dP[$m];
  83			$yP {$M . $p} = $yG;
  84		}
  85	} else {	# Meridians other than 0�.		
  86		# Calculate joints by line intersection from A and M
  87		($flag, $xJ {$M . 1}, $yJ {$M . 1}) = 
  88			LineIntersection ($xM, $yM, 2*$m/3, $xA, $yA, $m);
  89		if ($flag ne "Intersect") { die "$flag while calculating 1st joint at $m.\n"; }
  90		($flag, $xJ {$M . 2}, $yJ {$M . 2}) = 
  91			LineIntersection ($xM, $yM, 2*$m/3, $xJ {$M . 3}, $yJ {$M . 3}, $m/3);
  92		if ($flag ne "Intersect") { die "$flag while calculating 2nd joint at $m.\n"; }
  93		# Calculate length of meridian segments: These are calculated from point A 
  94		# to equator, even if this meridian starts at parallel 85�.
  95		$L {$M . 0} = Length ($xA, $yA, $xJ {$M . 1}, $yJ {$M . 1});
  96		$L {$M . 1} = Length ($xJ {$M . 1}, $yJ {$M . 1}, $xJ {$M . 2}, $yJ {$M . 2});
  97		$L {$M . 2} = Length ($xJ {$M . 2}, $yJ {$M . 2}, $xJ {$M . 3}, $yJ {$M . 3});
  98		$L {$M . 3} = $L {$M . 1} + $L {$M . 2};	# Torrid + Temperate
  99		$L {$M . 4} = $L {$M . 3} + $L {$M . 0};	# Total meridian length from A
 100		# Calculate parallel intersections
 101		$dP[$m] = $L {$M . 4} / 90;	# Distance between parallels
 102		foreach $p (1 .. 89) {
 103			$distance = $p * $dP [$m];
 104			if ($distance <= $L {$M . 2} ) {		# Torrid segment
 105				$xP{$M.$p} = 
 106					Interpolate ($distance, $L{$M.2}, $xJ{$M.3}, $xJ{$M.2});
 107				$yP{$M.$p} = 
 108					Interpolate ($distance, $L{$M.2}, $yJ{$M.3}, $yJ{$M.2});
 109			} elsif ($distance <= $L {$M . 3} ) {	# Temperate seg., past one joint
 110				$extra = $distance - $L {$M . 2};
 111				$xP{$M.$p} = 
 112					Interpolate ($extra, $L{$M.1}, $xJ{$M.2}, $xJ{$M.1});
 113				$yP{$M.$p} = 
 114					Interpolate ($extra, $L{$M.1}, $yJ{$M.2}, $yJ{$M.1});
 115			} else {					# Frigid segment, past two joints
 116				$extra = $distance - $L {$M . 3};
 117				$xP{$M.$p} = Interpolate ($extra, $L{$M.0}, $xJ{$M.1}, $xA);
 118				$yP{$M.$p} = Interpolate ($extra, $L{$M.0}, $yJ{$M.1}, $yA);
 119			}
 120		}	# End of foreach $p loop, that is, end of parallel calculation
 121	}	# End of calculating meridian joints
 122}	# End of loop for each meridian.
 123
 124#	Start point of every meridian was previously set to be point A. However, the drawn
 125#	line for the minor meridians will start at parallel 85�; this could not be set previously
 126#	because that parallel still hadn't been calculated.
 127foreach $m (0 .. 45) {
 128	$M = $m . ',' ;
 129	unless ($m % 5 == 0) {
 130		$xJ {$M . 0} = $xP {$M . 85};
 131		$yJ {$M . 0} = $yP {$M . 85};
 132	}
 133}
 134
 135# Output all the values.
 136$PointFormat = "%3s \t%11.4f \t%11.4f\n";
 137$LengthFormat = "%s \t%11.4f\n";
 138$HashKeyValue = "%4s \t%11.4f\n";
 139
 140$Skip = "Yes";
 141unless ($Skip) {	# Option to skip printing hashes in spreadsheet readable csv format
 142	open (CSV, ">Macros4\/Hashes4.csv");
 143	print CSV "\ndP\n";		# @dP is an array, not a hash; can't use sub PrintHash
 144	foreach $i (0 .. 45) { printf CSV "\t%10d",$i; } print CSV "\n";
 145	print CSV " ";
 146	foreach $m (0 .. 45) {
 147		if (defined ($dP[$m]) ) {printf CSV "\t%11.4f",$dP[$m];
 148		} else { printf CSV "\t undef"; }
 149	}
 150	print CSV "\n";
 151	print CSV "\nL\n"; PrintHash (4,%L);
 152	print CSV "\nxJ\n"; PrintHash (3,%xJ);
 153	print CSV "\nyJ\n"; PrintHash (3,%yJ);
 154	print CSV "\nxP\n"; PrintHash (89,%xP);
 155	print CSV "\nyP\n"; PrintHash (89,%yP);
 156	close (CSV);
 157	print 'Wrote values of arrays to file "Macros4/Hashes4.csv".', "\n";
 158}	# End of Skip printing hashes in Spreadsheet readable csv format
 159
 160$Skip = "Yes";
 161unless ($Skip) {	# Option to skip printing, hashes by their keys
 162	print "\n* * * Keys of hashes\n\n";
 163	print "* * * \%xJ\n";
 164	foreach $key (keys (%xJ)) {
 165		printf ($HashKeyValue, $key, $xJ{$key});
 166	}
 167	print "* * * \%yJ\n";
 168	foreach $key (keys (%yJ)) {
 169		printf ($HashKeyValue, $key, $yJ{$key});
 170	}
 171	print "* * * \%L\n";
 172	foreach $key (keys (%L)) {
 173		printf ($HashKeyValue, $key, $L{$key});
 174	}
 175	print "\n* * * COORDINATES\n";
 176	print "* * * \%xP\n";
 177	foreach $key (keys (%xP)) {
 178		printf ($HashKeyValue, $key, $xP{$key});
 179	}
 180	print "* * * \%yP\n";
 181	foreach $key (keys (%yP)) {
 182		printf ($HashKeyValue, $key, $yP{$key});
 183	}
 184}	# End of Skip of printing hashes by their keys
 185
 186$Skip = "Yes";
 187unless ($Skip) {	# Option to skip printing points and lengths
 188print "* * * Major points (point, x, y);\n";
 189printf ($PointFormat, "M", $xM, $yM);
 190printf ($PointFormat, "N", $xN, $yN);
 191printf ($PointFormat, "A", $xA, $yA);
 192printf ($PointFormat, "B", $xB, $yB);
 193printf ($PointFormat, "D", $xD, $yD);
 194printf ($PointFormat, "E", $xE, $yE);
 195printf ($PointFormat, "F", $xF, $yF);
 196printf ($PointFormat, "G", $xG, $yG);
 197# Main lengths
 198print "\n* * * Main lengths\n";
 199printf ($LengthFormat,'MG',$xG);
 200printf ($LengthFormat,'MN',$MN);
 201printf ($LengthFormat,'MA',$MA);
 202printf ($LengthFormat,'MB',$MB);
 203printf ($LengthFormat,'AG',$AG);
 204printf ($LengthFormat,'AB',$AB);
 205printf ($LengthFormat,'BD',$BD);
 206printf ($LengthFormat,'EF',$EF);
 207printf ($LengthFormat,'FG',$FG);
 208printf ($LengthFormat,'EFG',$EFG);
 209printf ($LengthFormat,'ABDE',2 * $EFG);
 210print "\n";
 211printf ($LengthFormat,'Frigid parallel spacing',$DeltaFrigid);
 212printf ($LengthFormat,'Equatorial meridian spacing',$DeltaMEq);
 213}	# End of Skip of printing major points and lengths
 214
 215# Make TurboCAD macros.
 216# Scaffold triangle macro
 217
 218$Skip = "Yes";
 219unless ($Skip) {	# Option to skip making macro files of half octant or one octant
 220StartMacroFile (MACRO, "Macros4\/Triangle.txt");
 221print MACRO ("Pen 9\r");
 222print MACRO ("Line ", Coordinate($xM, $yM), Coordinate ($xN, $yN), "\\\r");
 223print MACRO (Coordinate($xG, $yG), Coordinate ($xM, $yM), '[`;`]', "\r");
 224EndMacroFile (MACRO);
 225print 'Wrote file "Macros4/Triangle.txt".', "\n";
 226
 227# And the same thing in Gene's coordinate system
 228StartMacroFile (MACRO, "Macros4\/TriangleG.txt");
 229print MACRO ("Pen 9\r");
 230print MACRO ("Line ", Coordinate(MJtoG($xM,$yM)), Coordinate(MJtoG($xN,$yN)), "\\\r");
 231print MACRO (Coordinate(MJtoG($xG,$yG)), Coordinate(MJtoG($xM,$yM)), '[`;`]', "\r");
 232EndMacroFile (MACRO);
 233print 'Wrote file "Macros4/TriangleG.txt."', "\n";
 234
 235# Half-octant outline macro
 236StartMacroFile (MACRO, "Macros4\/Half-Octant.txt");
 237print MACRO ("Pen 12\r");
 238print MACRO ("Line ", Coordinate($xA, $yA), Coordinate ($xB, $yB), "\\\r");
 239print MACRO (Coordinate($xD, $yD), Coordinate ($xE, $yE), "\\\r");
 240print MACRO (Coordinate($xF, $yF), Coordinate ($xG, $yG), '[`;`]', "\r");
 241EndMacroFile (MACRO);
 242print 'Wrote file "Macros4/Half-Octant.txt".', "\n";
 243
 244# And the same thing in Gene's coordinate system
 245StartMacroFile (MACRO, "Macros4\/Half-OctantG.txt");
 246print MACRO ("Pen 12\r");
 247print MACRO ("Line ", Coordinate(MJtoG($xA,$yA)), Coordinate(MJtoG($xB, $yB)), "\\\r");
 248print MACRO (Coordinate(MJtoG($xD,$yD)), Coordinate(MJtoG($xE, $yE)), "\\\r");
 249print MACRO (Coordinate(MJtoG($xF,$yF)), Coordinate(MJtoG($xG, $yG)), '[`;`]', "\r");
 250EndMacroFile (MACRO);
 251print 'Wrote file "Macros4/Half-OctantG.txt".', "\n";
 252
 253# Macro of 8 equilateral triangles
 254StartMacroFile (MACRO, "Macros4\/MScaffold.txt");
 255print MACRO ("Pen 9\r");
 256MLinesMacro (MACRO,0,$xN,-$yN,1,$xM,$yM,1,$xN,$yN,1,$xN,-$yN,1,$xM,$yM,6,
 257	$xN,$yN,1,$xN,$yN,2,$xM,$yM,1);
 258MLinesMacro (MACRO,0,$xN,$yN,1,$xM,$yM,7,$xN,$yN,2,$xN,$yN,3,$xM,$yM,7);
 259MLinesMacro (MACRO,0,$xN,$yN,2,$xM,$yM,3,$xN,$yN,3,$xN,$yN,4,$xM,$yM,3);
 260MLinesMacro (MACRO,0,$xN,$yN,3,$xM,$yM,5,$xN,$yN,4);
 261EndMacroFile (MACRO);
 262print 'Wrote file "Macros4/MScaffold.txt".', "\n";
 263
 264# Macro of 8 Octants
 265StartMacroFile (MACRO, "Macros4\/MOctants.txt");
 266print MACRO ("Pen 12\r");
 267MLinesMacro (MACRO,6,	$xF,$yF,$xF,-$yF,$xE,-$yE,$xD,-$yD,$xB,-$yB,
 268				$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,$xF,$yF);
 269MLinesMacro (MACRO,1,	$xF,-$yF,$xE,-$yE,$xD,-$yD,$xB,-$yB,
 270				$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,$xF,$yF);
 271MLinesMacro (MACRO,2,	$xB,-$yB,$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,
 272				$xF,$yF,$xF,-$yF,$xE,-$yE,$xD,-$yD);
 273MLinesMacro (MACRO,7,	$xF,-$yF,$xE,-$yE,$xD,-$yD,$xB,-$yB,
 274				$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,$xF,$yF);
 275MLinesMacro (MACRO,8,	$xB,-$yB,$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,
 276				$xF,$yF,$xF,-$yF,$xE,-$yE,$xD,-$yD);
 277MLinesMacro (MACRO,3,	$xF,-$yF,$xE,-$yE,$xD,-$yD,$xB,-$yB,
 278				$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,$xF,$yF);
 279MLinesMacro (MACRO,5,	$xF,-$yF,$xE,-$yE,$xD,-$yD,$xB,-$yB,
 280				$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,$xF,$yF);
 281MLinesMacro (MACRO,4,	$xB,-$yB,$xA,$yA,$xB,$yB,$xD,$yD,$xE,$yE,
 282				$xF,$yF,$xF,-$yF,$xE,-$yE,$xD,-$yD);
 283EndMacroFile (MACRO);
 284print 'Wrote file "Macros4/MOctants.txt".', "\n";
 285
 286# Meridians macros
 287
 288StartMacroFile (MAJOR, "Macros4\/MajorMeridians.txt");
 289print MAJOR ("Pen 12\r");
 290StartMacroFile (MINOR, "Macros4\/MinorMeridians.txt");
 291print MINOR ("Pen 13\r");
 292MeridianMacro ();
 293EndMacroFile (MAJOR);
 294EndMacroFile (MINOR);
 295print 'Wrote files "Macros4/MajorMeridians.txt" and "Macros4/MinorMeridians.txt".', "\n";
 296
 297# Parallels macros
 298StartMacroFile (MAJOR, "Macros4\/ParallelsMajor.txt");
 299print MAJOR ("Pen 12\r");
 300StartMacroFile (MINOR, "Macros4\/ParallelsMinor.txt");
 301print MINOR ("Pen 13\r");
 302ParallelMacro (1,89,0,45);
 303EndMacroFile (MAJOR);
 304EndMacroFile (MINOR);
 305print 'Wrote files "Macros4/ParallelsMajor.txt" and "Macros4/ParallelsMinor.txt".', "\n";
 306}	# End of Skip of making macro files of half octant or one octant
 307
 308$Skip = "Yes";
 309unless ($Skip) {	# Option to skip making macro files of parallel points and of joints points
 310# Macro of points for the intersections of the parallels with the meridians; this will be
 311# a large file which will take quite awhile to plot; pens to be used are specified in the 
 312# subroutine. Meant only to check if the data seems right.
 313StartMacroFile (MACRO, "Macros4\/Points.txt");
 314PointsMacro ();
 315EndMacroFile (MACRO);
 316print 'Wrote file "Macros4/Points.txt".', "\n";
 317
 318# Macro of meridian joints and end points (4 points per meridian)
 319StartMacroFile (MACRO, "Macros4\/Joints.txt");
 320JointsMacro ();
 321EndMacroFile (MACRO);
 322print 'Wrote file "Macros4/Joints.txt".', "\n";
 323}	# End of Skip of making macro files
 324
 325# Macros for all eight whole octants
 326
 327$Skip = "Yes";
 328unless ($Skip) {	# Option to skip making macro files of eight octants
 329# Make arrays for Octant; I should have done this much earlier in the program
 330@xOct = ($xA,$xB,$xD,$xE,$xF,$xG);
 331@yOct = ($yA,$yB,$yD,$yE,$yF,$yG);
 332
 333# Make a copy of arrays and hashes that I will want to convert; the other ones
 334# will be changed to coordinates for different octant halves.
 335@xOctC = @xOct;	@yOctC = @yOct;
 336%xJC = %xJ;		%yJC = %yJ;
 337%xPC = %xP;	%yPC = %yP;
 338
 339# Make grid files
 340StartMacroFile (MINOR, "Macros4\/M-Grid-Minor.txt");
 341print MINOR ("Pen 15\r");
 342GridMacro (MINOR, -20000,20000,-10000,10000,100,1000);
 343EndMacroFile (MINOR);
 344
 345StartMacroFile (MAJOR, "Macros4\/M-Grid-Major.txt");
 346print MAJOR ("Pen 14\r");
 347GridMacro (MAJOR, -20000,20000,-10000,10000,1000,0);
 348EndMacroFile (MAJOR);
 349
 350StartMacroFile (MAJOR, "Macros4\/M-Grid-5Kred.txt");
 351print MAJOR ("Pen 10\r");
 352GridMacro (MAJOR, -20000,20000,-10000,10000,5000,0);
 353EndMacroFile (MAJOR);
 354
 355print 'Wrote files "M-Grid-Minor.txt", "M-Grid-Major.txt" and "M-Grid-5Kred.txt" in "Macros4/".',
 356	"\n";
 357# Draw octants, half at a time
 358StartMacroFile (MACRO, "Macros4\/M-Octants.txt");
 359print MACRO ("Pen 11\r");
 360foreach $Oct (1 .. 8) {		# For octant 1 and octant 2
 361	foreach $i (0..5) {		# For each point (A,B,D,E,F,G)
 362		($xOct[$i],$yOct[$i]) = MJtoG ($xOctC[$i],-$yOctC[$i],$Oct);	# West half
 363	}
 364	OctantMacro ();
 365	foreach $i (0..5) {
 366		($xOct[$i],$yOct[$i]) = MJtoG ($xOctC[$i],$yOctC[$i],$Oct);	# East half
 367	}
 368	OctantMacro ();
 369}
 370EndMacroFile (MACRO);
 371print 'Wrote file "Macros4/M-Octants.txt".', "\n";
 372
 373# Draw Meridians, half an octant at a time
 374StartMacroFile (MINOR, "Macros4\/M-Meridians-Minor.txt");
 375StartMacroFile (MAJOR, "Macros4\/M-Meridians-Major.txt");
 376print MINOR ("Pen 13\r");
 377print MAJOR ("Pen 12\r");
 378
 379foreach $Oct (1 .. 8) {		# For each octant
 380	# West half of octant
 381	foreach $m (0..45) {
 382		$M = $m . ',' ;
 383		foreach $i (0..3) {
 384			($xJ{$M.$i},$yJ{$M.$i}) = MJtoG ($xJC{$M.$i},-$yJC{$M.$i},$Oct);
 385		}
 386	}
 387	MeridianMacro ();
 388
 389	# East half of octant
 390	foreach $m (0..45) {
 391		$M = $m . ',' ;
 392		foreach $i (0..3) {
 393			($xJ{$M.$i},$yJ{$M.$i}) = MJtoG ($xJC{$M.$i},$yJC{$M.$i},$Oct);
 394		}
 395	}
 396	MeridianMacro ();
 397}
 398EndMacroFile (MINOR);
 399EndMacroFile (MAJOR);
 400print 'Wrote files "Macros4/M-Meridians-Minor.txt" and "Macros4/M-Meridians-Major.txt".', "\n";
 401
 402# Draw Parallels, half an octant at a time
 403StartMacroFile (MINOR, "Macros4\/M-Parallels-Minor.txt");
 404StartMacroFile (MAJOR, "Macros4\/M-Parallels-Major.txt");
 405print MINOR ("Pen 13\r");
 406print MAJOR ("Pen 12\r");
 407
 408foreach $Oct (1 .. 8) {		# For each octant
 409	# West half of octant
 410	foreach $m (0..45) {	# Convert coordinates
 411		$M = $m . ',' ;
 412		foreach $p (0..89) {
 413			($xP{$M.$p},$yP{$M.$p}) = MJtoG ($xPC{$M.$p},-$yPC{$M.$p},$Oct);
 414		}
 415	}
 416	# All parallels from 1� to 89� of latitude; skip parallel 0� &#x2013; octant boundary
 417	ParallelMacro (1, 89, 0, 45, ", Octant $Oct West");
 418
 419	# East half of octant
 420	foreach $m (0..45) {	# Convert coordinates
 421		$M = $m . ',' ;
 422		foreach $p (0..89) {
 423			($xP{$M.$p},$yP{$M.$p}) = MJtoG ($xPC{$M.$p},$yPC{$M.$p},$Oct);
 424		}
 425	}
 426	# All parallels from 1� to 89� of latitude; skip parallel 0� &#x2013; octant boundary
 427	ParallelMacro (1, 89, 0, 45, ", Octant $Oct East");
 428}
 429EndMacroFile (MINOR);
 430EndMacroFile (MAJOR);
 431print 'Wrote file "Macros4/M-Parallels-Minor.txt" and "Macros4/M-Parallels-Major.txt".', "\n";
 432
 433# Draw Parallels, half an octant at a time, but connecting only nodes of major meridians;
 434# that is, not connecting nodes of minor meridians
 435StartMacroFile (MINOR, "Macros4\/M-Parallels-Minor5.txt");
 436StartMacroFile (MAJOR, "Macros4\/M-Parallels-Major5.txt");
 437print MINOR ("Pen 13\r");
 438print MAJOR ("Pen 12\r");
 439
 440foreach $Oct (1 .. 8) {		# For each octant
 441	# West half of octant
 442	foreach $m (0..45) {	# Convert coordinates
 443		$M = $m . ',' ;
 444		foreach $p (0..89) {
 445			($xP{$M.$p},$yP{$M.$p}) = MJtoG ($xPC{$M.$p},-$yPC{$M.$p},$Oct);
 446		}
 447	}
 448	# All parallels from 1� to 89� of latitude; skip parallel 0� &#x2013; octant boundary
 449	ParallelMacro5 (1, 89, 0, 45, ", Octant $Oct West");
 450
 451	# East half of octant
 452	foreach $m (0..45) {	# Convert coordinates
 453		$M = $m . ',' ;
 454		foreach $p (0..89) {
 455			($xP{$M.$p},$yP{$M.$p}) = MJtoG ($xPC{$M.$p},$yPC{$M.$p},$Oct);
 456		}
 457	}
 458	# All parallels from 1� to 89� of latitude; skip parallel 0� &#x2013; octant boundary
 459	ParallelMacro5 (1, 89, 0, 45, ", Octant $Oct East");
 460}
 461EndMacroFile (MINOR);
 462EndMacroFile (MAJOR);
 463print 'Wrote file "Macros4/M-Parallels-Major5.txt" and "Macros4/M-Parallels-Minor5.txt".', "\n";
 464
 465}	# End of skip of making macro files of eight octants
 466
 467# Read some Coastal Data, convert to M-map coordinates, and output a TurboCAD macro.
 468# File read is of MAPGEN data format: two columns ASCII flat file with: longitude tab latitude;
 469# at the start of each segment there is a line containing only "# -b".
 470# This particular file covers longitude from -67� to -59� and latitude from 43� to 48�, that is,
 471# Nova Scotia's area. File is 5781.dat.gz, is zipped, and has 92,227 lines.
 472
 473$Skip = "Yes";
 474unless ($Skip) {	# Option to skip making macro files of Coastal Data
 475print "Going to read land data.\n";
 476$maxX=-99999;	$maxY=-99999;	$minX=999999;	$minY=999999;
 477StartMacroFile (LAND, "Macros4\/NSmacro.txt");
 478print LAND ("Pen 0\r");
 479open (DATA, "zcat 5781.dat.gz | ");
 480$nData = 0;
 481while (<DATA>) {
 482	$Line = $_ ;
 483	chomp($Line);
 484	$nData ++;
 485	if ($nData % 1000 == 0) { print "Read $nData lines of land data so far.\n"; }
 486	if ($Line ne "\# -b") {
 487		($Long,$Lat) = split(/\t/,$Line);
 488		# Convert to M-map x and y coordinates
 489		($x,$y) = Real2Map ($Long,$Lat);
 490		# Keep track of maximum and minimum values
 491		if ($x < $minX) {$minX = $x;}
 492		if ($x > $maxX) {$maxX = $x;}
 493		if ($y < $minY) {$minY = $y;}
 494		if ($y > $maxY) {$maxY = $y}
 495		# Make arrays of points for this segment to be used by LandMacro
 496		push (@Xs,$x);
 497		push (@Ys,$y);
 498	} elsif (defined(@Xs) ) {	# Do only if data points have already been read
 499		LandMacro ();	# Draw the segment of boundary with arrays @Xs and @Ys
 500		@Xs = ();	@Ys = ();	# Empty arrays, ready for next segment
 501	}
 502}
 503close (DATA);
 504if ($Line ne "\# -b") { # Draw last segment, if it hasn't been drawn
 505	LandMacro ();
 506}
 507EndMacroFile (LAND);
 508print "Read a total of $nData lines, and wrote file",' "Macros4/NSmacro.txt".', "\n";
 509print $minX, ' <= x <= ',$maxX," and ",$minY, ' <= y <= ', $maxY, "\n";
 510
 511}	# End of skipping making macro for Coastal Data
 512
 513$Skip = "";
 514unless ($Skip) {	# Option to skip making macro file of NS 1:2M of Coastal Data
 515print "Going to read land data.\n";
 516$maxX=-99999;	$maxY=-99999;	$minX=999999;	$minY=999999;
 517StartMacroFile (LAND, "Macros4\/NS-2M-macro4.txt");
 518print LAND ("Pen 0\r");
 519open (DATA, "zcat NS-2M.dat.gz | ");
 520$nData = 0;
 521while (<DATA>) {
 522	$Line = $_ ;
 523	chomp($Line);
 524	$nData ++;
 525	if ($nData % 1000 == 0) { print "Read $nData lines of land data so far.\n"; }
 526	if ($Line ne "\# -b") {
 527		($Long,$Lat) = split(/\t/,$Line);
 528		# Convert to M-map x and y coordinates
 529		($x,$y) = Real2Map ($Long,$Lat);
 530		# Keep track of maximum and minimum values
 531		if ($x < $minX) {$minX = $x;}
 532		if ($x > $maxX) {$maxX = $x;}
 533		if ($y < $minY) {$minY = $y;}
 534		if ($y > $maxY) {$maxY = $y}
 535		# Make arrays of points for this segment to be used by LandMacro
 536		push (@Xs,$x);
 537		push (@Ys,$y);
 538	} elsif (defined(@Xs) ) {	# Do only if data points have already been read
 539		# Shift coordinates to see if TurboCAD plots them better
 540		$dx = -7600;
 541		$dy = 3700;
 542		$n = @Xs -1;
 543		foreach $i (0..$n) {
 544			$Xs[$i] = $Xs[$i] - $dx;
 545			$Ys[$i] = $Ys[$i] - $dy;
 546		}
 547		LandMacro ();	# Draw the segment of boundary with arrays @Xs and @Ys
 548		@Xs = ();	@Ys = ();	# Empty arrays, ready for next segment
 549	}
 550}
 551close (DATA);
 552if ($Line ne "\# -b") { # Draw last segment, if it hasn't been drawn
 553	LandMacro ();
 554}
 555EndMacroFile (LAND);
 556print "Read a total of $nData lines, and wrote file",' "Macros4/NS-2M4-macro .txt".', "\n";
 557print $minX, ' <= x <= ',$maxX," and ",$minY, ' <= y <= ', $maxY, "\n";
 558print "Need plot area ", $maxX-$dx, " by ", $maxY-$dy,".\n";
 559
 560}	# End of skipping making macro for NS, 2:M Coastal Data
 561
 562
 563# -  -  -  -  -  -  -  -  -    SUBROUTINES         -  -  -  -  -  -  -  -  -  -  -  -  -
 564sub PrintHash {
 565	# Prints one of my two-dimension arrays in the form of a hash; it is assumed that
 566	# the indices are separated by a comma, the first indece varies from 0 to 45 and the
 567	# second indice varies from 0 to the first input parameter.
 568	# Values are printed to filehandle CSV.
 569	my ($nI, %Hash) = @_;
 570	my ($i,$m);
 571	foreach $i (0 .. 45) { printf CSV "\t%10d",$i; } print CSV "\n";
 572	foreach $i (0 .. $nI) {
 573		print CSV $i;
 574		foreach $m (0 .. 45) {
 575			if (defined ($Hash{$m.",".$i})) {printf CSV "\t%11.4f",$Hash{$m.",".$i};
 576			} else { printf CSV "\t undef"; }
 577		}
 578		print CSV "\n";
 579	}
 580}	# End of subroutine Hash
 581
 582sub LineIntersection {	# Written 2010-02-28
 583# Subroutine/function to calculate coordinates of point of intersection of two lines which
 584# are given by a point on the line and the line's slope angle in degrees.
 585#
 586# Return is an array of either one or three values:
 587#	- if the lines don't intersect or are the same, the return is "Parallel";
 588#	- if the lines intersect, the returns are: "Intersect" to mean that it worked, and 
 589#	the x and y coordinates of the point of intersection.
 590#	- if not enough arguments were given or the angles were too large, then return is
 591#	 "Error".
 592#
 593# Arguments should be 6, in this order:
 594#	x and y coordinates of point of first line; slope of first line in degrees;
 595#	x and y coordinates of point of second line; slope of second line in degrees;
 596#
 597# Equations used are from: slope of line = tangent angle = delta-y / delta-x, and the fact
 598#	 that intersection point x,y is on both lines.
 599
 600use Math::Trig;
 601my ($nArguments,$xp,$yp,$m1,$m2);
 602$nArguments=@_ ;
 603my ($x1,$y1,$angle1,$x2,$y2,$angle2) = @_ ;
 604
 605if ($nArguments != 6) {	# Not the exact number of arguments
 606	print "Sub LineIntersection requires 6 arguments but received $nArguments \n";
 607	return "Error";
 608} elsif ($angle1<-180 || $angle2<-180 || $angle1>180 || $angle2>180) { # Angles too large
 609	print "Sub LineIntersection: Angles should be between -180 and 180. \n";
 610	return "Error";
 611} elsif ($angle1==$angle2 || $angle1 == ($angle2+180) || $angle1 == ($angle2-180) ) {
 612	return "Parallel";		# Lines are either parallel or the same line
 613} elsif (	not defined($x1) || not defined($y1) || not defined($angle1) || 
 614		not defined($x2) || not defined($y2) || not defined($angle2) ) {
 615	print ("Undefined in LineInterpolation: $x1,$y1,$angle1,$x2,$y2,$angle2\n");
 616	return "Error";
 617} else {	# Lines intersect. Calculate point of intersection.
 618	if ($angle1 == 0 || $angle1 == 180 || $angle1 == -180) {		# Line 1 horizontal
 619		if ($angle2 == 90 || $angle2 == -90) {	# Line 1 horizontal, line 2 vertical
 620			return ("Intersect",$y1,$x2);
 621		} else {					# Line 1 horizontal, line 2 slanted
 622			$xp = $x2 - ($y2 - $y1) / tan(deg2rad($angle2));
 623			return ("Intersect",$xp,$y1);
 624		}
 625	} elsif ($angle1 == 90 || $angle1 == -90) {	# Calculate when line 1 is vertical
 626		if ($angle2 == 0 || $angle2 == 180 || $angle2 == -180) {
 627								# Line 1 vertical, line 2 horizontal
 628			return ("Intersect",$x1,$y2);
 629		} else {					# Line 1 vertical, line 2 slanted
 630			$yp = $y2 - ($x2 - $x1) * tan(deg2rad($angle2));
 631			return ("Intersect",$x1,$yp);
 632		}
 633	} elsif ($angle2 == 0 || $angle2 == 180 || $angle2 == -180) {
 634								# Line 1 slanted, line 2 horizontal
 635			$xp = $x1 - ($y1 - $y2) / tan(deg2rad($angle1));
 636			return ("Intersect",$xp,$y2);
 637	} elsif ($angle2 == 90 or $angle2 == -90) {	# Line 1 slanted, line 2 vertical
 638			$yp = $y1 - ($x1-$x2) * tan(deg2rad($angle1));
 639			return ("Intersect",$x2,$yp);
 640	} else {						# Line 1 slanted, line 2 slanted
 641			$m1 = tan(deg2rad($angle1));
 642			$m2 = tan(deg2rad($angle2));
 643			$xp = ($m1 * $x1 - $m2 * $x2 - $y1 + $y2) / ($m1 - $m2);
 644			$yp = $m1 * $xp - $m1 * $x1 + $y1;
 645			return ("Intersect",$xp,$yp);
 646	}	#End of all the vertical, horizontal, slanted combinations
 647}	# End of if-else &#x2013; lines intersect
 648} # End of Subroutine LineIntersection
 649
 650sub Length {
 651	# Input are x1,y1,x2,y2
 652	my ($x1,$y1,$x2,$y2) = @_;
 653	return sqrt( ($x1-$x2)**2 + ($y1-$y2)**2);
 654}
 655
 656sub Interpolate {
 657# Inputs are 4: length wanted, of total-segment-length, start, end;
 658# Total-segment-length is different from (end - start); end and start may be x-coordinates,
 659# or y-coordinates, while length takes into account the other coordinates.
 660# (Start - End ) / Length = (Wanted - Start) / NewLength
 661# Returns single value: Wanted.
 662	my ($NewLength, $Length, $Start, $End) = @_;
 663	my ($Wanted);
 664	$Wanted = $Start + ($End - $Start) * $NewLength / $Length;
 665$Skip = "Yes";
 666unless ($Skip) {	# Skip of not the printing of arguments and value being returned
 667	my ($nIn, $i);
 668	$nIn = @_ - 1;
 669#	print "Interpolate inputs are ", $nIn+1, ": ";
 670	foreach $i (0 .. $nIn) {
 671		if (defined ($_[$i]) ) {
 672			printf ("%11.4f, ", $_[$i]);
 673		} else {
 674			printf ("%11s, ", 'UNDEFINED!');
 675		}
 676	}
 677	if (defined ($Wanted) ) {
 678#		print "Returning: ";
 679		printf ("\>%11.4f",$Wanted);
 680	} else {
 681		printf ("\>%11s", 'UNDEFINED!');
 682	}
 683	print "\n";
 684}	# End of Skip printing arguments and value being returned
 685	return $Wanted;
 686}	# End of subroutine Interpolate
 687
 688# -  -  -  -  -  -  -  -  -    SUBROUTINES relating to TurboCAD    -  -  -  -  -  -  -  -  -  -  -  -  -
 689
 690sub StartMacroFile {
 691	my ($FileHandle, $ToOpen);
 692	$FileHandle = $_ [0];
 693	$ToOpen = '>' . $_[1];
 694	open ($FileHandle, $ToOpen);
 695	print $FileHandle ('#PenWas = %WT.P', "\r");
 696}
 697
 698sub EndMacroFile {
 699	my ($FileHandle);
 700	$FileHandle = $_ [0];
 701	print $FileHandle ('Pen #PenWas', "\r");
 702	print $FileHandle ("Bell\rBell\r");
 703	close ($FileHandle);
 704}
 705
 706sub Coordinate {
 707	# returns string with coordinate format for TurboCAD;
 708	# expects two arguments: x and y.
 709	return '[a,' . sprintf ("%11.4f,%11.4f",$_[0],$_[1]) . ',0,` `] '
 710}
 711
 712sub PointsMacro {
 713	# I will be lazy and use global variables instead of trying to pass hashes as 
 714	# parameters.
 715	# Plotting all parallel intersections, that have already been calculated: arrays
 716	# %xP, %yP.
 717	# Subroutine will prepare TurboCad commands to plot each point. This will be written
 718	# to filehandle MACRO
 719	my ($m, $p, $i, $s, $c);
 720	$s = 'PT Period Pen = #';
 721	print MACRO ("\#MiMi = 13\r\#MjMi = 12\r\#MjMj = 9\r");	# Pen variables
 722	print MACRO ($s, 'MjMj ', Coordinate($xA, $yA), '[`;`]', "\r" );	# Point at the pole
 723	foreach $m (0 .. 45) {
 724		foreach $p (0 .. 89) {
 725			$i = $m . "," . $p ;
 726			if (defined ($xP {$i}) && defined ($yP {$i})) {
 727				if ($m % 5 == 0 && $p % 5 == 0) {
 728					$c = 'MjMj';		# major meridian and major parallel
 729				} elsif ($m % 5 == 0 || $p % 5 == 0) {
 730					$c = 'MjMi';	# major parallel or major meridian
 731				} else {
 732					$c = 'MiMi';	# minor parallel with minor meridian
 733				}
 734				print MACRO ($s,$c,' ', Coordinate($xP{$i},$yP{$i}), '[`;`]', "\r" );
 735			}
 736		}
 737	}
 738}		# End of subroutine PointsMacro.
 739
 740sub JointsMacro {
 741	# I will be lazy and use global variables instead of trying to pass hashes as 
 742	# parameters.
 743	# Plotting all joint intersections, that have already been calculated: arrays
 744	# %xJ, %yJ.
 745	# Subroutine will prepare TurboCad commands to plot each point. This will be written
 746	# to filehandle MACRO.
 747	# This subroutine is almost the same as PointsMacro; differences are:  different
 748	# arrays and array indeces; not plotting "Point at the pole"; only two pens used;
 749	# and there is no check for whether the variable is defined.
 750
 751	my ($m, $p, $i, $s, $c);
 752	$s = 'PT Period Pen = #';
 753	print MACRO ("\#Major = 9\r\#Minor = 12\r");	# Pen variables
 754	foreach $m (0 .. 45) {
 755		if ($m % 5 == 0) {
 756			$c = 'Major';		# Major meridian
 757		} else {
 758			$c = 'Minor';		# Minor meridian
 759		}
 760		foreach $p (0 .. 3) {
 761			$i = $m . "," . $p ;
 762			print MACRO ($s,$c,' ', Coordinate($xJ{$i},$yJ{$i}), '[`;`]', "\r" );
 763		}
 764	}
 765}		# End of subroutine JointsMacro.
 766
 767
 768sub MeridianMacro {
 769	# I will be lazy and use global variables instead of trying to pass hashes as 
 770	# parameters.
 771	# It is assumed that every meridian multiple of 5 will be written to file MAJOR;
 772	# other meridians are written to file MINOR; Meridian numbers are from 0 to 45;
 773	# will use hashes %xJ and %yJ, each of which is indexed as "i,j", where i= 0 .. 45,
 774	# and j = 0 .. 4, for start polar point, two joints and equator point.
 775	my ($m, $M, $s, $FileHandle);
 776	foreach $m (0 .. 45) {
 777		$M = $m . ',' ;
 778		if ( ($m % 5) == 0 ) {		# Major meridian
 779			$FileHandle = MAJOR;
 780		} else {				# Minor meridian
 781			$FileHandle = MINOR;
 782		}
 783		print $FileHandle ('LI ', Coordinate($xJ {$M . 0}, $yJ {$M . 0})); # First point
 784		for ($s = 1; $s <= 3; $s++ ) {	# for each line segment
 785			print $FileHandle ("\\\r",Coordinate($xJ {$M . $s},$yJ {$M . $s}) );
 786		}	# End of doing one meridian with 3 segments
 787		print $FileHandle ('[`;`]', "\r");		# End drawing of this multi-line
 788	}	# End of for(m)-loop
 789}	# End of Subroutine MeridianMacro
 790
 791sub ParallelArcs {
 792	# This subroutine expects first and last parallels to do, and start and end meridians
 793	# for the circular arc. If there is a 5th argument, it is a string to append to
 794	# the comment line at the start of each parallel (e.g. ", Octant 1 West").
 795	# It is assumed that centre of circle is point A, with global coordinates $xA, $yA.
 796	# I will be lazy and use global variables instead of trying to pass hashes as 
 797	# parameters. Will use hashes %xP and %yP, each of which is indexed as "m,p".
 798	# It is assumed that every parallel multiple of 5 will be written to file MAJOR;
 799	# other meridians are written to file MINOR.
 800	my ($m, $p, $P, $FileHandle, $pStart, $pEnd, $mStart, $mEnd, $Memo);
 801	($pStart, $pEnd, $mStart, $mEnd, $Memo) = @_;
 802	if (not defined($Memo) ) {$Memo = "";}
 803
 804	#	Make sure start is less than end, to be able to use foreach
 805	if ($pStart > $pEnd) {
 806		$pStart = $pEnd;	$pEnd = $_[0];
 807	}
 808	foreach $p ($pStart .. $pEnd) {
 809		$P = ',' . $p ;
 810		if ( ($p % 5) == 0 ) {		# Major parallel
 811			$FileHandle = MAJOR;
 812		} else {				# Minor parallel
 813			$FileHandle = MINOR;
 814		}
 815		print $FileHandle ('{ Parallel ', $p, $Memo,' }', "\r");
 816		print $FileHandle ('AR ', Coordinate ($xA, $yA), "\\\r");
 817		print $FileHandle (Coordinate ($xP {$mStart . $P}, $yP {$mStart . $P}) );
 818		print $FileHandle (Coordinate ($xP {$mEnd . $P}, $yP {$mEnd . $P} ) );
 819		print $FileHandle (' [`;`]', "\r");
 820	}
 821}	# End of subroutine ParallelArcs
 822
 823
 824
 825sub ParallelMacro {
 826	# This subroutine expects first and last parallels to do, and start and end meridians
 827	# for the segmented parallel. If there is a 5th argument, it is a string to append to
 828	# the comment line at the start of each parallel (e.g. ", Octant 1 West").
 829	# I will be lazy and use global variables instead of trying to pass hashes as 
 830	# parameters. Will use hashes %xP and %yP, each of which is indexed as "m,p".
 831	# It is assumed that every parallel multiple of 5 will be written to file MAJOR;
 832	# other parallels are written to file MINOR.
 833	my ($m, $p, $P, $FileHandle, $pStart, $pEnd, $mStart, $mEnd, $Memo, $mS, $mE);
 834	($pStart, $pEnd, $mStart, $mEnd, $Memo) = @_;
 835	if (not defined($Memo) ) {$Memo = "";}
 836
 837	#	Make sure starts are less than ends, to be able to use foreach
 838	if ($pStart > $pEnd) {
 839		$pStart = $pEnd;	$pEnd = $_[0];
 840	}
 841	if ($mStart > $mEnd) {
 842		$mStart = $mEnd;	$mEnd = $_[2];
 843	}
 844	$mS = $mStart + 1;	$mE = $mEnd - 1;
 845
 846	# For each parallel to be done
 847	foreach $p ($pStart .. $pEnd) {
 848		$P = ',' . $p ;
 849		if ( ($p % 5) == 0 ) {		# Major parallel
 850			$FileHandle = MAJOR;
 851		} else {				# Minor parallel
 852			$FileHandle = MINOR;
 853		}
 854		print $FileHandle ('{ Parallel ', $p, $Memo, ' }', "\r");
 855		# Start point of line
 856		print $FileHandle ('LI ', Coordinate ($xP {$mStart . $P}, $yP {$mStart . $P} ) );
 857		foreach $m ($mS .. $mE) {	# For each meridian
 858			print $FileHandle ("\\\r", Coordinate ($xP {$m . $P}, $yP {$m . $P} ) );
 859			if (($m % 15) == 0) {
 860				# Every 15 segments or so, terminate and restart line command
 861				print $FileHandle ('[`;`]', "\rLI ", 
 862					Coordinate ($xP {$m . $P}, $yP {$m . $P} ) );
 863			}
 864		} 	# End of foreach meridian
 865		# Last point of parallel
 866		print $FileHandle (Coordinate ($xP{$mEnd.$P}, $yP{$mEnd.$P}), '[`;`]', "\r");
 867	}	#End of foreach parallel
 868}	# End of subroutine ParallelMacro
 869
 870sub ParallelMacro5 {
 871	# This subroutine expects first and last parallels to do, and start and end meridians
 872	# for the segmented parallel. It draws line segments for parallels by connecting nodes
 873	# of major meridians, that is, every 5th meridian, skipping in between nodes.
 874	# If there is a 5th argument, it is a string to append to the comment line at the start 
 875	# of each parallel (e.g. ", Octant 1 West").
 876	# I will be lazy and use global variables instead of trying to pass hashes as 
 877	# parameters. Will use hashes %xP and %yP, each of which is indexed as "m,p".
 878	# It is assumed that every parallel multiple of 5 will be written to file MAJOR;
 879	# other meridians are written to file MINOR.
 880	# This subroutine writes one parallel per LINE command (5� times 9 segments=45�).
 881
 882	my ($m, $p, $P, $FileHandle, $pStart, $pEnd, $mStart, $mEnd, $Memo, $mS, $mE);
 883	($pStart, $pEnd, $mStart, $mEnd, $Memo) = @_;
 884	if (not defined($Memo) ) {$Memo = "";}
 885
 886	#	Make sure starts are less than ends, to be able to use foreach
 887	if ($pStart > $pEnd) {
 888		$pStart = $pEnd;	$pEnd = $_[0];
 889	}
 890	if ($mStart > $mEnd) {
 891		$mStart = $mEnd;	$mEnd = $_[2];
 892	}
 893	$mS = $mStart + 1;	$mE = $mEnd - 1;
 894
 895	# For each parallel to be done
 896	foreach $p ($pStart .. $pEnd) {
 897		$P = ',' . $p ;
 898		if ( ($p % 5) == 0 ) {		# Major parallel
 899			$FileHandle = MAJOR;
 900		} else {				# Minor parallel
 901			$FileHandle = MINOR;
 902		}
 903		print $FileHandle ('{ Parallel ', $p, $Memo, ' }', "\r");
 904		# Start point of line
 905		print $FileHandle ('LI ', Coordinate ($xP {$mStart . $P}, $yP {$mStart . $P} ) );
 906		foreach $m ($mS .. $mE) {	# For each meridian
 907			if (($m % 5) == 0) {	# only do for major meridians
 908				print $FileHandle ("\\\r", Coordinate($xP{$m.$P},$yP{$m.$P} ) );
 909			}
 910		} 	# End of foreach meridian
 911		# Last point of parallel
 912		print $FileHandle (Coordinate ($xP{$mEnd.$P}, $yP{$mEnd.$P}), '[`;`]', "\r");
 913	}	#End of foreach parallel
 914}	# End of subroutine ParallelMacro5
 915
 916sub MLinesMacro {
 917	# Macro to draw concatenated lines, with coordinates from different octants.
 918	# First argument is file handle. The second argument must be a number: 
 919	# If it is 0, the following arguments are sets of x,y, and octant number.
 920	# If it is a number between 1 and 8, that is the number of the octant for all
 921	# subsequent pairs of x,y coordinates.
 922	# All coordinates are converted to M-map system.
 923	my ($FileHandle,$Octant, @a) = @_ ;
 924	my ($nArg, $i,$x,$y,$x2,$y2,$n);
 925	$nArg = @a - 1;
 926	$n = 0;	# Flag for starting point
 927	if ($Octant == 0) {
 928		for ($i = 0; $i < $nArg; $i=$i+3) {
 929			($x2,$y2) = MJtoG ($a[$i], $a[$i+1], $a[$i+2]);
 930			if ($n == 0) {	# Starting point; update flag
 931				print $FileHandle ('LI ', Coordinate ($x2,$y2) );
 932				$n ++ ;
 933			} else {	
 934				# Do line continuation, and this point
 935				print $FileHandle ("\\\r",Coordinate ($x2,$y2) );
 936			}
 937		}
 938	} else {
 939		for ($i = 0; $i < $nArg; $i=$i+2) {
 940			($x2,$y2) = MJtoG ($a[$i], $a[$i+1], $Octant);
 941			if ($n == 0) {	# Starting point; and update flag
 942				print $FileHandle ('LI ', Coordinate ($x2,$y2) );
 943				$n ++ ;
 944			} else {	
 945				# Do line continuation, and this point
 946				print $FileHandle ("\\\r",Coordinate ($x2,$y2) );
 947			}
 948		}
 949	}
 950	print $FileHandle (' [`;`]', "\r");
 951}	# End of subroutine MLinesMacro
 952
 953sub OctantMacro {
 954	# Macro to draw concatenated lines. Being lazy and using global arrays
 955	# @xOct and @yOct, and writing to file handle MACRO.
 956	my ($FileHandle,@x, @y, $n, $i);
 957	$FileHandle = MACRO;
 958	@x = @xOct;	@y = @yOct;
 959	$n = @xOct - 1;
 960	print $FileHandle ('LI ', Coordinate ($x[0],$y[0]) );
 961	foreach $i (1 .. $n) {
 962		# Do line continuation, and this point
 963		print $FileHandle ("\\\r",Coordinate ($x[$i],$y[$i]) );
 964	}
 965	print $FileHandle (' [`;`]', "\r");
 966}	# End of subroutine LinesMacro
 967
 968sub GridMacro {
 969	# Arguments are: filehandle, left, right, bottom, top, spacing, spacing of next
 970	# larger grid or 0 to not skip major lines.
 971	# Notice: every $skip line is not plotted; it will be part of the grid of next larger size.
 972	my ($filehandle, $left, $right, $bottom, $top, $delta, $skip) = @_ ;
 973	my ($i, $a, $b, $c);
 974	$a = 'LI [a,';
 975	$b = ',0,` `] [a,';
 976	$c = ',0,` `] [`;`]' . "\r";
 977	$i = $bottom;
 978	# If bottom is not a multiple of delta, start plotting at the first multiple of delta
 979	if ($i % $delta !=0) { $i = $i + $delta - ($i % $delta); }
 980	print $filehandle '{ Horizontal lines }', "\r" ;
 981	while ($i <= $top) {	# horizontal lines
 982		if ($skip != 0 && $i % $skip == 0) {
 983			$i = $i + $delta;
 984			 next;
 985		}
 986		print $filehandle $a,$left,',',$i,$b,$right,',',$i,$c;
 987		$i = $i + $delta;
 988	}
 989	$i = $left;
 990	# If left is not a multiple of delta, start plotting at the first multiple of delta
 991	if ($i % $delta !=0) { $i = $i + $delta - ($i % $delta); }
 992	print $filehandle '{ Vertical lines }', "\r" ;
 993	while ($i <= $right) {	# vertical lines
 994		if ($skip != 0 && $i % $skip == 0) {
 995			$i = $i + $delta;
 996			 next;
 997		}
 998		print $filehandle $a,$i,',',$top,$b,$i,',',$bottom,$c;
 999		$i = $i + $delta;
1000	}
1001}	# End of subroutine GridMacro
1002
1003sub LandMacro {
1004	# Will be lazy and read global arrays @Xs and @Ys, and write to filehandle LAND.
1005	my ($nPoints, $FileHandle, $n);
1006	$nPoints = @Xs - 2;
1007	$FileHandle = LAND;
1008	print $FileHandle ('LI ', Coordinate ($Xs [0], $Ys [0]) );
1009	foreach $n (1 .. $nPoints) {
1010		print $FileHandle ("\\\r", Coordinate ($Xs[$n], $Ys [$n] ) );
1011		if (($n % 15) == 0) {
1012			# Every 15 segments or so, terminate and restart line command
1013			print $FileHandle ('[`;`]', "\rLI ", Coordinate ($Xs[$n], $Ys [$n] ) );
1014		}
1015	}
1016	print $FileHandle (Coordinate ($Xs[$nPoints+1], $Ys[$nPoints+1]), '[`;`]', "\r");
1017}		# End of subroutine LandMacro
1018
1019# -  -  -  -  -  -  -  -  -    SUBROUTINES For Coordinate conversion   -  -  -  -  -  -  -  -  -  -  -  -
1020
1021sub MJtoG {
1022# Subroutine to convert (that is, do coordinate transformation of) x and y coordinates
1023# from Mary Jo's half-octant on its side to Gene's leaning, single octant coordinates, or
1024# to Gene's M-map (eight-octants) coordinates.
1025# 
1026# Subroutine returns converted x and y coordinates.
1027#
1028# Arguments are: 
1029# - x and y coordinates of point to convert.
1030# - Third argument, if given, is the Octant to convert to:
1031#	- no argument or 0 &#x2013; Gene's single-octant system, with y-axis on its left;
1032#	- 1, 2, 3 or 4 &#x2013; convert to M-map coordinates, respectively to first, second, third or
1033#		 fourth northern octant, from the left;
1034#	- 5, 6, 7 or 8 &#x2013; convert to M-map coordinates, respectively to fourth, first, second or
1035#		 third southern octant from the left.
1036#
1037# - In MJ's coordinates, point M is the origin, at (0,0), points M, A and G are on the positive
1038# x-axis, and point G is at (10000, 0).
1039# - In G's system, point M, L, J and P are on the positive y-axis, and point G is on the positive
1040# x-axis; in this system, point G is at coordinates (5000, 0).
1041# - From MJ's system to G's, there is a +60� rotation, and also a translation.
1042# - The M-map coordinate system is like G's system, except that the y-axis is 10000mm
1043# to the right, that is, the x-coordinates for the start octant are 10000mm smaller.
1044#
1045# I got the equations for rotation and translation from my pocketbook "The Universal
1046# Encyclopedia of Mathematics, with a Foreword by James R. Newman", �1964 by
1047# George Allen and Unwin, Ltd.; translated from original German language edition,
1048# pages 152, 153.
1049#
1050# This subroutine, in order to not have to recompute them each time the routine is called,
1051# uses a few global variables from the main program:
1052# $cos60 = cos (deg2rad(60));
1053# $sin60 = sin (deg2rad(60));
1054# $yTranslate = 10000 * $sin60;
1055
1056	my ($nArgs, $xnew, $ynew);
1057	my ($x, $y, $Octant) = @_ ;
1058	if (not defined ($Octant) ) { $Octant = 0; }
1059	if ($Octant == 0) {
1060		($xnew, $ynew) = Rotate ($x, $y, 60);
1061	} elsif ($Octant == 1) {
1062		($xnew, $ynew) = Rotate ($x, $y, 120);
1063		$xnew = $xnew - 10000;
1064	} elsif ($Octant == 2) {
1065		($xnew, $ynew) = Rotate ($x, $y, 60);
1066		$xnew = $xnew - 10000;
1067	} elsif ($Octant == 3) {
1068		($xnew, $ynew) = Rotate ($x, $y, 120);
1069		$xnew = $xnew + 10000;
1070	} elsif ($Octant == 4) {
1071		($xnew, $ynew) = Rotate ($x, $y, 60);
1072		$xnew = $xnew + 10000;
1073	} elsif ($Octant == 5) {
1074		($xnew, $ynew) = Rotate ((20000-$x), $y, 60);
1075		$xnew = $xnew + 10000;
1076	} elsif ($Octant == 6) {
1077		($xnew, $ynew) = Rotate ((20000-$x), $y, 120);
1078		$xnew = $xnew - 10000;
1079	} elsif ($Octant == 7) {
1080		($xnew, $ynew) = Rotate ((20000-$x), $y, 60);
1081		$xnew = $xnew - 10000;
1082	} elsif ($Octant == 8) {
1083		($xnew, $ynew) = Rotate ((20000-$x), $y, 120);
1084		$xnew = $xnew + 10000;
1085	} else {
1086		print "Error converting to M-map coordinates; there is no $Octant octant!\n";
1087		return ($x,$y);
1088	}
1089	$ynew = $ynew + $yTranslate;
1090	return ($xnew, $ynew);
1091
1092}	# End of subroutine MJtoM, which converts coordinates to octants on M-map
1093
1094sub Rotate {
1095	# Receives 3 arguments: x, y, and angle by which to rotate the coordinate system;
1096	# Expects that the axes will be rotated either 60� or 120�.
1097	# Returns new x and y values.
1098	# Uses $cos60, $sin60 global variables from main program, to minimize computations.
1099	my ($x, $y, $angle) = @_ ;
1100	my ($xnew, $ynew);
1101	if ( $angle == 60 ) {
1102		$xnew = $x * $cos60 + $y * $sin60;
1103		$ynew = -$x * $sin60 + $y * $cos60;
1104	} elsif ( $angle == 120 ) {
1105		$xnew = -$x * $cos60 + $y * $sin60;
1106		$ynew = -$x * $sin60 - $y * $cos60;
1107	}
1108	return $xnew, $ynew;
1109}	# End of subroutine Rotate
1110
1111sub Real2Map {
1112	# Arguments are real world longitude and latitude for one point, in decimal degrees.
1113	# West longitudes and south latitudes have negative values.
1114	# Returns x and y values in M-Map coordinates.
1115	#
1116	# Subroutine uses the following global variables, without modifying them:
1117	# $DeltaMEq, $EF, $AG, $xA, $yA, $xE, $yE, $xF, $yF, $xG, $xM, $yM.
1118	# Also, it calls subroutines LineIntersection, Length, Interpolate and MJtoG.
1119
1120	# @South are southern octants corresponding to northern octants 1, 2, 3 and 4; the 0
1121	# is just a place holder to facilitate correspondence.
1122	my (@South) = (0,6,7,8,5);
1123	# $m and $p are the meridian and parallel numbers in scaffold half-octant setting;
1124	# $Octant is the M-map octant of the real point; $Sign is for east or west side of
1125	# scaffold octant; $x and $y are coordinates in scaffold half-octant setting; $xNew
1126	# and $yNew are the coordinates converted to M-Map system, which are returned.
1127	my ($Octant, $m, $p, $Sign, $x, $y, $xNew, $yNew);
1128	# Local arrays @xJ, @yJ and @L are equivalent to the global hashes %xJ, %yJ, %L,
1129	# but just for this particular meridian. The other variables are temporary.
1130	my (@xJ, @yJ, @L, $distance, $extra, $flag);
1131	# Input values
1132	my ($Long, $Lat) = @_ ;
1133	# Determine the correct octant; Octant 1 is +160� to -110�; octant 4 is 70� to 160�
1134	$Octant = int ( ($Long + 200) / 90 ) + 1;
1135	# Make longitude fit within scaffold/template triangle, and determine if y value should
1136	# be positive or negative.
1137	$m = ( ($Long +200) - (90*($Octant - 1))) - 45;
1138	if ($m < 0) {
1139		$Sign = -1;	$m = -$m;
1140	} else {
1141		$Sign = 1;
1142	}
1143	# Fix the octant number, if necessary
1144	if ($Octant == 5) { $Octant = 1; }
1145	if ($Lat < 0) {
1146		$Octant = $South [$Octant];
1147		$p = -$Lat;
1148	} else {
1149		$p = $Lat;
1150	}
1151#	print "$Long, $Lat convert to $m, $p, $Sign, $Octant.\n";
1152	# Now need to calculate x and y coordinates for $m,$p in scaffold triangle.
1153	if ($m == 0) {	# Centre line meridian is straight, and easy to calculate
1154		$y = 0;
1155		$x = $xG - $p * $AG / 90;
1156	} else {	# Complicated jointed meridian
1157		# Calculate coordinates of joints
1158		$xJ [0] = $xA;	$yJ [0] = $yA;	# Polar end of meridian
1159		$distance = $m * $DeltaMEq;
1160		# Equatorial end of meridian
1161		if ($distance <= $yF) {
1162			$xJ [3] = $xG;
1163			$yJ [3] = $distance;
1164		} else {	# Past joint at point F
1165			$extra = $distance - $yF;
1166			$xJ [3] = Interpolate ($extra, $EF, $xF, $xE);
1167			$yJ [3] = Interpolate ($extra, $EF, $yF, $yE);
1168		}
1169		# Calculate joints by line intersection from A and M
1170		($flag, $xJ [1], $yJ [1]) = LineIntersection ($xM, $yM, 2*$m/3, $xA, $yA, $m);
1171		if ($flag ne "Intersect") { die "$flag while calculating 1st joint at $m.\n"; }
1172		($flag, $xJ [2], $yJ [2]) = 
1173			LineIntersection ($xM, $yM, 2*$m/3, $xJ [3], $yJ [3], $m/3);
1174		if ($flag ne "Intersect") { die "$flag while calculating 2nd joint at $m.\n"; }
1175		# Calculate length of meridian segments.
1176		$L [0] = Length ($xA, $yA, $xJ [1], $yJ [1]);
1177		$L [1] = Length ($xJ [1], $yJ [1], $xJ [2], $yJ [2]);
1178		$L [2] = Length ($xJ [2], $yJ [2], $xJ [3], $yJ [3]);
1179		$L [3] = $L [1] + $L [2];		# Torrid + Temperate
1180		$L [4] = $L [3] + $L [0];		# Total meridian length
1181		# Calculate parallel intersection &#x2013; location of this latitude on this meridian
1182		$distance = $p *  $L [4] / 90;	# total length corresponds to 90�
1183		if ($distance <= $L [2] ) {		# Torrid segment
1184			$x = Interpolate ($distance, $L [2], $xJ [3], $xJ [2]);
1185			$y = Interpolate ($distance, $L [2], $yJ [3], $yJ [2]);
1186		} elsif ($distance <= $L [3] ) {	# Temperate seg., past one joint
1187			$extra = $distance - $L [2];
1188			$x = Interpolate ($extra, $L [1], $xJ [2], $xJ [1]);
1189			$y = Interpolate ($extra, $L [1], $yJ [2], $yJ [1]);
1190		} else {				# Frigid segment, past two joints
1191			$extra = $distance - $L [3];
1192			$x = Interpolate ($extra, $L [0], $xJ [1], $xJ [0]);
1193			$y = Interpolate ($extra, $L [0], $yJ [1], $yJ [0]);
1194		}
1195	}	# Now have $x and $y coordinates in scaffold half-octant coordinates
1196	# Convert to M-map coordinates
1197	($xNew, $yNew) = MJtoG ($x,$Sign*$y, $Octant);
1198#	print "and x,y: $x, $y, converted to $xNew, $yNew.\n";
1199	return $xNew, $yNew;
1200}	# End of subroutine Real2Map
1201
1202# -  -  -  -  -  -  -  -  -    SUBROUTINES         -  -  -  -  -  -  -  -  -  -  -  -  -