NavList:

   

Reply

I agree with Peter Fogg that there's something amiss with Dan Allen's
program, in that it finds itself making a sudden dog-leg right at the end
ov the voyage.

I also agree with Peter's result, which he has achieved by some process of
successive approximation.

However, it's not necessary to do it that way. It can be done in a single
shot, remembering that a great-circle crosses the equator in two places
exactly 90� away, in longitude, from the longitude of each vertex. The
vertex is the point on a great circle which is nearest to the North (or
South) pole.

As usual, Bowditch provides the information about how to do it. Once the
initial course C has been calculated (which Dan and Peter have done without
problems) passing through the lat and long of Valparaiso, that defines the
great circle which will arrive at San Francisco. There's no further need to
consider San Francisco's coordinates.

Mine is the 2-volume (1981) edition of Bowditch. In vol 1, the relevant
diagrams are fig. 1016, a and b, on page 605.

The formulae are on page 604, in which we can stir together equations 25
and 27, to give the cos of the difference in long from Valparaiso to the
equator crossing as cos course / sqrt(1- (cos lat * sin course)^2).

lat is the latitude of Valparaiso, and course is the initial great-circle
course from Valparaiso to San Francisco.

Once the long of the Equator crossing has been established (I make it
W95.2687�, the same as Peter does ) then the result can be easily checked
by working out the great-circle course from that point to Valparaiso and to
San Francisco, which should differ by exactly 180� (they do).

George Huxtable.

========================

Peter Fogg said-

>a) The initial course is 321.9d and distance 5117.9 nm.
>b) At the equator I get a longitude of W95d 16' and course 327.8d
>This is very close to Dan's 95.26867, which is W95d 16m 7.2s, or W95d
>16.12'
>I have no idea why the course is different. To find the equator I had to
>slice the segments into smaller and smaller lengths, they are all 327.8d
>around those waters. Dan's courses seem to make sense until the last one
>which may be just a quirk of his program.
>
>Dan Allen wrote:
>
>> 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"
>> }

------------------------------

george@huxtable.u-net.com
George Huxtable, 1 Sandy Lane, Southmoor, Abingdon, Oxon OX13 5HX, UK.
Tel. 01865 820222 or (int.) +44 1865 820222.
------------------------------

   

Reply