NavList:

   

Reply

Again with apologies for the rough and ready presentation.  Same
explanation as in the past; if I wait 'til it's all cleaned up,
it'll probably never get posted.

Goal here is to produce an Assumed Position-like location for use in
the latitude by GMT and lunar distance method.  We might call it the
First order position (FOP), since it's based on the observed lunar
distances without refraction corrections.  (It can't include
refraction since we need altitude to calculate refraction, and we need
location to calculate altitude, and location is what we're trying to
find.)

We start by writing the equation of a cone whose vertex is at the
earth's center and whose axis points to the North Pole.  The cone
angle is the observed lunar distance to the first star center to
center.  We then rotate the equation of the cone so the axis points at
the first star of interest.  We next translate the cone to the center
of the moon.  Call this equation a(x,y,z).  Repeat for the second lunar
distance-star.  Call this equation b(x,y,z).  Finally, since we want
our FOP to be on the surface of the earth, we write
c(x,y,z)=x^2+y^2+z^2-1.  The point we seek lies on the surface of both
cones and on the surface of the earth.  We now have three equations
with three unknowns.  Since we have  higher order terms in x, y, and z,
we have a number of solutions, most imaginary.  We do have two real
solutions however, and these are the two of interest.  As has been
pointed out, one is on the side of the earth facing the moon and the
other on the opposite side.  Taking the one facing the moon, we now
have the starting point AP=FOP, for the calculation using
interpolation/extrapolation from four corners, or using the
azimuth/intercept method.

For anyone interested, below is a FORTRAN program to find the cone
equations, and a MAXIMA script to solve the system of equations.  The
example is that given by Frank.  (With thanks again!)  The answers at
the end are:             71.2W         38.8W

Be warned that errors/bugs may certainly still remain.


    real rotm(3,3)
    pi=4.*atan(1.0)
    print*,pi
    kang=0
c true ld1
    theta=38.+(13.1+16.1)/60
    aw=0.
    np=1
    tol=1.e-7
    print*,'kang, theta, aw, np, tol'
    print*,kang, theta, aw, np, tol
      call aptcois (kang, theta, aw, np, tol, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, nerr)
    print*,'qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr'
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr
999 format(10f8.5,i2)
    n1=2
