Download or view forensicSeismology.frink in plain text format
// Solution to the "Forensic Seismology" geocache:
//
// http://www.geocaching.com/seek/cache_details.aspx?guid=68083e90-f516-46de-a2c9-55709818a7ed
// Import the sun.frink file which contains great circle calculations
// and definitions for many of the functions below.
use navigation.frink
// Coordinates of seismic stations
lata = DMS[39,30.840] North
longa = DMS[104,44.109] West
latb = DMS[39, 31.710] North
longb = DMS[104, 47.901] West
latc = DMS[39, 27.302] North
longc = DMS[104, 46.119] West
// Find the distances between seismic stations
[dab, bab] = earthDistanceAndBearing[lata, longa, latb, longb]
[dbc, bbc] = earthDistanceAndBearing[latb, longb, latc, longc]
dca = earthDistance[latc, longc, lata, longa]
println["dab: $dab"]
println["dbc: $dbc"]
println["dca: $dca"]
// Find bearings between seismic stations.
println["Bearing from a to b: " + (bab->DM)]
bac = earthBearing[lata, longa, latc, longc]
println["Bearing from a to c: " + (bac->DM)]
println["Bearing from b to c: " + (bbc->DM)]
// Times of arrival
ta1 = .1673 s
tb1 = .2890 s
tc1 = .5507 s
ta2 = .4240 s
tb2 = .6065 s
tc2 = .9990 s
// Velocities of waves
v1 = 6.0 km/s
v2 = 4.0 km/s
// This function was one I derived for finding the distance between
// the meteorite and the seismic station, given times of arrival for
// the P-waves and S-waves and their velocities.
d[t1, t2, v1=v1, v2=v2] := -((t1-t2) v1 v2)/(v1-v2)
// Find distances to seismic stations from the meteorite.
// I call the meteorite "point D"
dad = d[ta1, ta2]
dbd = d[tb1, tb2]
dcd = d[tc1, tc2]
println["dad: $dad"]
println["dbd: $dbd"]
println["dcd: $dcd"]
// Cosine rule for triangles... find the angle DAB between the meteor,
// point A and point B.
adab = arccos[(dad^2 + dab^2 - dbd^2)/(2 dad dab)]
println["Angle DAB: " + (adab -> DM)]
// You can either add or subtract this angle to the angle bab (the bearing
// from point a to point b.) I tried both ways, and only subtracting gave a
// self-consistent solution when all points were considered.
bad = bab - adab
println["Bearing from A to D: " + (bad -> DM)]
// Use calculation to find the resultant lat/long of the meteorite given
// the initial point (that of seismic station a), the distance, and the
// bearing of the meteorite.
[latd, longd] = resultantLatLong[lata, longa, dad, bad]
println["D is at lat: " + (latd->DM)]
println["D is at long: " + (longd->DM)]
// Reverse calculations. These distances and bearings should match the
// ones calculated above.
println["\nVerifying:"]
[dad, bad] = earthDistanceAndBearing[lata, longa, latd, longd]
println["Bearing from A to D: " + (bad -> DM)]
println["Distance from A to D: " + dad]
[dbd, bbd] = earthDistanceAndBearing[latb, longb, latd, longd]
println["Bearing from B to D: " + (bbd -> DM)]
println["Distance from B to D: " + dbd]
[dcd, bcd] = earthDistanceAndBearing[latc, longc, latd, longd]
println["Bearing from C to D: " + (bcd -> DM)]
println["Distance from C to D: " + dcd]
// Karen's solutions (assuming UTM zone 13, easting 519926, northing 4372660)
latk = DMS[39,30.199999] North
longk = DMS[104,46.093333] West
[diffDist, diffBearing] = earthDistanceAndBearing[latd, longd, latk, longk]
println["Difference in solutions is " + (diffDist->"ft")]
println["Bearing from Alan's solutions to Karen's: " + (diffBearing->DM)]
Download or view forensicSeismology.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19966 days, 20 hours, 47 minutes ago.