// Solution for the "WWII" geocache (GCYANP)
// http://www.geocaching.com/seek/cache_details.aspx?guid=a94cdfa3-a467-40e2-974c-37cb57adba95
// "WAYPOINT #1
// There is a 14k+ mountain peak in CO that is named the same as the
// battleship where the Japanese signed their WWII surrender to the US."
// "If you stood on the summit of this mountain at 7pm UT on the day of the
// Japanese surrender signing, you would observe the Sun at an altitude of
// ____ degrees and at an azimuth (E of N) of ______ degrees."
//
// The ship will be trivial to anyone who has seen the fabulous Steven Seagal
// epic "Under Siege".
//
// Location of the summit of the mountain
// http://www.14ers.com/php14ers/qpick.php?parmpeak=36
summitLat =DMS[38,56,50]
summitLong=DMS[106,22,40]
summitHeight=14067 feet
// Grab in my sun/moon prediction library
use sun.frink
surrenderDate = #September 2, 1945 7:00 PM UT#
println["Local time is " + (surrenderDate -> Mountain)]
// Calculate refracted, parallax-corrected apparent position of sun.
[azimuth, altitude] = refractedSunAzimuthAltitude[surrenderDate, summitLat, summitLong]
// Convert Meeus' odd coordinate system to normal coordinates
azimuth = (azimuth + 180 degrees) mod circle
println["Altitude is: " + format[altitude, "degrees", 5]]
println["Azimuth is: " + format[azimuth, "degrees", 5]]
// Meta-Calculation:
// It's hard to guess what temperature and atmospheric pressure they assumed
// for the day, so use my defaults. The calculations show that the sun
// would be at an altitude of about 59 degrees (with refraction) so refraction
// error should hopefully be sort of low. Let's calculate it here.
println["\nMeta-Calculation of possible refraction discrepancy:"]
[airlessAz, airlessAlt] = airlessSunAzimuthAltitude[surrenderDate, summitLat, summitLong]
// Convert from Meeus odd coordinate system
airlessAz = (airlessAz + 180 degrees) mod circle
println["Airless Altitude is: " + format[airlessAlt, "degrees", 5]]
println["Airless Azimuth is: " + format[airlessAz, "degrees", 5]]
// Calculate refraction angle using my defaults.
refractionAngle = refractionAngle[airlessAlt]
println[]
println["Refraction angle is: " + format[refractionAngle, "degrees", 5]]
println["If they didn't correct for refraction, error is " + format[refractionAngle earthradius, "feet", 1] + "."]
// Running this shows that the refraction angle is small (about 0.01 degrees)
// so it shouldn't be a problem if we round to the nearest degree, but it is
// possibly a problem if we round to the nearest 0.1 degree, and certainly
// if we round to the nearest 0.01 degree.
// The problem solution doesn't state if/how they rounded the value. This is
// potentially a show-stopper.
// After more research, the USNO page seems to show one decimal place
// after the decimal point for altitude and azimuth. Can we assume that's
// what they were using? It's too big of an assumption to make carelessly;
// .1 degree on the earth's surface is (in Frink notation:)
// .1 degree earthradius -> miles
//
// that's still almost 7 miles! Far too large an area for me to hike.
// In addition, it's hard to know if they really corrected for the parallax
// of the summit, as opposed to just using the reference geoid. The problem
// statement says "from the summit," so let's assume so. Unfortunately,
// my parallax model only currently calculates from the reference geoid,
// so let's see how bad that might throw us off.
// Based on the link they gave to the USNO, http://aa.usno.navy.mil/ ,
// which doesn't have inputs for altitude, I might guess that they didn't
// correct parallax for the actual altitude of the summit, but rather just
// used the reference geoid or even something simpler. The USNO's notes
// don't say anything about parallax or refraction, but my previous encounters
// with USNO predictions show that they do have some refraction model close
// to my defaults (but theirs do a crazy step-function to zero as soon as the
// centerline of the sun or moon crosses the horizon.) Hmmm... let's
// calculate the magnitude of possible error due to parallax.
println["\nMeta-Calculation of possible parallax discrepancy:"]
// The maximum magnitude of the parallax error would occur if the sun were
// at the horizon, but since it's high, parallax is reduced.
parallax = parallaxAngleAlt[sundist, airlessAlt]
println["Total parallax angle is: " + format[parallax, "degrees", 5]]
// Running this, the nominal parallax angle is about 0.00126 degrees, which
// is small (but nonzero). The error due to not adjusting for geodetic
// elevation would be:
parallaxErrorFactor = summitHeight / earthradius
println["Parallax error factor is: " + format[parallaxErrorFactor, 1, 5]]
println["If they didn't correct for parallax at all, error is " + format[parallax earthradius, "feet", 1] + "."]
println["If they corrected for parallax for the geoid, but not for the mountain summit, error is " + format[parallax earthradius parallaxErrorFactor, "feet", 1] + "."]
// So again, multiplying this is small, at least for the sun. (It would be
// significant for the moon, which is 389 times closer!)
// We're still stuck by not knowing how much they rounded the alt/azimuth
// figures. Without knowing this, we might be off by as much as a degree
// (about 70 miles) in lat or long!
//
// Clarification has been requested. We'll see what they say.
roundAlt = round[altitude, 0.1 degrees]
roundAz = round[azimuth, 0.1 degrees]
println["Rounded altitude is: " + format[roundAlt, "degrees", 5]]
println["Rounded azimuth is : " + format[roundAz, "degrees", 5]]
// "Actual latitude (N) of Waypoint #1 is the Sun's altitude
// minus 19.23007 degrees."
W1lat = roundAlt - 19.23007 degrees
// "Actual longitude (W) of Waypoint #1 is the Sun's azimuth
// minus 72.63358 degrees."
W1long = roundAz - 72.63358 degrees
println[]
println["Waypoint 1 latitude is: " + format[W1lat, "degrees", 5]]
println["Waypoint 1 longitude is: " + format[W1long, "degrees", 5]]
// "WAYPOINT #2
// "At Waypoint #1 you will find a Rx pill bottle containing "adjusted" UTM
// coordinates for WP #2 and a specific factor ________________used in
// calculating the coordinates of Waypoint 2._______________ N and
// _____________________ E"
// Hmmm... by "unadjusted" do they mean what that usually means for UTM
// coordinates, that the easting is shifted by 500000 meters and is relative
// to the bounding meridians and not to the central meridian? Or do they
// mean that they're "normal" UTM coordinates and the "unadjusted" means that
// they just haven't added the numbers below? Must ask.
// "WP 1 is
// located in a small park named after a location where terrible events took
// place during WWII. These events were memorialized by a poem set to music. A
// number in the music title is used in calculating WP2 coordinates.
// The UTM values found at WP1 must be adjusted as follows to obtain the
// actual WP2 coordinates:
println[]
println["Waypoint 2"]
println["a. The mountain peak elevation, in feet, " + (summitHeight->feet) + ", multiplied by 0.97588."]
a2 = (summitHeight/feet) * 0.97588
println["result is $a2"]
// Do they round the Julian day? To the nearest day, or do they use
// all the digits? Are they including the 7 PM UT in the JD calculation?
// I guess it doesn't matter too much, because these are offsets to UTM
// coordinates and thus just have units of meters.
b2 = JD[surrenderDate] / days
println["\nAdd the UT Julian date of the Japanese surrender signing: "]
println[format[b2, 1, 5] + ", to result A."]
b3 = a2 + b2
println["b. result is " + format[b3, 1, 5]]