c star1 gha=3.035 dec=28.01167
    th1=(90.-28.01167)
    n2=3
    th2=-3.035
    n3=1
    th3=0.
    print*,'n1, th1, n2, th2, n3, th3, kang, tol'
    print'(i5,f12.8)',n1, th1, n2, th2, n3, th3, ku, tol
      call aptrots (n1, th1, n2, th2, n3, th3, kang, tol,
     &                    rotm, nerr)
    print*,'rotm,nerr'
    print'(3f10.6)',rotm
    print*,nerr
    inv=0
    ax=0.
    ay=0.
    az=0.
      call aptrois (rotm, inv, ax, ay, az, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr
    print*,'moon pos'
    r=58.852

c moon1 dist=58.852 earth radii gha=47.1 dec=27.64833
    ph=pi/180*(90-27.64833)
    th=pi/180*(-47.1)
    print*,'r,ph,th',r,ph,th
c x=rsin ph cos th  y=r sin phi sin th  z=rcos ph ph=90-dec th=360-gha
    dx=r*sin(ph)*cos(th)
    dy=r*sin(ph)*sin(th)
    dz=r*cos(ph)
c      subroutine apttris (nopt, ax, ay, az, qc, qx, qy, qz,
c     &                    qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
    nopt=1
    print*,'nopt= ',nopt
      call apttris (nopt, dx, dy, dz, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
    print*,'after apttris1, translate,  nopt= ',nopt
    print*,qc,qx,qy,qz,qxy
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr

    print*,'BEGIN STAR2'
c star2 gha 49.635  dec=16.525
c moon2 gha 46.38 dec=27.645
c pos= 69.6 dec=38.6
c true ld2 11.51992
    theta=11.51992/2.
    theta=11.51992
    theta=11.665
    theta=11.+(39.9-16.1)/60.
    print*,'kang, theta, aw, np, tol'
    print*,kang, theta, aw, np, tol
      call aptcois (kang, theta, aw, np, tol, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, nerr)
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr
    n1=2
    th1=(90.-16.525)
    n2=3
    th2=-49.635
    n3=1
    th3=0.
    print*,'n1, th1, n2, th2, n3, th3, kang, tol'
    print'(i5,f12.8)',n1, th1, n2, th2, n3, th3, ku, tol
      call aptrots (n1, th1, n2, th2, n3, th3, kang, tol,
     &                    rotm, nerr)
    print*,'rotm,nerr'
    print'(3f10.6)',rotm
    print*,nerr
      call aptrois (rotm, inv, ax, ay, az, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr
    ph=pi/180.*(90.-27.645)
    th=pi/180.*(-46.38)
    dx=r*sin(ph)*cos(th)
    dy=r*sin(ph)*sin(th)
    dz=r*cos(ph)
    print*,'dx,dy,dz',dx,dy,dz
    nopt=1
    print*,'nopt= ',nopt
      call apttris (nopt, dx, dy,dz, qc, qx, qy, qz,
     &                    qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
    print*,'after apttris, translate,  nopt= ',nopt
    print*,qc,qx,qy,qz,qxy
    print999,qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr

    end


C:\G77\bin>a.exe
  3.14159274
 kang, theta, aw, np, tol
 0  38.4866676  0. 1  1.00000001E-007
 qc, qx, qy, qz,qxy, qyz, qzx, qxx, qyy, qzz, nerr
 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.61270
0.61270-0.38730
 0
 n1, th1, n2, th2, n3, th3, kang, tol
    2 61.98833084
    3 -3.03500009
    1  0.00000000
    0  0.00000010
 rotm,nerr
  0.468993 -0.024866 -0.882852
  0.052946  0.998597  0.000000
  0.881614 -0.046743  0.469651
 0
 0.00000 0.00000 0.00000 0.00000 0.08242 0.04391-0.82810-0.16454
0.61052 0.39213
 0
 moon pos
 r,ph,th  58.8520012  1.08824193 -0.822050095
 nopt=  1
 after apttris1, translate,  nopt=  1
  15.5745239  37.4410362  42.5059509  9.64573669  0.0824193433
15.5745237.4410442.50595 9.64574 0.08242 0.04391-0.82810-0.16454
0.61052 0.39213
 0
 BEGIN STAR2
 kang, theta, aw, np, tol
 0  11.3966665  0. 1  1.00000001E-007
 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.96095
0.96095-0.03905
 0
 n1, th1, n2, th2, n3, th3, kang, tol
    2 73.47499847
    3-49.63499832
    1  0.00000000
    0  0.00000010
 rotm,nerr
  0.184215 -0.216720 -0.958696
  0.761934  0.647655  0.000000
  0.620904 -0.730463  0.284434
 0
 0.00000 0.00000 0.00000 0.00000 0.90709 0.41554-0.35321 0.57543
0.42738 0.88005
 0
 dx,dy,dz  35.9654045 -37.7410126  27.3068504
 nopt=  1
 after apttris, translate,  nopt=  1
  2.90124512  2.48841476 -11.711647 -19.6767006  0.90709424
 2.90125 2.48841**************** 0.90709 0.41554-0.35321 0.57543
0.42738 0.88005
 0




%i248) batch("/mnt/floppy/simpline3");

batching #p/mnt/auto/floppy/simpline3.mac
(%i249) a(x, y, z) := 0.39213 z z + 0.61052 y y - 0.16454 x x - 0.8281
z x
       + 0.04391 y z + 0.082419 x y + 9.6457 z + 42.5059 y + 37.441 x +
15.5745
(%o249) a(x, y, z) := 0.39213 z z + 0.61052 y y - 0.16454 x x - 0.8281
z x
       + 0.04391 y z + 0.082419 x y + 9.6457 z + 42.5059 y + 37.441 x +
15.5745
(%i250) b(x, y, z) := 0.88005 z z + 0.42738 y y + 0.57543 x x - 0.35321
z x
      + 0.41554 y z + 0.90709 x y - 19.6767 z - 11.7116 y + 2.48841 x +
2.90125
(%o250) b(x, y, z) := 0.88005 z z + 0.42738 y y + 0.57543 x x - 0.35321
z x
      + 0.41554 y z + 0.90709 x y - 19.6767 z - 11.7116 y + 2.48841 x +
2.90125
(%i251)               c(x, y, z) := - 1 + z z + y y + x x
(%o251)               c(x, y, z) := - 1 + z z + y y + x x
(%i252)        xx : solve([a(x, y, z), b(x, y, z), c(x, y, z)])

`rat' replaced 15.5745 by 31149//2000 = 15.5745

`rat' replaced 37.441 by 37441//1000 = 37.441

`rat' replaced -0.16454 by -2579//15674 = -0.164540002552

`rat' replaced 42.5059 by 28819//678 = 42.50589970501475

`rat' replaced 0.082419 by 1277//15494 = 0.08241900090358

`rat' replaced 0.61052 by 4074//6673 = 0.61052000599431

`rat' replaced 9.6457 by 27143//2814 = 9.645700071073206

`rat' replaced -0.8281 by -6725//8121 = -0.8280999876862

`rat' replaced 0.04391 by 1813//41289 = 0.0439100002422

`rat' replaced 0.39213 by 2003//5108 = 0.39212999216915

`rat' replaced 2.90125 by 2321//800 = 2.90125

`rat' replaced 2.48841 by 13741//5522 = 2.488409996378124

`rat' replaced 0.57543 by 10806//18779 = 0.57543000159753

`rat' replaced -11.7116 by -25543//2181 = -11.7116001834021

`rat' replaced 0.90709 by 4657//5134 = 0.90708998831321

`rat' replaced 0.42738 by 4817//11271 = 0.42738000177447

`rat' replaced -19.6767 by -17650//897 = -19.6767001114827

`rat' replaced -0.35321 by -6415//18162 = -0.3532099988988

`rat' replaced 0.41554 by 5225//12574 = 0.41554000318117

`rat' replaced 0.88005 by 17601//20000 = 0.88005
(%o252) [[z = 186.9795605059265 - 71.03037545338685 %i,
y = 156.6552474191666 %i + 78.65861212645633,
x = 107.8904598096931 %i + 8.888126366487127],
[z = 71.03037545338685 %i + 186.9795605059265,
y = 78.65861212645621 - 156.6552474191664 %i,
x = 8.888126366487402 - 107.8904598096932 %i],
[z = - 62.46498747201671 %i - 14.20858464799212,
y = 20.0682027972187 - 60.980812743406 %i,
x = 4.009354556116049 %i + 83.86293856926758],
[z = 62.46498747201671 %i - 14.20858464799212,
y = 60.98081274340601 %i + 20.06820279721871,
x = 83.86293856926757 - 4.009354556116057 %i],
[z = 16.55536532251499 - 3.291714645530333 %i,
y = 11.28740475395483 %i - 2.988682285874108,
x = - 13.62692118616466 %i - 6.474683747745093],
[z = 3.291714645530333 %i + 16.55536532251499,
y = - 11.28740475395483 %i - 2.988682285874106,
x = 13.62692118616466 %i - 6.474683747745092],
[z = 0.62787950383934, y = - 0.73687926570958, x = 0.25055201974282],
[z = - 0.22476558190844, y = 0.45003531073446, x = - 0.86426198386331]]
(%i253)                        yy : seventh(xx)
(%o253) [z = 0.62787950383934, y = - 0.73687926570958, x =
0.25055201974282]
(%i254)                       zf : rhs(first(yy))
(%o254)                        0.62787950383934
(%i255)                      yf : rhs(second(yy))
(%o255)                       - 0.73687926570958
(%i256)                       xf : rhs(third(yy))
(%o256)                        0.25055201974282
                                        yf
                               180 atan(--)
                                        xf
(%i257)                     ev(------------, numer)
                                   %pi
(%o257)                       - 71.22105587662567
                                 180 acos(zf)
(%i258)                  ev(90 - ------------, numer)
                                     %pi
(%o258)                        38.8938488718472
(%i259)                         yy : eighth(xx)
(%o259) [z = - 0.22476558190844, y = 0.45003531073446, x = -
0.86426198386331]
(%i260)                       zf : rhs(first(yy))
(%o260)                       - 0.22476558190844
(%i261)                      yf : rhs(second(yy))
(%o261)                        0.45003531073446
(%i262)                       xf : rhs(third(yy))
(%o262)                       - 0.86426198386331
                                        yf
                               180 atan(--)
                                        xf
(%i263)                     ev(------------, numer)
                                   %pi
(%o263)                       - 27.50672813187197
                                 180 acos(zf)
(%i264)                  ev(90 - ------------, numer)
                                     %pi
(%o264)                       - 12.98909392363139
(%i265)


--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to NavList@navlist.net
To , send email to NavList-@navlist.net
-~----------~----~----~----~------~----~------~--~---

   

Reply