```#! /usr/bin/perl -w
# 2011-06-30 This is program "MegamapMaker-prep" by Mary Jo Graça
```
# Go back to Technical Annexes, Cahill-Keyes Multi-Scale Megamap
```use Math::Trig;
# Calculate values that should be global, to minimize repeated calculations
(\$sin60, \$cos60, \$yTranslate, @Prelims) = Preliminary();

# BLOCK (1) : read file of longitude,latitude coordinates and convert to x,y coordinates
\$Skip[1] = "YES";
# BLOCK (2): prepare file of graticule data for Octant 3: "OctantTemplate.txt"
\$Skip[2] = "YES";
# BLOCK (3): Output values calculated in sub Preliminary to standard output (i.e. screen)
\$Skip[3] = "YES";
# BLOCK (4): Print coordinates of joints of whole-numbered meridians, to standard output.
# This is for template half-octant in MJ's coordinate system.
\$Skip[4] = "YES";
# BLOCK (5): Print coordinates of whole-numbered meridian-parallel points, to standard
# output. This is for template half-octant in MJ's coordinate system.
\$Skip[5] = "YES";

print "Type:\n- 1 to convert MAPGEN file (longitude, latitude) to Cahill-Keyes x,y coordinates;";
print "\n- 2 to prepare file \"OctantTemplate.txt\" needed to plot graticule;";
print "\n- 3 to output values calculated in sub Preliminary to standard output (i.e. screen);";
print "\n- 4 to output MJ template coordinates of joints of whole-numbered meridians to standard output;";
print "\n- 5 to output MJ template coordinates of whole-numbered meridian-parallel points to standard output;";
print "\n- or anything else to do nothing.\n";

# BLOCK (1) : read file of longitude,latitude coordinates and convert to x,y coordinates

unless (\$Skip[1]) {  # Option to skip making x,y file of Coastal Data
print "For which octant do you want to prepare data (1 to 8)?\n";
\$WantedOctant = <STDIN>;
chomp (\$WantedOctant);
if (!(\$WantedOctant =~ /[1-8]/ && length(\$WantedOctant) == 1) ) {
die "Octant must be 1, 2, 3, 4, 5, 6, 7 or 8!\n";
}
print "Type in name of MAPGEN file with coastal data\n";
\$FileIn = <STDIN>;
chomp (\$FileIn);
print "Type in name for output file to receive x,y coordinates.\n";
\$FileOut = <STDIN>;
chomp (\$FileOut);
# *** NOTE: Adding the path for Gene's directory with coastal data.
# *** THIS MUST BE CHANGED OR COMMENTED OUT TO RUN ON ANOTHER COMPUTER
#   \$FileIn = "/media/D:/My Documents/CKOG/COASTLINE-DATA-2/" . \$FileIn;
#   \$FileOut = "/media/D:/My Documents/CKOG/COASTLINE-DATA-2/" . \$FileOut;
# Output will be a plain text file with the following for each segment:
# line with "L,0" (without quotes), followed by lines of x,y (comma included) values.
#
# Read some Coastal Data, convert to M-map coordinates, and write to output file.
# File read is a zipped or unzipped file of MAPGEN data format: two columns ASCII flat
# file with: longitude tab latitude; at the start of each segment there is a line containing
# only "# -b".
# Variable \$Map can be set to "M" to output coordinates in M-map system, or to
# any other value, to output coordinates in Gene's one-octant system.
\$Map = "M";
# Set up a few variables
\$maxX=-99999;  \$maxY=-99999;  \$minX=999999;  \$minY=999999;
\$oldLong = 99999;  \$oldLat = 99999;
# Start output file
\$FileOut2 = "> " . \$FileOut;
# Open output file to receive x,y coordinates
open (FILEXY, \$FileOut2) || die "Sorry, I couldn't open output file, \"\$FileOut\".\n";
# Open coastal data file, taking into account whether or not it is compressed as .zip
if ( \$FileIn =~ /.\.zip\$/i || \$FileIn =~ /.\.gz\$/i) {
# If file ends in a letter plus a period plus "zip" or "ZIP" or gz
\$FileIn2 = "zcat \"" . \$FileIn . "\" | ";  # Compressed file
} else {
\$FileIn2 = \$FileIn;  # normal, uncompressed file
}
open (DATA, \$FileIn2) || die "Sorry, I couldn't open input file, \"\$FileIn2\".\n";
\$nData = 0;
print "This might take awhile. I'll let you know every 1,000 lines that I read.\n\n";
while (<DATA>) {
\$Line = \$_ ;
chomp(\$Line);  # Remove new-line character from end of line of data read
if (\$Line =~ /\r\$/ ) {
chop(\$Line);  # removing a return character, if there is one
}
\$nData ++;
if (\$nData % 1000 == 0) { print "Read \$nData lines of land data so far.\n"; }
if (index(\$Line,'#') < 0) {  # Process line with longitude, latitude
(\$Long,\$Lat) = split(/\t/,\$Line);
if ( ! (\$Long == \$oldLong && \$Lat == \$oldLat) ) {
# If this point is a repeat of the previous point on this segment, this section is
# not run, and this point is neither converted nor included in the output file.
\$oldLong = \$Long;  \$oldLat = \$Lat;
# COORDINATE CONVERSION
# Convert real longitude, latitude to template half-octant meridian and parallel
(\$m, \$p, \$Sign, \$Octant) = LLtoMP (\$Long, \$Lat);
if (\$Octant != \$WantedOctant) {
if ((\$Lat >= 0 && \$WantedOctant < 5) || (\$Lat <= 0 && \$WantedOctant > 4) ) {
# Point is in the correct norther or southern hemisphere
@Pairs = (0,6,7,8,5,4,1,2,3);
@WestLong = (0,160,-110,-20,70,70,160,-110,-20);
@EastLong = (0,-110,-20,70,160,160,-110,-20,70);
if (\$p == 0 && \$Octant == \$Pairs[\$WantedOctant] ) {
\$Octant = \$WantedOctant; # On boundary with its southern or northern octant
} elsif (\$p == 90) { # Check the poles
\$Octant = \$WantedOctant;
} elsif (\$m == 45) { # Check eastern and western boundaries
if (\$Long == \$WestLong[\$WantedOctant]) {
\$Octant = \$WantedOctant;
\$Sign = -1;
} elsif (\$Long == \$EastLong[\$WantedOctant]) {
\$Octant = \$WantedOctant;
\$Sign = 1;
}
}  # End of if \$p == 0 or \$p == 45
}  # End of checking if it is on the correct hemisphere
# If none of the above fixed the octant, then this point is not on the wanted octant,
# not even right on its boundary.
if (\$Octant != \$WantedOctant) {
die "OOPS!: On line \$nData, Point (\$Long,\$Lat) would be point \$Sign*\$m,\$p " .
"in octant \$Octant, not the wanted octant \$WantedOctant. Quitting!\n";
}
} # End of correcting octant number
# Convert template meridian, parallel to template half-octant x, y coordinates
(\$x, \$y) = MPtoXY (\$m, \$p, @Prelims);
# Convert template x, y coordinates to M-map or G's x and y coordinates.
# Variable \$Map was set previously in this "Skip" block.
if (\$Map eq "M") { # M-map coordinates
(\$xNew, \$yNew) = MJtoG (\$x, \$Sign*\$y, \$Octant, \$sin60, \$cos60, \$yTranslate);
} else { # G's single octant coordinates
# Using zero for octant number, which the Sub takes as Gene's octant coordinates
(\$xNew, \$yNew) = MJtoG (\$x, \$Sign*\$y, 0, \$sin60, \$cos60, \$yTranslate);
}
# Keep track of maximum and minimum values
if (\$xNew < \$minX) {\$minX = \$xNew;}
if (\$xNew > \$maxX) {\$maxX = \$xNew;}
if (\$yNew < \$minY) {\$minY = \$yNew;}
if (\$yNew > \$maxY) {\$maxY = \$yNew;}
# Make arrays of points for this segment to be written to output file
push (@Xs,\$xNew);
push (@Ys,\$yNew);
}  # End of skipping if the last point read was the same as the previous one
} elsif (defined(@Xs) ) {  # End of segment; write only if data points were read
# Write this segment of coastline with arrays @Xs and @Ys
\$nPoints = @Xs - 1;
print FILEXY ("L,0\n" );
foreach \$i (0 .. \$nPoints) {
# Values are multiplied by 100 to convert to 100ths of mm, and y-value is made
# negative because OOo Draw uses y positive downwards.
printf FILEXY ("%.0f,%.0f\n", \$Xs[\$i] * 100, -\$Ys[\$i] * 100);
}

@Xs = ();  @Ys = ();  # Empty arrays, ready for next segment
\$oldLong = 99999;  \$oldLat = 99999;  # Reset previous values to impossible values
}
}
close (DATA);
if (\$Line ne "\# -b") { # Draw last segment, if it hasn't been drawn
# Write this segment of coastline with arrays @Xs and @Ys
\$nPoints = @Xs - 1;
print FILEXY ("L,0\n" );
foreach \$i (0 .. \$nPoints) {
# Values are multiplied by 100 to convert to 100ths of mm, and y-value is made
# negative because OOo Draw uses y positive downwards.
printf FILEXY ("%.0f,%.0f\n", \$Xs[\$i] * 100, -\$Ys[\$i] * 100);
}
}
close (FILEXY);
print "Read a total of \$nData lines, and wrote file \"\$FileOut\". \n";
print \$minX, ' <= x <= ',\$maxX," and ",\$minY, ' <= y <= ', \$maxY, "\n";

}  # End of skipping making x,y file of Coastal Data

# BLOCK (2): prepare file of graticule data for Octant 3: "OctantTemplate.txt"

unless (\$Skip[2]) { # Option to skip writing Octant 3 with graticule for OOo macro
(\$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE, \$xF, \$yF, \$xG, \$yG,
\$xM, \$yM, \$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R, \$DeltaMEq) = @Prelims;
# The output file will have: minor meridians, minor parallels, major meridians,
# major parallels, octant's boundary, in that order. This requires saving the code in
# variables \$Major and \$Minor until all the meridians and parallels have been calculated.

# Prepare meridians
# Meridians multiple of 5° are drawn from point A with color 2; other meridians are
# drawn from parallel 85° with color 3. Array's index is: 0 = polar start (point A or
# parallel 85°); 1 = frigid joint; 2 = torrid joint; 3 = equator.
# Meridian 0° has no joints; it goes from A to G.
\$Minor = "";
\$Major = "L,2\n";
(\$xNew, \$yNew) = MJtoG (\$xA, \$yA, 3, \$sin60, \$cos60, \$yTranslate);
\$Major .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
(\$xNew, \$yNew) = MJtoG (\$xG, \$yG, 3, \$sin60, \$cos60, \$yTranslate);
\$Major .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);

# Other than meridian 0°, do both positive and negative meridians (i.e. 5° and -5°).
# Meridian 45° not done because it would be covered by the octant's boundary.
foreach \$m (1 .. 44) {
(\$x[3], \$y[3], \$x[2], \$y[2], \$x[1], \$y[1], \$Lt, \$Lm) =
Joints (\$m, \$xA, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq);
if (\$m %5 == 0) {  # Every 5th meridian starts at point A
(\$x[0], \$y[0]) = (\$xA, \$yA);
\$Color = 2;
} else {  # Minor meridians start at 85° of latitude
(\$x[0], \$y[0]) = MPtoXY (\$m, 85, @Prelims);
\$Color = 3;
}
# Do the meridian in the positive half-octant
\$temp = "L,\$Color\n";
foreach \$i (0..3) {
(\$xNew, \$yNew) = MJtoG (\$x[\$i], \$y[\$i], 3, \$sin60, \$cos60, \$yTranslate);
\$temp .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
}
# Do the meridian in the negative half-octant
\$temp .= "L,\$Color\n";
foreach \$i (0..3) {
(\$xNew, \$yNew) = MJtoG (\$x[\$i], -\$y[\$i], 3, \$sin60, \$cos60, \$yTranslate);
\$temp .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
}
if (\$m %5 == 0) {  # Major meridian
\$Major .= \$temp;
} else {  # Minor meridian
\$Minor .= \$temp;
}
} # End of meridians
# Prepare parallels; parallel 0° not done because it is on the octant's boundary.
foreach \$p (1 .. 89) {  # Calculate all values
foreach \$m (0 .. 45) {
(\$x[\$m][\$p], \$y[\$m][\$p]) = MPtoXY (\$m, \$p, @Prelims);
}
}
foreach \$p (1 .. 89) {  # Convert to octant 3 and write out values
\$temp1 = ""; \$temp2 = "";
# Positive half-octant
foreach \$m (0..45) {
(\$xNew, \$yNew) = MJtoG (\$x[\$m][\$p], \$y[\$m][\$p], 3, \$sin60, \$cos60, \$yTranslate);
\$temp1 .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
}
# Negative half-octant
foreach \$m (0..45) {
(\$xNew, \$yNew) = MJtoG (\$x[\$m][\$p], -\$y[\$m][\$p], 3, \$sin60, \$cos60, \$yTranslate);
\$temp2 .= sprintf ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
}
if (\$p %5 == 0) { # Major parallel, drawn with color 2
\$Major .= "L,2\n" . \$temp1 . "L,2\n" . \$temp2
} else { # Minor parallel, drawn with color 3
\$Minor .= "L,3\n" . \$temp1 . "L,3\n" . \$temp2
}
}  # End of parallels
# Open file, write out meridians and parallels, minor first and then the major ones.
open (MACRO, ">OctantTemplate.txt");  # Name for output file with OOo macro
print MACRO \$Minor, \$Major;
# Octant's outline (drawn with color 1)
print MACRO "L,1\n";
@x = (\$xA, \$xB, \$xD, \$xE, \$xF, \$xG, \$xF, \$xE, \$xD, \$xB, \$xA);
@y = (\$yA, \$yB, \$yD, \$yE, \$yF, \$yG, -\$yF, -\$yE, -\$yD, -\$yB, \$yA);
foreach \$i (0..10) {
(\$xNew, \$yNew) = MJtoG (\$x[\$i], \$y[\$i], 3, \$sin60, \$cos60, \$yTranslate); # 3 for octant 3
printf MACRO ("%.0f,%.0f\n", 100*\$xNew,-100*\$yNew);
}
close (MACRO);
}  # End of skip writing Octant 3 with graticule for OOo Basic macro

# BLOCK (3): Output values calculated in sub Preliminary to standard output (i.e. screen)

unless (\$Skip[3]) {	# Option to skip printing output of sub Preliminary
(\$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE, \$xF, \$yF, \$xG, \$yG,
\$xM, \$yM, \$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R, \$DeltaMEq) = @Prelims;
print "\n\nPoints";
foreach \$i qw(A B C D E F G M T) {printf "\t%s", \$i; }
print "\nx";
foreach \$i (\$xA,\$xB,\$xC,\$xD,\$xE,\$xF,\$xG,\$xM,\$xT) {printf "\t%.4f", \$i};
print "\ny";
foreach \$i (\$yA,\$yB,\$yC,\$yD,\$yE,\$yF,\$yG,\$yM,\$yT) {printf "\t%.4f", \$i};
print "\n\nLengths";
foreach \$i qw(AG AB BD GF BDE GFE R DeltaMEq) {printf "\t%s", \$i; }
print "\n";
foreach \$i (\$AG,\$AB,\$BD,\$GF,\$BDE,\$GFE,\$R,\$DeltaMEq) {printf "\t%.4f", \$i};
print "\n\n";
}  # End of skip printing output of sub Preliminary

# BLOCK (4): Print coordinates of joints of whole-numbered meridians, to standard output.
# This is for template half-octant in MJ's coordinate system.

unless (\$Skip[4]) {	# Option to skip calculating and printing out whole-numbered meridians
# Meridians multiple of 5° are drawn from point A; other meridians are drawn from
# parallel 85°. Array's second index is: 0 = polar start (point A or parallel 85°);
# 1 = frigid joint; 2 = torrid joint; 3 = equator.
(\$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE, \$xF, \$yF, \$xG, \$yG,
\$xM, \$yM, \$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R, \$DeltaMEq) = @Prelims;
# Meridian 0° has no joints; will make elements 1 and 2 equal to equatorial point, point G.
(\$x[0][0], \$y[0][0]) = (\$xA, \$yA);
(\$x[0][1], \$y[0][1]) = (\$xG, \$yG);
(\$x[0][2], \$y[0][2]) = (\$xG, \$yG);
(\$x[0][3], \$y[0][3]) = (\$xG, \$yG);
foreach \$m (1 .. 45) {
# Meridians multiple of 5° are drawn from point A; other meridians are drawn from
# parallel 85°.
if (\$m %5 == 0) {  # Every 5th meridian starts at point A
(\$x[\$m][0], \$y[\$m][0]) = (\$xA, \$yA);
} else {  # Minor meridians start at 85° of latitude
(\$x[\$m][0], \$y[\$m][0]) = MPtoXY (\$m, 85, @Prelims);
}
(\$x[\$m][3], \$y[\$m][3], \$x[\$m][2], \$y[\$m][2], \$x[\$m][1], \$y[\$m][1], \$Lt, \$Lm) =
Joints (\$m, \$xA, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq);
}
print "\n\nJoints\nx";
foreach \$m (0 .. 45) { printf "\t%d", \$m; } print "\n";
foreach \$i (0 .. 3) {
print \$i;
foreach \$m (0 .. 45) {
if (defined (\$x[\$m][\$i])) {printf "\t%.4f", \$x[\$m][\$i];
} else { printf "\t undef"; }
}
print "\n";
}
print "\ny";
foreach \$m (0 .. 45) { printf "\t%d", \$m; } print "\n";
foreach \$i (0 .. 3) {
print \$i;
foreach \$m (0 .. 45) {
if (defined (\$y[\$m][\$i])) {printf "\t%.4f", \$y[\$m][\$i];
} else { printf "\t undef"; }
}
print "\n";
}
print "\n";
}  # End skipping calculation and printing out of whole-numbered meridians

# BLOCK (5): Print coordinates of whole-numbered meridian-parallel points, to standard
# output. This is for template half-octant in MJ's coordinate system.

unless (\$Skip[5]) { # Option to skip calculating all whole-numbered meridian-parallel points
foreach \$p (0 .. 90) {
foreach \$m (0 .. 45) {
(\$x[\$m][\$p], \$y[\$m][\$p]) = MPtoXY (\$m, \$p, @Prelims);
}
}
print "\n\nPoints\nx";
foreach \$m (0 .. 45) { printf "\t%d", \$m; } print "\n";
foreach \$p (0 .. 90) {
print \$p;
foreach \$m (0 .. 45) {
if (defined (\$x[\$m][\$p])) {printf "\t%.4f", \$x[\$m][\$p];
} else { printf "\t undef"; }
}
print "\n";
}
print "\ny";
foreach \$m (0 .. 45) { printf "\t%d", \$m; } print "\n";
foreach \$p (0 .. 90) {
print \$p;
foreach \$m (0 .. 45) {
if (defined (\$y[\$m][\$p])) {printf "\t%.4f", \$y[\$m][\$p];
} else { printf "\t undef"; }
}
print "\n";
}
}  # End skipping calculation of all whole-numbered meridian-parallel points

print "All done.\n";
} else {
print "Nothing done.\n";
}

# -  -  -  -  -  -  -  -  -  S U B R O U T I N E S  -  -  -  -  -  -  -  -
sub Preliminary{
# Calculates and returns an array with 29 values (0 to 28), in this order:
# (for use of subs MJtoG and Rotate) \$sin60, \$cos60, \$yTranslate,
# (x and y coordinates of points) \$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE,
# \$xF, \$yF, \$xG, \$yG, \$xM, \$yM, \$xT, \$yT, (lengths) \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE,
# \$R, \$DeltaMEq.

use Math::Trig;
# Values that will be returned
my (\$sin60, \$cos60, \$yTranslate);
my (\$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE, \$xF, \$yF, \$xG, \$yG, \$xM, \$yM);
my (\$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R, \$DeltaMEq);
# Variables temporary to this sub
my (\$xN, \$yN, \$MB, \$MN, \$xU, \$yU, \$k, \$xV, \$yV);

# Some constants for use by subs MJtoG and Rotate, which do coordinate axis
# transformation. Angle of rotation is 60°.  Point G is (10000,0) in MJ, and (5000,0) in G
\$yTranslate = 10000 * \$sin60;

# Given input
\$xM = 0;      \$yM = 0;      # Point M is the origin of the axes
\$xG = 10000;      \$yG = 0;      # Point G, at center of base of octant
\$xA = 940;      \$yA = 0;      # Point A at apex of octant
# Other points and lengths of interest, relating to scaffold triangle and half-octant
\$xN = \$xG;      \$yN = \$xG * tan (deg2rad(30));      # Point N, point of triangle MNG
(\$xB, \$yB) = LineIntersection (\$xM, \$yM, 30, \$xA, \$yA, 45);      # Point B
\$AG = \$xG - \$xA;
\$AB = Length (\$xA, \$yA, \$xB, \$yB);
\$MB = Length (\$xM, \$yM, \$xB, \$yB);
\$MN = Length (\$xM, \$yM, \$xN, \$yN);
# Calculate point D, considering that length DN = MB
\$xD = Interpolate (\$MB, \$MN, \$xN, \$xM);      # D is away from N as B is away from M
\$yD = Interpolate (\$MB, \$MN, \$yN, \$yM);
\$xF = \$xG;
\$yF = \$yN - \$MB;
# Distance from point E to point N = distance from point A to point M = xA; calculate E
\$xE = \$xN - \$xA * sin (deg2rad(30));
\$yE = \$yN - \$xA * cos (deg2rad(30));
\$GF = \$yF;
\$BD = Length (\$xB, \$yB, \$xD, \$yD);
\$BDE = \$BD + \$AB;      # Length AB = length DE
\$GFE = \$AB + \$GF;      # Length AB = length FE
\$DeltaMEq = \$GFE / 45;      # 45 meridian spacings along equator for half an octant
# Calculate Point T: First calculate point U = (30°,73°). Radius to circular arc of 73° =
# 15° x 104mm/° + 2° x 100 mm/° = 1760 mm.
\$xU = \$xA + 1760 * cos (deg2rad(30));
\$yU = 1760 * sin (deg2rad(30));
# Point T is at intersection of line BD with line from point U perpendicular to BD.
# Since line BD is 30° from horizontal, perpendicular line is -60° from horizontal.
(\$xT, \$yT) = LineIntersection ( \$xU, \$yU, -60, \$xB, \$yB, 30);

# To calculate point C, must first calculate point V = (29°, 15°).
# First calculate joints of meridian 29°
(\$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$Lt, \$Lm) =
Joints (29, \$xA, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq);
# Next need point on parallel 73° for this meridian 29°; really, only \$Lf is needed.
(\$xP73, \$yP73, \$Lf) = Parallel73 (29, \$xA, \$xT, \$yT, \$xJf, \$yJf);
# Do something with \$xP73 and \$yP73, only so that the compiler doesn't complain
# that they were used only once; I really don't need them now.
\$xP73 = 1 * \$xP73; \$yP73 = 1 * \$yP73;
# Parallels are equally spaced between the equator and latitude 73° in this zone.
# Both torrid and frigid joints are at latitudes lower than 73° in this region.
# To find point V, calculate length from equator to parallel 15°, along meridian.
# (\$Lt + \$Lm + \$Lf) = length from equator to parallel 73° on meridian 29°.
\$L = 15 * (\$Lt + \$Lm + \$Lf) / 73;
if (\$L <= \$Lt) {
# Measure length along the torrid segment, from the equator
\$xV = Interpolate (\$L, \$Lt, \$xJe, \$xJt);
\$yV = Interpolate (\$L, \$Lt, \$yJe, \$yJt);
} else {
# Measure length along the middle segment, from the torrid joint
\$L = \$L - \$Lt;
\$xV = Interpolate (\$L, \$Lm, \$xJt, \$xJf);
\$yV = Interpolate (\$L, \$Lm, \$yJt, \$yJf);
}
# Point C is the center of circular arc for parallel 15° with ends at points D and V, and,
# therefore, it is equidistant from both. Radius, R = CD = CV. Thus:
# \$R^2 = (\$xD - \$xC)^2 + (\$yD - \$yC)^2 = (\$xV - \$xC)^2 + (\$yV - \$yC)^2
# Point C is also on line MD, which has angle 30° with horizontal. M = (0 mm, 0 mm).
# Thus, \$yC / \$xC = tan(deg2rad(30)) = 1 / sqrt(3) ; letting \$k = sqrt(3), last equation is
# equivalent to \$xC = \$k * \$yC. Replacing this in the first equation and solving for \$yC,
# yields:
\$k = sqrt(3);
\$yC = (\$xV * \$xV + \$yV  * \$yV - \$xD * \$xD - \$yD * \$yD) /
(2 * (\$k * \$xV + \$yV  - \$k * \$xD - \$yD ) );
\$xC = \$k * \$yC;
\$R = Length (\$xC, \$yC, \$xD, \$yD);

# Return values needed by main program
return (\$sin60, \$cos60, \$yTranslate, \$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE,
\$xF, \$yF, \$xG, \$yG, \$xM, \$yM, \$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R,
\$DeltaMEq);
} # End of sub Preliminary

sub Equator {
# Sub calculates equatorial point for a meridian, and returns (\$xJe, \$yJe).
# Input is the wanted meridian, \$m, and the following values calculated in sub Preliminary:
# \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq
use Math::Trig;
my (\$m,  \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq) = @_;  # Input arguments
my (\$xJe, \$yJe);  # Values to be returned
my (\$L);  # Variable used just within this sub

# Calculate point Je, the Intersection of meridian with equator, as in zone (d)
\$L = \$DeltaMEq * \$m;
if (\$L <= \$GF) {
\$xJe = \$xG;
\$yJe = \$L
} else {
# Past point F; find point on line FE, a distance L from point G, along equator.
# Length FE = length AB
\$L = \$L - \$GF;    # Part of length on segment FE
\$xJe = Interpolate (\$L, \$AB, \$xF, \$xE);
\$yJe = Interpolate (\$L, \$AB, \$yF, \$yE);
}
return (\$xJe, \$yJe);
}  # End sub Equator

sub Joints {
# Sub calculates equatorial, torrid and frigid joints for given meridian, and lengths of
# middle segments. Returns are array: (\$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$Lt, \$Lm).
# \$xJe and \$yJe are calculated by calling sub Equator.
# Input is the wanted meridian, \$m, and the following values calculated in sub Preliminary:
# \$xA, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq
use Math::Trig;
my (\$m, \$xA, @Prelims) = @_;  # Input arguments
# Parse the input arguments
my (\$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq) = @Prelims ;

my (\$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$Lt, \$Lm);  # Values to be returned
my (\$L);  # Variable just within this sub

# Calculate point Je, the Intersection of meridian with equator
(\$xJe, \$yJe) = Equator (\$m, @Prelims);

# Calculate torrid joint, Jt, the intersection of line of angle (\$m/3) starting at point Je
# with line of angle (2/3 * \$m) starting at point M = (0 mm, 0 mm)
(\$xJt, \$yJt) = LineIntersection (0, 0, 2*\$m/3, \$xJe, \$yJe, \$m/3);

# Calculate frigid joing, Jf, the intersection of line of angle (\$m) starting at point A
# with line of angle (2/3 * \$m) starting at point M. Point A = (\$xA mm, 0 mm).
(\$xJf, \$yJf) = LineIntersection (\$xA, 0, \$m, 0, 0, 2*\$m/3);

# Calculate lengths of torrid segment, \$Lt = Je to Jt, and of middle segment, \$Lm = Jt to Jf
\$Lt = Length (\$xJe, \$yJe, \$xJt, \$yJt);
\$Lm = Length (\$xJt, \$yJt, \$xJf, \$yJf);

return (\$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$Lt, \$Lm);
} # End sub Joints

sub Parallel73 {
# Sub calculates parallel 73° for a meridian, and length from that point to the frigid joint.
# Note: if the point is on the middle segment, the length, Lf, to the frigid joint is given as
# a negative number; this only happens for some of the meridians between 44° and 45°.
# Returns are (\$xP73, \$yP73, \$Lf).
# Input is the wanted meridian, \$m; \$xA, \$xT, and \$yT, calculated in sub Preliminary;
# \$xJf, and \$yJf, the frigid joint, calculated in sub Joints.
use Math::Trig;
my (\$m, \$xA, \$xT, \$yT, \$xJf, \$yJf) = @_;  # Input arguments
my (\$xP73, \$yP73, \$Lf);  # Values to be returned
my (\$x, \$y);  # Values used only in this sub
# Calculate point P73 = (\$m, 73°) and length \$Lf = distance from Jf to P73 (negative if
# on middle segment).
if (\$m <= 30) {
# Circular arc portion:
# Radius to circular arc of 73° = 15° x 104mm/° + 2° x 100 mm/° = 1760 mm.
\$xP73 = \$xA + 1760 * cos (deg2rad(\$m));
\$yP73 = 1760 * sin (deg2rad(\$m));
# Calculate length \$Lf = distance from point Jf to point P73
\$Lf = Length (\$xJf, \$yJf, \$xP73, \$yP73);
} else {
# Straight portion of parallel 73°. Calculate point P73, at the intersection of line UT
# (angled -60° with the horizontal) with frigid segment of meridian \$m, which is
# angled +\$m ° and passes through point . Point U = (30°, 73°) was
# used to calculate point T, in sub Preliminary.
(\$xP73, \$yP73) = LineIntersection(\$xT, \$yT, -60, \$xJf, \$yJf, \$m);

# Calculate length \$Lf, from point Jf to point P73
\$Lf = Length (\$xJf, \$yJf, \$xP73, \$yP73);
if (\$m > 44) {
# Point P73 is on middle meridian segment for some of these meridians; check if it is.
# Calculate intersection of line UT with middle segment, angled +(2/3*\$m)°.
(\$x, \$y) = LineIntersection (\$xT, \$yT, -60, \$xJf, \$yJf, (2/3*\$m) );
if (\$x > \$xP73) {
# Correct intersection is on middle segment; correct point and length
\$xP73 = \$x;
\$yP73 = \$y;
\$Lf = - Length (\$xJf, \$yJf, \$xP73, \$yP73);  # Recalculating length and making it negative
}  # End of correction
}  # End of checking if it is on middle segment
}
return (\$xP73, \$yP73, \$Lf);
}  # End sub Parallel73

sub MPtoXY {
# Sub converts half-octant meridian,parallel to x,y coordinates.
# Arguments are meridian, parallel, and array output by sub Preliminary not including
# its first 3 values.
# Sub returns (x,y).
use Math::Trig;
my (\$m, \$p, @Prelims) = @_;  # Input arguments
# Parse the input arguments
my (\$xA, \$yA, \$xB, \$yB, \$xC, \$yC, \$xD, \$yD, \$xE, \$yE, \$xF, \$yF, \$xG, \$yG,
\$xM, \$yM, \$xT, \$yT, \$AG, \$AB, \$BD, \$GF, \$BDE, \$GFE, \$R, \$DeltaMEq) = @Prelims ;
my (\$x, \$y);  # Variables to be returned

# Extra variables used in this sub
my (\$L, \$xP73, \$yP73, \$xP75, \$yP75, \$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$f73, \$f75,
\$Lt, \$Lm, \$Lf, \$L73, \$xPm, \$yPm, \$flag);

if (\$m == 0) {  # Zones (a) and (b) on the center line of octant, on the base of
# scaffold half-triangle
\$y = 0;
if (\$p >= 75) {  # Zone (a) – frigid center line
# Parallel spacing is 104 mm/° from 75° to 90° of latitude, measured from point A (pole)
\$x = \$xA + (90 - \$p) * 104;

} else {  # Zone(b) – torrid/temperate center line
# Parallel spacing is 100 mm/° from 0° to 75° of latitude; measure from equator, that is
# from point G.
\$x = \$xG - \$p * 100;
}

} elsif (\$p >= 75) {  # Zone (c) – polar region with circular parallels spaced 104 mm/°
# Meridian \$m starts at point A and makes angle \$m (in degrees) with line AG
\$L = 104 * (90 - \$p);
\$x = \$xA + \$L * cos ( deg2rad(\$m) );
\$y = \$L * sin( deg2rad(\$m) );

} elsif (\$p == 0) {  # Zone (d) – equator
(\$x, \$y) = Equator (\$m, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq);

} elsif (\$p >= 73 && \$m <= 30) {  # Zone (e) – frigid region with circular parallels
# spaced 100 mm/° between 73° and 75° of latitude; meridian \$m starts at point A
# and makes angle \$m with line AG.
# Length from A to parallel 75° = 1560 mm = 104 mm/° x (90° - 75°).
\$L = 1560 + (75 - \$p) * 100;
\$x = \$xA + \$L * cos ( deg2rad(\$m) );
\$y = \$L * sin ( deg2rad(\$m) );

} elsif (\$m == 45) {  # Outer boundary of octant, zones (f), (g), and (h)
if (\$p <= 15) {  # Zone (f) – torrid zone of outer boundary, that is, along line ED
# E = 0° and D = 15° of latitude. Parallels are equally spaced within this zone.
\$x = Interpolate (\$p, 15, \$xE, \$xD);
\$y = Interpolate (\$p, 15, \$yE, \$yD);

} elsif (\$p <= 73) {  # Zone (g) – temperate zone of outer boundary, that is, along DT.
# D = 15°; T = 73°. Parallels are equally spaced within this zone. 73° - 15° = 58°
\$L = \$p - 15;
\$x = Interpolate (\$L, 58, \$xD, \$xT);
\$y = Interpolate (\$L, 58, \$yD, \$yT);

} else {  # Zone (h) – frigid supple zone of outer boundary
# Calculate point P75 = (45 °, 75°), point at parallel 75° on this meridian
# Length from A to parallel 75° = 1560 mm = 104 mm/° x (90° - 75°).
\$xP75 = \$xA + 1560 * cos (deg2rad(45));
\$yP75 = 1560 * sin (deg2rad(45));

# Calculate length \$Lf = parallel 73 (which is point T) to frigid joint (point B)
\$Lf = Length (\$xT, \$yT, \$xB, \$yB);
# Calculate length from \$Lf75 = distance from frigid joint (point B) to point P75
\$Lf75 = Length (\$xB, \$yB, \$xP75, \$yP75);
# Length from P75 to P73 covers 2°
\$L = (75 - \$p) * (\$Lf75 + \$Lf) / 2;  # Distance from parallel 75° to parallel p°
if (\$L <= \$Lf75) {
# Wanted latitude is on frigid segment, parallel 75° (P75) to B
\$x = Interpolate (\$L, \$Lf75, \$xP75, \$xB);
\$y = Interpolate (\$L, \$Lf75, \$yP75, \$yB);
} else {
# Wanted latitude is on segment B to T
\$L = \$L - \$Lf75;
\$x = Interpolate (\$L, \$Lf, \$xB, \$xP73);
\$y = Interpolate (\$L, \$Lf, \$yB, \$yP73);
}
}  # End of zones (f), (g), and (h); more specifically, end of zone (h)
} else {  # Zones (i), (j), (k), and (l) which require more complicated calculations
# Need to calculate meridian joints and segment lengths for this meridian.
(\$xJe, \$yJe, \$xJt, \$yJt, \$xJf, \$yJf, \$Lt, \$Lm) =
Joints (\$m, \$xA, \$xE, \$yE, \$xF, \$yF, \$xG, \$AB, \$GF, \$DeltaMEq);

# Calculate point P73 = ( \$m, 73°), point at latitude 73° on this meridian and distance
# from that point to frigid joint, \$Lf. These may later be modified for zones (j), (k) and (l).
(\$xP73, \$yP73, \$Lf) = Parallel73 (\$m, \$xA, \$xT, \$yT, \$xJf, \$yJf);

if (\$m <= 29) {  # Zone (i) – torrid and temperate areas of central two-thirds of octant
# Parallels are equally spaced between the equator and latitude 73° in this zone.
# Both torrid and frigid joints are at latitudes lower than 73° in this region.
# Calculate length from equator to point (\$m, \$p), along this meridian \$m.
\$L = \$p * (\$Lt + \$Lm + \$Lf) / 73;
if (\$L <= \$Lt) {
# Measure length along the torrid segment, from the equator
\$x = Interpolate (\$L, \$Lt, \$xJe, \$xJt);
\$y = Interpolate (\$L, \$Lt, \$yJe, \$yJt);
} elsif (\$L <= (\$Lt + \$Lm)) {
# Measure length along the middle segment, from the torrid joint
\$L = \$L - \$Lt;
\$x = Interpolate (\$L, \$Lm, \$xJt, \$xJf);
\$y = Interpolate (\$L, \$Lm, \$yJt, \$yJf);
} else {
# Measure length along the frigid segment
\$L = \$L - \$Lt - \$Lm;
\$x = Interpolate (\$L, \$Lf, \$xJf, \$xP73);
\$y = Interpolate (\$L, \$Lf, \$yJf, \$yP73);
}    # end of area (i)

} else {  # Supple zones (j), (k), and (l): 29° < \$m < 45° and 0° < \$p < 73

if (\$p >= 73) {  # Zone (j) – frigid supple zone
# Calculate point P75 = (\$m °, 75°), point at parallel 75° on this meridian
# Length from A to parallel 75° = 1560 mm = 104 mm/° x (90° - 75°).

\$xP75 = \$xA + 1560 * cos (deg2rad(\$m));
\$yP75 = 1560 * sin (deg2rad(\$m));
# Calculate length from \$Lf75 = distance from frigid joint, Jf, to point P75
\$Lf75 = Length (\$xJf, \$yJf, \$xP75, \$yP75);
# Length from P75 to P73 covers 2°; remember that Lf, from P73 to Jf is sometimes
# negative for a few meridians between 44° and 45°.
\$L = (75 - \$p) * (\$Lf75 - \$Lf) / 2;  # Distance from parallel 75° to parallel p°
if (\$L <= \$Lf75) {
# Wanted latitude is on frigid segment
\$x = Interpolate (\$L, \$Lf75, \$xP75, \$xJf);
\$y = Interpolate (\$L, \$Lf75, \$yP75, \$yJf);
} else {
# Wanted latitude is on middle segment
\$L = \$L - \$Lf75;
\$x = Interpolate (\$L, -\$Lf, \$xJf, \$xP73);
\$y = Interpolate (\$L, -\$Lf, \$yJf, \$yP73);
}

} else {  # Zones (k) and (l)
# Calculate point P15 = (m, 15°), that is, point on this meridian at latitude 15°, which
# is at intersection of meridian with circular arc of center C and radius R. Also
# calculate length L15 = distance from equator (Je) to P15.
# Try middle segment first, since most, if not all, parallel 15° points are in this segment
(\$flag, \$xP15, \$yP15) = CircleLineIntersection (\$xC, \$yC, \$R, \$xJt, \$yJt, \$xJf, \$yJf);
if (\$flag == 1) {    # Found the intersection point in middle segment
\$L15 = \$Lt + Length (\$xJt, \$yJt, \$xP15, \$yP15);
} else {    # Intersection point is in torrid segment
(\$flag, \$xP15, \$yP15) = CircleLineIntersection (\$xC, \$yC, \$R, \$xJe, \$yJe, \$xJt, \$yJt);
if (\$flag==0) {    # Hmmm... no intersection!
die " No line-circular arc intersection for M \$m, at parallel 15! Terminating.\n";
}
\$L15 = \$Lt - Length (\$xJt, \$yJt, \$xP15, \$yP15);
}

if (\$p <= 15) {  # Zone (k) – torrid supple zone
# Parallels equally spaced between equator and 15°
\$L = \$p * \$L15 / 15;
if (\$L <= \$Lt) {  # Point is in torrid segment
\$x = Interpolate (\$L, \$Lt, \$xJe, \$xJt);
\$y = Interpolate (\$L, \$Lt, \$yJe, \$yJt);
} else {  # Point is in middle segment
\$L = \$L - \$Lt;
\$x = Interpolate (\$L, \$Lm, \$xJt, \$xJf);
\$y = Interpolate (\$L, \$Lm, \$yJt, \$yJf);
}

} else {  # Zone (l) – middle supple zone
# Parallels equally spaced between 15° and 73°.
# Will measure from the equator. (\$Lt+\$Lm+\$Lf) = equator to P73. 58° = 73° - 15°
\$L = \$L15 + (\$p - 15) * ((\$Lt + \$Lm + \$Lf) - \$L15) / 58;
if (\$L <= \$Lt) {  # On torrid segment
\$x = Interpolate (\$L, \$Lt, \$xJe, \$xJt);
\$y = Interpolate (\$L, \$Lt, \$yJe, \$yJt);
} elsif (\$L <= \$Lt + \$Lm) {  # On middle segment
\$L = \$L - \$Lt;
\$x = Interpolate (\$L, \$Lm, \$xJt, \$xJf);
\$y = Interpolate (\$L, \$Lm, \$yJt, \$yJf);
} else {  # On frigid segment
\$L = \$L - \$Lt - \$Lm;
\$x = Interpolate (\$L, \$Lf, \$xJf, \$xP73);
\$y = Interpolate (\$L, \$Lf, \$yJf, \$yP73);
}

} # end zones (k) and (l)
} # end zones (j), (k), and (l)
} # end zones (i), (j), (k), and (l)
} # end all zones
return (\$x, \$y);
} # End of sub MPtoXY

sub LineIntersection {  # Written 2010-02-28; modified 2010-11-28
# Subroutine/function to calculate coordinates of point of intersection of two lines which
# are given by a point on the line and the line's slope angle in degrees.
#
# 2010-11-28 – Modified to assume that the lines do intersect, neither is either horizontal
#  or vertical, the arguments are the correct number and are all defined, and the
#  angles are within [-180,180].
#  Unlike on the previous version, this one has no checks and doesn't return a flag.
#
# Return is an array of two values, the x and y coordinates of the point of intersection.
#
# Arguments should be 6, in this order:
#  x and y coordinates of point of first line; slope of first line in degrees;
#  x and y coordinates of point of second line; slope of second line in degrees;
#
# Equations used are from: slope of line = tangent angle = delta-y / delta-x, and the fact
#  that intersection point x,y is on both lines.
use Math::Trig;
my (\$nArguments,\$xp,\$yp,\$m1,\$m2);
\$nArguments=@_ ;
my (\$x1,\$y1,\$angle1,\$x2,\$y2,\$angle2) = @_ ;

\$xp = (\$m1 * \$x1 - \$m2 * \$x2 - \$y1 + \$y2) / (\$m1 - \$m2);
\$yp = \$m1 * \$xp - \$m1 * \$x1 + \$y1;
return (\$xp,\$yp);
} # End of sub LineIntersection

sub Length {
# Input are x1,y1,x2,y2
my (\$x1,\$y1,\$x2,\$y2) = @_;
return sqrt( (\$x1-\$x2)**2 + (\$y1-\$y2)**2);
} # End of sub Length

sub Interpolate {
# Inputs are 4: length wanted, of total-segment-length, start, end;
# Total-segment-length is different from (end - start); end and start may be x-coordinates,
# or y-coordinates, while length takes into account the other coordinates.
# (End - Start ) / Length = (Wanted - Start) / NewLength
# Returns single value: Wanted.
my (\$NewLength, \$Length, \$Start, \$End) = @_;
my (\$Wanted);
\$Wanted = \$Start + (\$End - \$Start) * \$NewLength / \$Length;
\$Skip = "YES";
unless (\$Skip) {  # Option to skip printing arguments
foreach \$i (\$NewLength, \$Length, \$Start, \$End, \$Wanted) {printf "\t%.5f", \$i;}
print "\n";
} # End of skipping printing arguments and result
return \$Wanted;
}  # End of sub Interpolate

sub CircleLineIntersection {
# Subroutine to calculate intersection of circle with line segment. Equations from
# http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
# Arguments are 7, in the following order:
# Circle given as x-center, y-center, radius; line segment given as (x1,y1), (x2,y2).
# If line segment does not intersect circle, return is 0; else, return is 1,x,y of
# point of intersection; it is assumed that circle only intersects line segment at one
# point. If you want a subroutine for other purposes, read that website.
my (\$n);
\$n = @_;
if ( \$n != 7) {
print "Sub CircleLineIntersection requires 7 arguments but got \$n.\n";
return 0;
}
my (\$xc,\$yc,\$r,\$x1,\$y1,\$x2,\$y2) = @_;
my (\$u1,\$u2,\$a,\$b,\$c,\$d,\$x,\$y);
# Check if there is a point of intersection
\$a = (\$x2-\$x1)**2 + (\$y2-\$y1)**2;
\$b = 2 * ( (\$x2-\$x1) * (\$x1-\$xc) + (\$y2-\$y1) * (\$y1-\$yc) );
\$c = \$xc**2 + \$yc**2 + \$x1**2 + \$y1**2 - 2 * (\$xc*\$x1 + \$yc*\$y1) - \$r**2;
\$d = \$b**2 - 4*\$a*\$c;	# Determinant
if (\$a == 0) {
# print "In sub CircleLineIntersection: line given is just one point!\n";
return 0;
}elsif (\$d < 0) {
# Determinant is negative: circle does not intersect the line, much less the
# segment
# print "In sub CircleLineIntersection: line doesn't intersect circle.\n";
return 0;
}
# \$u1 and \$u2 are the roots to a quadratic equation
\$u1 = (-\$b + sqrt(\$d) ) / (2*\$a);	# + of +/- of the solution to the quadratic equation
\$u2 = (-\$b - sqrt(\$d) ) / (2*\$a);	# - of +/- of the solution to the quadratic equation

# Check if there is an intersection and if it is within the line segment (not only the line)
# If \$u1=\$u2, line is tangent to the circle; if \$u1 != \$u2, line intersects circle at two
# points; however, point or points of intersection are within the line segment only if
# the root is within interval [0,1].
if (0 <= \$u1 && \$u1 <= 1) {  # This root is on the line segment; use it
\$x = \$x1 + \$u1 * (\$x2 - \$x1);
\$y = \$y1 + \$u1 * (\$y2 - \$y1);
return 1,\$x,\$y;
} elsif (0 <= \$u2 && \$u2 <= 1) {
# 1st root was not on line segment but 2nd one is; use it
\$x = \$x1 + \$u2 * (\$x2 - \$x1);
\$y = \$y1 + \$u2 * (\$y2 - \$y1);
return 1,\$x,\$y;
} else {  # neither root is on line segment
# print "In sub CircleLineIntersection: line segment doesn't intersect circle.\n";
return 0;
}
}  # End of sub CircleLineIntersection

# -  -  -  -  -  -  -  -  -    SUBROUTINES For Coordinate conversion   -  -  -  -  -  -  -  -  -  -

sub LLtoMP {
# Arguments are real world longitude and latitude for one point, in decimal degrees.
# West longitudes and south latitudes have negative values.
# Returns corresponding meridian (\$m) and parallel (\$p) in MJ's template half-octant,
# sign for meridian (-1 for western half octant and +1 for eastern one), and octant
# number. Returned values of \$m and \$p are always positive.

my (\$Long, \$Lat) = @_ ;  # Input values

# \$m and \$p are the meridian and parallel numbers in template half-octant setting;
# \$Octant is the M-map octant of the real point; \$Sign is for east or west side of
# template octant;
my (\$m, \$p, \$Sign, \$Octant);  # Values to be returned

# @South are southern octants corresponding to northern octants 1, 2, 3 and 4; the 0
# is just a place holder to facilitate correspondence.
my (@South) = (0,6,7,8,5);  # Variables used only in this sub

# Determine the correct octant; Octant 1 is +160° to -110°; octant 4 is 70° to 160°
\$Octant = int ( (\$Long + 200) / 90 ) + 1;
# Make longitude fit within template half-octant, and determine if y value should
# be positive or negative.
\$m = ( (\$Long + 200) - (90*(\$Octant - 1))) - 45;
if (\$m < 0) {
\$Sign = -1;
\$m = -\$m;
} else {
\$Sign = 1;
}
# Fix the octant number, if necessary
if (\$Octant == 5) { \$Octant = 1; }
if (\$Lat < 0) {
\$Octant = \$South [\$Octant];
\$p = -\$Lat;
} else {
\$p = \$Lat;
}
return (\$m, \$p, \$Sign, \$Octant);
}  # End sub LLtoMP

sub MJtoG {
# Subroutine to convert (that is, do coordinate transformation of) x and y coordinates
# from Mary Jo's half-octant on its side to Gene's leaning, single octant coordinates, or
# to Gene's M-map (eight-octants) coordinates.
#
# Subroutine returns converted x and y coordinates.
#
# Arguments are:
# - x and y coordinates of point to convert.
# - Third argument is the Octant to convert to:
#   - 0 – Gene's single-octant system, with y-axis on its left;
#   - 1, 2, 3 or 4 – convert to M-map coordinates, respectively to first, second, third or
#     fourth northern octant, from the left;
#   - 5, 6, 7 or 8 – convert to M-map coordinates, respectively to fourth, first, second or
#     third southern octant from the left.
# - \$sin60, \$cos60, \$yTranslate – values calculated once, in sub Interpolate, to minimize
#     computations. (\$sin60 = sin 60°, \$cos60 = cos 60°, \$yTranslate = 10,000 * sin 60°).
#
# - In MJ's coordinates, point M is the origin, at (0,0), points M, A and G are on the positive
# x-axis, and point G is at (10000, 0).
# - In G's system, point M, L, J and P are on the positive y-axis, and point G is on the
# positive x-axis; in this system, point G is at coordinates (5000, 0).
# - From MJ's system to G's, there is a +60° rotation, and also a translation.
# - The M-map coordinate system is like G's system, except that the y-axis is 10000mm
# to the right, that is, the x-coordinates for the start octant are 10000mm smaller.
#
# I got the equations for rotation and translation from my pocketbook "The Universal
# Encyclopedia of Mathematics, with a Foreword by James R. Newman", ©1964 by
# George Allen and Unwin, Ltd.; translated from original German language edition,
# pages 152, 153.

my (\$nArgs, \$xnew, \$ynew);
my (\$x, \$y, \$Octant, \$sin60, \$cos60, \$yTranslate) = @_ ;
if (not defined (\$Octant) ) { \$Octant = 0; }
if (\$Octant == 0) {
(\$xnew, \$ynew) = Rotate (\$x, \$y, 60, \$sin60, \$cos60);
} elsif (\$Octant == 1) {
(\$xnew, \$ynew) = Rotate (\$x, \$y, 120, \$sin60, \$cos60);
\$xnew = \$xnew - 10000;
} elsif (\$Octant == 2) {
(\$xnew, \$ynew) = Rotate (\$x, \$y, 60, \$sin60, \$cos60);
\$xnew = \$xnew - 10000;
} elsif (\$Octant == 3) {
(\$xnew, \$ynew) = Rotate (\$x, \$y, 120, \$sin60, \$cos60);
\$xnew = \$xnew + 10000;
} elsif (\$Octant == 4) {
(\$xnew, \$ynew) = Rotate (\$x, \$y, 60, \$sin60, \$cos60);
\$xnew = \$xnew + 10000;
} elsif (\$Octant == 5) {
(\$xnew, \$ynew) = Rotate ((20000-\$x), \$y, 60, \$sin60, \$cos60);
\$xnew = \$xnew + 10000;
} elsif (\$Octant == 6) {
(\$xnew, \$ynew) = Rotate ((20000-\$x), \$y, 120, \$sin60, \$cos60);
\$xnew = \$xnew - 10000;
} elsif (\$Octant == 7) {
(\$xnew, \$ynew) = Rotate ((20000-\$x), \$y, 60, \$sin60, \$cos60);
\$xnew = \$xnew - 10000;
} elsif (\$Octant == 8) {
(\$xnew, \$ynew) = Rotate ((20000-\$x), \$y, 120, \$sin60, \$cos60);
\$xnew = \$xnew + 10000;
} else {
print "Error converting to M-map coordinates; there is no \$Octant octant!\n";
return (\$x,\$y);
}
\$ynew = \$ynew + \$yTranslate;
return (\$xnew, \$ynew);

}  # End of sub MJtoG, which converts coordinates to octants on M-map

sub Rotate {
# Receives 5 arguments: x, y, angle by which to rotate the coordinate system, and
# sin 60° and cos 60°. The last two are calculated once in sub Preliminary, to minimize
# computations.
# Expects that the axes will be rotated either 60° or 120°.
# Returns new x and y values.
my (\$x, \$y, \$angle, \$sin60, \$cos60) = @_ ;
my (\$xnew, \$ynew);
if ( \$angle == 60 ) {
\$xnew = \$x * \$cos60 + \$y * \$sin60;
\$ynew = -\$x * \$sin60 + \$y * \$cos60;
} elsif ( \$angle == 120 ) {
\$xnew = -\$x * \$cos60 + \$y * \$sin60;
\$ynew = -\$x * \$sin60 - \$y * \$cos60;
} else {
print "Sub Rotate expected angle = 60 or 120 but received \$angle.!\n";
}
return \$xnew, \$ynew;
}  # End of sub Rotate
# -  -  -  -  -  -  -  -  -       End SUBROUTINES         -  -  -  -  -  -  -  -  -  -  -  -  -```