forensicSeismology.frink

View or download 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)]


View or download 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 17596 days, 22 hours, 5 minutes ago.