Go back to B.J.S. Cahill Resource Page
Cahill Butterfly Map 1909
Cahill 1909
Go back to Gene Keyes home page

Cahill-Keyes M-layout world map silhouette including Antarctica
Cahill-Keyes 1975

Page 2 of 6

Cahill-Keyes Octant Graticule:
Principles and Specifications

with Perl programs and OpenOffice.org 2.0 macros
for 1/1,000,000 Megamap


Gene Keyes
2010-08-20


9) Constructing the Cahill-Keyes Octant Graticule
and Megamap with Perl programming
and OpenOffice.org Draw macros



Notes by Mary Jo Graça:


I started programming a whole octant, using Gene Keyes' coordinate system, that is, with a scaffold triangle upright and tilted left in a 10,000 x 12,000 grid. But I decided that having to consider both positive and negative meridians and positive and negative angles of lines was much more complicated than necessary. So I chose to work with a different coordinate system and with only half the octant. I named this, MJ's coordinate system (MJ for Mary Jo), in an area 10,000 x 5,773.5027 (the coordinates of point N).

Only the positive half of the octant is calculated. This is the "template" half octant. It is laid on its side, with point M at the origin of the set of axes, and point G on the x-axis, at point (0; 10,000):

Fig. 6. Half-octant as used by Mary Jo. This is merely a jpeg; lines in a printed version from the original odg or pdf would not be as wavy. GK
half-octant

(by MJ, cont.)

This has several advantages:

• All the angles are measured in a counterclockwise direction from the positive x-axis; this being the way I am used to thinking of angles, sines and cosines, there is less opportunity for mistakes.

• All the angles for meridians are positive. This makes intersection of the meridian segments easier to calculate

• To calculate the other half of the octant, one has only to use the negative of the y-coordinates

• Calculation of the coordinates for the southern octant are also easy: for a given latitude-longitude, the x-coordinate is 20,000 minus the x-coordinate for the northern octant, and the y-coordinate equals that of the northern octant


Given the above, all calculations for meridian joints, and meridian-parallel intersection can be done for the positive half-octant, and then easily converted for the remaining half or twin octant.

When coordinates are desired for Gene's one-octant or the 8 octants of the world map this is the procedure:

• For a land point, it is determined which octant it is on, based on its latitude and longitude; then, the absolute value of the latitude is taken, and the longitude is converted to the corresponding meridian on the "template half-octant"

• Using these template latitude and longitude, the point's position in MJ's coordinates is calculated

• If the point is on the negative side or the southern octant, its coordinates are converted to have the correct x and y coordinates

• Then a formula is applied to rotate and translate from MJ's coordinate system to Gene's, based on its octant

Notes by Gene Keyes

"HalfOctant8", the first of two Perl programs by Mary Jo, outputs a complete data set of all the x-y coordinates for all one degree geocells of a horizontal half octant (including the very narrow polar geocells from 90 to 85°), plus meridian joints, octant border, etc. The data format is a text-file of tab-separated values [shown as .csv, comma-separated values], which is convertible to a spreadsheet. That complete table is presented here in HTML on page 3. Its x-y coordinates differ from my hand-calculated, whole-octant values, but the segment lengths and geocell dimensions remain the same.

Then on p. 4 is Mary Jo's second Perl program, "OOmacroMaker4", to actually prepare macros which can draw the Cahill-Keyes Octant Graticule, and Megamap, in OpenOffice.org, using the output of this first Perl program.


