 Check out the bookstore at IRBS.com |
|
|
|
|||||
|
||||||
Subject: Barrie Hudson Challenge
From: Dan Allen (danallen46@XXX.XXX)
Date: Sat Oct 19 2002 - 01:20:08 EDT
On Friday, October 18, 2002, at 05:12 PM, bhudson wrote:
> Hi,
> Now that we can handle the Cos formula and GC sailing here is a poser!
> A ship sails from a position off Valpraiso Lat 33°01'S, Long 72°10'W
> to a
> position off San Francisco Lat 37°50'N, Long 123°14'W
> a) Find the initial course and GC distance.
> b) Find the position of the ship as it crosses the equator and the
> course it
> crosses the equator.
> Barrie Hudson
I got the following results:
DistNMI Latitude Longitude Course
---------------------------------------
0 -33.01667 72.16667 321.94335 Valpariso
1000 -19.42858 82.97198 326.76300
2000 -5.30354 92.05548 328.72709
2371 0.00000 95.26867 328.87595 Equator
2559 2.66837 96.88125 328.83839 Midpoint
3000 8.94774 100.72401 328.44879
4000 22.98133 110.10570 325.84358
5000 36.33635 121.63659 320.08411
5118 37.83333 123.23333 26.75367 San Francisco
---------------------------------------
which were generated by the following Awk program.
# Usage: awk -f hudson.awk
BEGIN { # all arguments and results are in decimal degrees
CONVFMT = OFMT = "%10.5f"
PI = 4*atan2(1,1)
DTOR = PI/180
RTOD = 180/PI
# Near Valpariso, Chile
lat1 = -(33 + 1/60)
lon1 = 72+10/60
# Near SF, CA
lat2 = 37+50/60
lon2 = 123+14/60
# Let's go!
d = GCDistance(lat1,lon1,lat2,lon2)
c = GCCourse(lat1,lon1,lat2,lon2)
print
print "DistNMI Latitude Longitude Course"
print "---------------------------------------"
for (i = 0; i <= d; i += 1000) {
s = GCPoint(lat1,lon1,c,i)
printf("%6d %s%s\n",i,s,i == 0 ? " Valpariso" : "")
if (i == 2000) {
EquatorCrossing(lat1,lon1,lat2,lon2)
printf("%6.0f %s Midpoint\n",d/2,GCPoint(lat1,lon1,c,d/2))
}
}
s = GCPoint(lat1,lon1,c,d)
printf("%6.0f %s San Francisco\n",d,s)
print "---------------------------------------"
}
function Abs(x) { return x < 0 ? -x : x }
function Floor(x) { return x < 0 ? int(x) - 1 : int(x) }
function Round(x) { return Floor(x+0.5) } # Everyday round > 0
function Mod(x,y) { return x - y * Floor(x/y) }
function Sin(x) { return sin(x*DTOR) }
function Cos(x) { return cos(x*DTOR) }
function Tan(x) { return Sin(x)/Cos(x) }
function ASin(x) { return atan2(x,sqrt(1 - x * x))*RTOD }
function ACos(x) { return atan2(sqrt(1 - x * x),x)*RTOD }
function ATan2(y,x){ return atan2(y,x)*RTOD }
function GCDistance(lat1,lon1,lat2,lon2) {
return 60*2*ASin(sqrt((Sin((lat1-lat2)/2))^2 +
Cos(lat1)*Cos(lat2)*(Sin((lon1-lon2)/2))^2))
}
function GCCourse(lat1,lon1,lat2,lon2) {
return
Mod(ATan2(Sin(lon1-lon2)*Cos(lat2),Cos(lat1)*Sin(lat2)-
Sin(lat1)*Cos(lat2)*Cos(lon1-lon2)),360)
}
function GCPoint(lat1,lon1,c,distance, d,dlon) { # lat,lon need to be
globals
d = distance/60
lat = ASin(Sin(lat1)*Cos(d) + Cos(lat1)*Sin(d)*Cos(c))
dlon = ATan2(Sin(c)*Sin(d)*Cos(lat1),Cos(d)-Sin(lat1)*Sin(lat))
lon = Mod(lon1 - dlon + 180,360) - 180
return lat " " lon " " GCCourse(lat,lon,lat2,lon2) # lat2, lon2 need
to be globals
}
function GCLatitude(lat1,lon1,lat2,lon2,lon3) {
return
ATan2(Sin(lat1)*Cos(lat2)*Sin(lon3-lon2)-Sin(lat2)*Cos(lat1)*Sin(lon3-
lon1),
Cos(lat1)*Cos(lat2)*Sin(lon1-lon2))
}
function GCLongitude(lat1,lon1,lat2,lon2,lat3,
dlon,lon,l12,lon3_1,lon3_2,A,B,C) {
l12 = lon1 - lon2
A = Sin(lat1)*Cos(lat2)*Cos(lat3)*Sin(l12)
B = Sin(lat1)*Cos(lat2)*Cos(lat3)*Cos(l12) -
Cos(lat1)*Sin(lat2)*Cos(lat3)
C = Cos(lat1)*Cos(lat2)*Sin(lat3)*Sin(l12)
lon = ATan2(B,A)
if (Abs(sqrt(A^2+B^2) - C) < 0.0000001)
dlon = 0
else
dlon = ACos(C/sqrt(A^2+B^2))
lon3_1 = Mod(lon1+dlon+lon+180,360) - 180
lon3_2 = Mod(lon1-dlon+lon+180,360) - 180
if (lon3_1 > lon1 && lon3_1 < lon2)
return lon3_1
else
return lon3_2
}
function EquatorCrossing(lat1,lon1,lat2,lon2, d,lonV) {
lonV = GCLongitude(lat1,lon1,lat2,lon2,0)
d = GCDistance(lat1,lon1,0,lonV)
printf("%6d ",d)
print GCPoint(lat1,lon1,GCCourse(lat1,lon1,0,lonV),d), " Equator"
}
|