10) Perl Program "HalfOctant8" by Mary Jo Graça
to calculate one-degree coordinates
for Cahill-Keyes Octant Graticule Template
(including scaffold triangle and octant itself)


 (Done as a half-octant with MJ's coordinates; see explanation above)



#! /usr/bin/perl -w
# Program HalfOctant8. Modified from program HalfOctant7 on 2010-04-18.
# Program outputs all the coordinates as a text file with tab separated values
# (csv format) which can be read by a spreadsheet program such as Excel or
# OpenOffice.org.
#
# Parts relating to TurboCAD macros, including functions for coordinate conversion,
# were deleted on 2010-08-05, since we've given up on using our old version of TurboCAD.
#
# Parallels in this version are as follows:
# - Latitude 75° to 89° — parallels are circular arcs with center at point A; parallel spacing is
# 104.
# - Latitude 74° — circular arc, with center at A and radius 100 more than parallel 75°, from
# meridian 0° to meridian 30°; half way between parallels 73° and 75° from
# meridians 31° to 45°.
# - Latitude 73° — circular arc, with center at A and radius 100 more than parallel 74°, from
# meridian 0° to meridian 30°; straight line from meridian 30° to meridian 45°,
# perpendicular to line BD.
# - Latitudes 1° to 72°, meridians 0° to 29° — equally spaced along meridian from equator to
# 73° of latitude.
# - Latitude 15°, meridian 45° — at the tropic joint, that is, at point D.
# - Latitude 15°, meridians 29° to 45° — circular arc; center chosen on the extension of line
# M-B-D-N, the side of the octant (corresponding to meridian 45°). Center is
# determined by: drawing chord from meridian 29° to 45°; drawing a line
# perpendicular to this chord, through its mid-point; finding intersection with
# line through point M (0,0) with angle of 30° (that is, with line M-B-D-N).
# - Latitudes 1° to 14°, meridians 30° to 45° — equally spaced along meridian, from equator
# to parallel 15°.
# - Latitudes 16° to 72°, meridians 30° to 45° — equally spaced along meridian, from parallel
# 15° to parallel 73°.
#
# Multidimensional arrays are done with hashes, with the keys used for i,j indexes.

use Math::Trig; # Call for library of trigonometric functions

# Given input
$xM = 0; $yM = 0; # Point M is the origin of the axes
$MA = 940; # Indent of octant's apex, with respect to triangle's
$xG = 10000; $yG = 0; # Point G, at center of base of octant
$DeltaFrigid = 104; # Distance for each degree of latitude close to the pole

# Other points and lengths of special interest
$xA = $MA; $yA = 0; # Point A, at the indent
$xN = $xG; $yN = $xG * tan (deg2rad(30)); # Point N
($flag, $xB, $yB) = LineIntersection ($xM,$yM,30,$xA,$yA,45); # Point B
if ($flag ne "Intersect") { die "$flag while calculating point B.\n"; }
$AG = $xG - $xA;
$AB = Length ($xA, $yA, $xB, $yB);
$EF = $AB; # Sometimes I prefer AB, others EF
$MB = Length ($xM, $yM, $xB, $yB);
$MN = Length ($xM, $yM, $xN, $yN);
$xD = Interpolate ($MB, $MN, $xN, $yM); # D is away from N as B is away from M
$yD = Interpolate ($MB, $MN, $yN, $yM);
$xF = $xG;
$yF = $yN - $MB;
$xE = $xN - $MA * sin (deg2rad(30));
$yE = $yN - $MA * cos (deg2rad(30));
$FG = $yF;
$BD = Length ($xB, $yB, $xD, $yD);
$EFG = $EF + $FG;

# Calculate polar start, two joints, and equatorial end, for each meridian; hashes %xJ, %yJ.
# Also calculate a few other values for other hashes, on the way.

$DeltaMEq = $EFG / 45; # 45 meridian spacings along equator for half an octant
$DeltaX5 = 5 * $DeltaFrigid;

foreach $m (0 .. 45) { # Loop for each meridian.
$M = $m . ',' ; # part of hash key (array row number)

# Assign Frigid distance of parallels per degree of latitude.
$dP {$M . 0} = $DeltaFrigid;

# Start point of meridians
if ($m % 5 == 0) { # Every fifth meridian starts at A
$xJ {$M . 0} = $xA; $yJ {$M . 0} = $yA;
} else { # Minor meridians start at 85° lat, that is 5° from A
$xJ {$M . 0} = $xA + $DeltaX5 * cos (deg2rad($m));
$yJ {$M . 0} = $yA + $DeltaX5 * sin (deg2rad($m));
}

# End point of meridians, that is, equidistant intersections with equator;
# These are also the meridian crossings for the zero-th parallel
$distance = $m * $DeltaMEq;
if ($distance <= $yF) {
$xJ {$M . 3} = $xP {$M . 0} = $xG;
$yJ {$M . 3} = $yP {$M . 0} = $distance;
} else { # Past joint at point F
$extra = $distance - $yF;
$xJ {$M . 3} = $xP {$M . 0} = Interpolate ($extra, $EF, $xF, $xE);
$yJ {$M . 3} = $yP {$M . 0} = Interpolate ($extra, $EF, $yF, $yE);
}

# Calculate both joints of meridians

if ($m == 0) { # Straight meridian at center of octant
# This is a good place to calculate a bunch of things for center line, AG.

# No joints; give coordinates of point G for the joints
@xJ {$M . 1, $M . 2} = ($xG, $xG);
@yJ {$M . 1, $M . 2} = ($yG, $yG);
# Length of meridian segments; this whole meridian is a single segment
@L {$M . 0, $M . 1, $M . 2, $M . 3, $M . 4} = ($AG, 0, 0, $AG, $AG);
# May as well calculate parallel intersections along this straight line
foreach $p (75 .. 89) { # (foreach doesn't like first value greater than last)
$xP {$M . $p} = $xA + (90 - $p) * $DeltaFrigid;
$yP {$M . $p} = $yA;
}
# Remaining parallels are equidistant over remaining meridian length
$delta = ($xG - $xP {$M . 75}) / 75;
foreach $p (1 .. 74) { # Parallel 0° was done with end point of meridians
$xP {$M . $p} = $xG - ($p * $delta);
$yP {$M . $p} = $yG;
}
# Distance between parallels
@dP {$M . 1, $M . 2, $M . 3} = ($delta, $delta, $delta);
@L {$M . 5} = $AG - $xP {$M . 73};

} else { # Meridians other than 0°.
# Calculate joints by line intersection from A and M
($flag, $xJ {$M . 1}, $yJ {$M . 1}) =
LineIntersection ($xM, $yM, 2*$m/3, $xA, $yA, $m);
if ($flag ne "Intersect") { die "$flag while calculating 1st joint at $m.\n"; }
($flag, $xJ {$M . 2}, $yJ {$M . 2}) =
LineIntersection ($xM, $yM, 2*$m/3, $xJ {$M . 3}, $yJ {$M . 3}, $m/3);
if ($flag ne "Intersect") { die "$flag while calculating 2st joint at $m.\n"; }
# Calculate length of meridian segments: These are calculated from point A
# to equator, even if this meridian starts at parallel 85°.
$L {$M . 0} = Length ($xA, $yA, $xJ {$M . 1}, $yJ {$M . 1});
$L {$M . 1} = Length ($xJ {$M . 1}, $yJ {$M . 1}, $xJ {$M . 2}, $yJ {$M . 2});
$L {$M . 2} = Length ($xJ {$M . 2}, $yJ {$M . 2}, $xJ {$M . 3}, $yJ {$M . 3});
$L {$M . 3} = $L {$M . 1} + $L {$M . 2}; # Torrid + Temperate
$L {$M . 4} = $L {$M . 3} + $L {$M . 0}; # Total meridian length from A
} # End of calculating meridian joints
} # End of loop for each meridian.

# Do each parallel.
# Do parallels for Frigid region: Parallel 75° to 89°
foreach $p (75 .. 89) {
foreach $m (1 .. 45) {
$M = $m . ',' ;
$distance = (90 - $p) * $dP {$M . 0}; # $dP {$M . 0} = $DeltaFrigid
$xP {$M . $p} = Interpolate ($distance, $L {$M . 0}, $xA, $xJ {$M . 1});
$yP {$M . $p} = Interpolate ($distance, $L {$M . 0}, $yA, $yJ {$M . 1});
}
} # End of Frigid region parallels
# Do supple zone: Parallels 73° and 74°
$distance75 = (90 - 75) * $DeltaFrigid; # Distance from pole to parallel 75
foreach $p (73 .. 74) { # Do only up to meridian 30°; rest is supple zone
foreach $m (1 .. 30) {
$M = $m . ',' ;
# For these meridians, the "supple zone" parallel distance is the same as
# for meridian 0°, since these parallels are circular arcs for this stretch.
$dP {$M . 1} = $dP {0 . ',' . 1};
$distance = $distance75 + (75 - $p) * $dP {$M . 1};
$xP {$M . $p} = Interpolate ($distance, $L {$M . 0}, $xA, $xJ {$M . 1});
$yP {$M . $p} = Interpolate ($distance, $L {$M . 0}, $yA, $yJ {$M . 1});
if ($p == 73) {
$L {$M . 5} = $L {$M . 4} - $distance; # Distance 0° to 73°
}
}
}
# Find intersection meridian 45° with parallel 73°, by a straight line from meridian 30°,
# perpendicular to line BD (or MN). Then calculate supple distances and rest of parallels
# 73° and 74°
$M = 45 . ',' ;
($flag, $xP {$M . 73}, $yP {$M . 73} ) =
LineIntersection ($xP {30 . ',' . 73}, $yP {30 . ',' . 73}, -60, $xB, $yB, 30);
if ($flag ne "Intersect") { die "$flag while calculating 45,73.\n"; }

# Supple zone parallel distance is half of parallel 73° to 75°; at meridian 45°, it goes over
# a joint at point B.
$L {$M . 5} = $L {$M . 4} - Length ($xP {$M . 73}, $yP {$M . 73}, $xB, $yB) - $AB;
$dP {$M . 1} = ( Length ($xP {$M . 73}, $yP {$M . 73}, $xB, $yB) +
Length ($xB, $yB, $xP {$M . 75}, $yP {$M . 75}) ) / 2;
$xP {$M . 74} = Interpolate ( ($distance75 + $dP {$M . 1}), $AB, $xA, $xB);
$yP {$M . 74} = Interpolate ( ($distance75 + $dP {$M . 1}), $AB, $yA, $yB);
foreach $m (31 .. 44) {
$M = $m . ',' ;
# Intersect meridian with straight parallel line from meridian 30° (to 45°)
($flag, $xP {$M . 73}, $yP {$M . 73}) =
LineIntersection ($xP {30 . ',' . 73}, $yP {30 . ',' . 73}, -60, $xA, $yA, $m);
if ($flag ne "Intersect") { die "$flag while calculating $m,73.\n"; }
$L {$M . 5} = $L {$M . 4} - Length ($xP {$M . 73}, $yP {$M . 73}, $xA, $yA);

# Supple zone parallel distance is half of parallel 73° to 75°; at these meridians it
# does not go over a meridian joint; this only happens on meridian 45° (or some
# meridians greater than 44°, if we did them for less than 1° interval).
$dP {$M . 1} = Length($xP {$M.73}, $yP {$M.73}, $xP {$M.75}, $yP {$M.75}) / 2;
$xP {$M . 74} = Interpolate (($distance75+$dP{$M.1}), $L{$M.0}, $xA, $xJ {$M.1});
$yP {$M . 74} = Interpolate (($distance75+$dP{$M.1}), $L{$M.0}, $yA, $yJ {$M.1});
}
# Calculate parallel crossings from 1° to 72° of latitude, for meridians 1° to 29°, plus 45°;
# meridians 30° to 44° are dealt with separately.
foreach $m (1 .. 29) {
$M = $m . ',' ;
# For these meridians, the distance between parallels is the same for both the
# Temperate and Torrid zones, and is 1 / 73 [for parallels 0° to 73°] of
# (total length - distance pole to parallel 75° - twice the supple area parallel
# distance) [twice the supple zone is for parallel 75° to 73°].
# In other words, 1 / 73 of distance from equator to parallel 73°, along meridian line.
$dP {$M . 2} = $dP {$M . 3} = ($L {$M . 4} - $distance75 - 2 * $dP {$M . 1}) / 73;
foreach $p (1 .. 72) {
$distance = $p * $dP {$M . 2};
if ($distance <= $L {$M . 2} ) { # Torrid segment
$xP{$M.$p} = Interpolate ($distance, $L{$M.2}, $xJ{$M.3}, $xJ{$M.2});
$yP{$M.$p} = Interpolate ($distance, $L{$M.2}, $yJ{$M.3}, $yJ{$M.2});
} elsif ($distance <= $L {$M . 3} ) { # Temperate segment, past one joint
$extra = $distance - $L {$M . 2};
$xP{$M.$p} = Interpolate ($extra, $L{$M.1}, $xJ{$M.2}, $xJ{$M.1});
$yP{$M.$p} = Interpolate ($extra, $L{$M.1}, $yJ{$M.2}, $yJ{$M.1});
} else { # Frigid segment, past two joints
$extra = $distance - $L {$M . 3};
$xP{$M.$p} = Interpolate ($extra, $L{$M.0}, $xJ{$M.1}, $xA);
$yP{$M.$p} = Interpolate ($extra, $L{$M.0}, $yJ{$M.1}, $yA);
}
} # End of foreach $p loop
} # End of parallels from latitude 1° to 72°, between meridians 1° and 29°

# Calculate parallel crossings at meridian 45°: point D is parallel 15°; parallels are equally
# spaced from 0° to 15° latitude, and also equally spaced from 15° to 73° latitude.
$M = 45 . ',' ;
$xP {$M . 15} = $xD; $yP {$M . 15} = $yD;
$dP {$M . 2} = Length ($xD, $yD, $xP {$M . 73}, $yP {$M . 73} ) / (73 - 15);
$dP {$M . 3} = $EF / 15; # $EF = length E to F = length E to D
foreach $p (1 .. 15) {
$distance = $p * $dP {$M . 3};
$xP {$M . $p} = Interpolate ($distance, $EF, $xE, $xD);
$yP {$M . $p} = Interpolate ($distance, $EF, $yE, $yD);
}
foreach $p (16 .. 72) {
$distance = ($p - 15) * $dP {$M . 2};
$xP {$M . $p} = Interpolate ($distance, $L {$M . 1},$xD, $xB);
$yP {$M . $p} = Interpolate ($distance, $L {$M . 1},$yD, $yB);
}

# Parallel 15°, and parallels 1° to 72°, between meridians 29° and 45°.
# Parallel 15° is a circular arc from meridian 29° to meridian 45°; center of circle on x-axis.
$P = ',15';
# Find center of chord for parallel 15°, between meridians 29° and 45°
$xmid = ($xP{29 . $P} + $xP{45 . $P}) / 2;
$ymid = ($yP{29 . $P} + $yP{45 . $P}) / 2;
# Slope of this chord is delta-y/delta-x; slope of perpendicular line is -delta-x/delta-y
$mmid = - ($xP{29 . $P} - $xP{45 . $P}) / ($yP{29 . $P} - $yP{45 . $P});
$midAngle=rad2deg(atan($mmid));
# Choose center to be point at intersection of this line with side of octant (that is, side
# which has angle 30° and corresponds to meridian 45°). [ ($xM,$yM) = (0,0) ]
($flag,$xC,$yC) = LineIntersection ($xM,$yM,30,$xmid,$ymid,$midAngle);
if ($flag ne "Intersect") { die "$flag while calculating center of circle for P 15°.\n"; }
$R = Length($xC,$yC,$xP{29 . $P},$yP{29 . $P}); # Radius of circle
$R2 = Length($xC,$yC,$xP{45 . $P},$yP{45 . $P}); # Radius for comparison; checks
print "Chord has ends at ",$xP{29 . $P},",",$yP{29 . $P}," and ",$xP{45 . $P}, ",";
print $yP{45 . $P},"\nMid-point of chord is at $xmid,$ymid;\n";
print "Center of circle is at $xC,$yC;\n";
print "Slope of ray through mid-point is $mmid and corresponding angle is $midAngle.\n";
print "Radius of circular arc is $R or $R2.\n";

# Do all the calculations of intersections of parallels with meridians 30° to 44°
foreach $m (30 .. 44) {
$M = $m . ',';
# Find intersection of circle with meridian; try torrid segment first
($flag,$xP{$m.$P},$yP{$m.$P}) = CircleLineIntersection($xC,$yC,$R,
$xJ{$M . 3},$yJ{$M . 3},$xJ{$M . 2},$yJ{$M . 2});
if ($flag==1) { # found an intersection point
# Distance Equator to parallel 15°.
# warn "Parallel 15° at meridian $m is in torrid zone.\n";
$dist15 = Length ($xP{$m.$P},$yP{$m.$P},$xJ{$M . 3},$yJ{$M . 3});
$dP{$M . 3} = $dist15 / 15;
$dP{$M . 2} = ($L{$M . 5} - $dist15) / 58;
} else { # intersection is on temperate segment
($flag,$xP{$m.$P},$yP{$m.$P}) = CircleLineIntersection($xC,$yC,$R,
$xJ{$M . 2},$yJ{$M . 2},$xJ{$M . 1},$yJ{$M . 1});
if ($flag==0) { # Hmmm.... no intersection!
print "No line-circular arc intersection for M $m, at parallel 15!";
die;
}
# Distance Equator to parallel 15°.
$dist15=Length($xP{$m.$P},$yP{$m.$P},$xJ{$M. 2},$yJ{$M. 2})+$L{$M. 2};
$dP{$M . 3} = $dist15 / 15;
$dP{$M . 2} = ($L{$M . 5} - $dist15) / 58;
}
# Calculate parallel crossings for latitudes 1° to 14°.
# Note that there might be a joint below parallel 15°.
foreach $p (1 .. 15) {
$distance = $p * $dP {$M . 3};
if ($distance <= $L {$M . 2} ) { # Torrid segment
$xP{$M.$p} = Interpolate ($distance, $L{$M.2}, $xJ{$M.3}, $xJ{$M.2});
$yP{$M.$p} = Interpolate ($distance, $L{$M.2}, $yJ{$M.3}, $yJ{$M.2});
} else { # Temperate segment, past one joint
$extra = $distance - $L {$M . 2};
$xP{$M.$p} = Interpolate ($extra, $L{$M.1}, $xJ{$M.2}, $xJ{$M.1});
$yP{$M.$p} = Interpolate ($extra, $L{$M.1}, $yJ{$M.2}, $yJ{$M.1});
}
} # End of foreach $p loop for parallels from latitude 1° to 15°
# Calculate parallel crossings for latitudes 16° to 72°, from meridian 30° to 44°;
# they may or may not all be above Tropical joint on the meridian.
foreach $p (16 .. 72) {
$distance = ($p - 15) * $dP {$M . 2} + $dist15; # Distance from Equator
if ($distance <= $L {$M . 2} ) { # Torrid segment
$xP{$M.$p} = Interpolate ($distance, $L{$M.2}, $xJ{$M.3}, $xJ{$M.2});
$yP{$M.$p} = Interpolate ($distance, $L{$M.2}, $yJ{$M.3}, $yJ{$M.2});
} elsif ($distance <= $L {$M . 3} ) { # Temperate segment, past one joint
$extra = $distance - $L {$M . 2};
$xP{$M.$p} = Interpolate ($extra, $L{$M.1}, $xJ{$M.2}, $xJ{$M.1});
$yP{$M.$p} = Interpolate ($extra, $L{$M.1}, $yJ{$M.2}, $yJ{$M.1});
} else { # Frigid segment, past two joints
$extra = $distance - $L {$M . 3};
$xP{$M.$p} = Interpolate ($extra, $L{$M.0}, $xJ{$M.1}, $xA);
$yP{$M.$p} = Interpolate ($extra, $L{$M.0}, $yJ{$M.1}, $yA);
}
} # End of foreach $p loop for parallels 16° to 72°
} # End of parallels from latitude 1° to 72°, at meridians 30° and 44°

# Output all the values in spreadsheet readable csv format (actually, it is tab separated
# format instead of comma separated format).

open (CSV, ">Macros8\/Hashes8.csv") ||
die "Sorry, I couldn\'t create " . ">Macros8\/Hashes8.csv" . "\n ";
print CSV "\ndP\n"; PrintHash (3,%dP);
print CSV "\nL\n"; PrintHash (5,%L);
print CSV "\nxJ\n"; PrintHash (3,%xJ);
print CSV "\nyJ\n"; PrintHash (3,%yJ);
print CSV "\nxP\n"; PrintHash (89,%xP);
print CSV "\nyP\n"; PrintHash (89,%yP);
# Print coordinates of points of half-octant and scaffold triangle
print CSV "\nPoints";
foreach $i qw(A B D E F G M N Center) {printf CSV "\t%s",$i; }
print CSV "\nx";
foreach $i ($xA,$xB,$xD,$xE,$xF,$xG,$xM,$xN,$xC) {printf CSV "\t%11.4f", $i};
print CSV "\ny";
foreach $i ($yA,$yB,$yD,$yE,$yF,$yG,$yM,$yN,$yC) {printf CSV "\t%11.4f", $i};
# Print lengths of line segments of half-octant and scaffold triangle
print CSV "\n\nLengths";
foreach $i qw(AB BD DE EF FG AG MG MN GN MA MB EFG ABDE R)
{printf CSV "\t%s",$i; }
print CSV "\n";
# length DE = length AB = length EF; also, length GN = y of point N
foreach $i ($AB,$BD,$EF,$EF,$FG,$AG,$xG,$MN,$yN,$MA,$MB,$EFG,2*$EFG,$R)
{printf CSV "\t%11.4f", $i}; print CSV "\n";
close (CSV);
print 'Wrote values of arrays to file "Macros8/Hashes8.csv".', "\n";


# - - - - - - - - - SUBROUTINES - - - - - - - - - - - - -
sub PrintHash {
# Prints one of my two-dimension arrays in the form of a hash; it is assumed that
# the indexes are separated by a comma, the first index varies from 0 to 45 and the
# second index varies from 0 to the first input parameter.
# Values are printed to filehandle CSV.
my ($nI, %Hash) = @_;
my ($i,$m);
foreach $i (0 .. 45) { printf CSV "\t%10d",$i; } print CSV "\n";
foreach $i (0 .. $nI) {
print CSV $i;
foreach $m (0 .. 45) {
if (defined ($Hash{$m.",".$i})) {printf CSV "\t%11.4f",$Hash{$m.",".$i};
} else { printf CSV "\t undef"; }
}
print CSV "\n";
}
} # End of sub PrintHash

sub LineIntersection { # Written 2010-02-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.
#
# Return is an array of either one or three values:
# - if the lines don't intersect or are the same, the return is "Parallel";
# - if the lines intersect, the returns are: "Intersect" to mean that it worked, and
# the x and y coordinates of the point of intersection.
# - if not enough arguments were given or the angles were too large, then return is
# "Error".
#
# 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) = @_ ;

if ($nArguments != 6) { # Not the exact number of arguments
print "Sub LineIntersection requires 6 arguments but received $nArguments \n";
return "Error";
} elsif ($angle1<-180 || $angle2<-180 || $angle1>180 || $angle2>180) { # Angles too large
print "Sub LineIntersection: Angles should be between -180 and 180. \n";
return "Error";
} elsif ($angle1==$angle2 || $angle1 == ($angle2+180) || $angle1 == ($angle2-180) ) {
return "Parallel"; # Lines are either parallel or the same line
} elsif ( not defined($x1) || not defined($y1) || not defined($angle1) ||
not defined($x2) || not defined($y2) || not defined($angle2) ) {
print ("Undefined in LineInterpolation: $x1,$y1,$angle1,$x2,$y2,$angle2\n");
return "Error";
} else { # Lines intersect. Calculate point of intersection.
if ($angle1 == 0 || $angle1 == 180 || $angle1 == -180) { # Line 1 horizontal
if ($angle2 == 90 || $angle2 == -90) { # Line 1 horizontal, line 2 vertical
return ("Intersect",$y1,$x2);
} else { # Line 1 horizontal, line 2 slanted
$xp = $x2 - ($y2 - $y1) / tan(deg2rad($angle2));
return ("Intersect",$xp,$y1);
}
} elsif ($angle1 == 90 || $angle1 == -90) { # Calculate when line 1 is vertical
if ($angle2 == 0 || $angle2 == 180 || $angle2 == -180) {
# Line 1 vertical, line 2 horizontal
return ("Intersect",$x1,$y2);
} else { # Line 1 vertical, line 2 slanted
$yp = $y2 - ($x2 - $x1) * tan(deg2rad($angle2));
return ("Intersect",$x1,$yp);
}
} elsif ($angle2 == 0 || $angle2 == 180 || $angle2 == -180) {
# Line 1 slanted, line 2 horizontal
$xp = $x1 - ($y1 - $y2) / tan(deg2rad($angle1));
return ("Intersect",$xp,$y2);
} elsif ($angle2 == 90 or $angle2 == -90) { # Line 1 slanted, line 2 vertical
$yp = $y1 - ($x1-$x2) * tan(deg2rad($angle1));
return ("Intersect",$x2,$yp);
} else { # Line 1 slanted, line 2 slanted
$m1 = tan(deg2rad($angle1));
$m2 = tan(deg2rad($angle2));
$xp = ($m1 * $x1 - $m2 * $x2 - $y1 + $y2) / ($m1 - $m2);
$yp = $m1 * $xp - $m1 * $x1 + $y1;
return ("Intersect",$xp,$yp);
} #End of all the vertical, horizontal, slanted combinations
} # End of if-else — lines intersect
} # 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.
# (Start - End ) / Length = (Wanted - Start) / NewLength
# Returns single value: Wanted.
my ($NewLength, $Length, $Start, $End) = @_;
my ($Wanted);
$Wanted = $Start + ($End - $Start) * $NewLength / $Length;
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

Additional output:

Chord has ends at 8396.03286341536,2945.73365712891 and 7775.93612044263,4489.43881233888
Mid-point of chord is at 8085.984491929,3717.5862347339;
Center of circle is at 2672.81110719657,1543.14821223296;
Slope of ray through mid-point is 0.401693769616511 and corresponding angle is 21.885020809006.
Radius of circular arc is 5892.58120021185 or 5892.58120021185.
Wrote values of arrays to file "Macros8/Hashes8.csv".

Go to page 3, HTML Table of Half-Octant Coordinates


Go back to page 1
Go back to B.J.S. Cahill Resource Page

Go back to Gene Keyes home page