| |
DK2ZA > TECHNIK 05.08.06 21:09l 294 Lines 7687 Bytes #999 (99) @ WW
BID : J3JIXIDB0FOR
Read: DL1LCA DL8FBH GUEST OE7FMI
Subj: re: Help Formulas PA3BNX
Path: DB0FHN<DB0FOR
Sent: 060805/1949z @:DB0FOR.#BAY.DEU.EU [Forchheim/JN59MO] DP6.00 $:J3JIXIDB0FO
From: DK2ZA @ DB0FOR.#BAY.DEU.EU (Helmut)
To: TECHNIK @ WW
Hello Lodewijk,
I read your "search for help with formulas" and found the
problem very interesting. I could not find a solution in the
internet, so I tried it myself.
The following programm seems to be ok except for the case when one
Point is at the North- or Southpole.
It always gives two solutions since the problem really has two solutions.
Here are two examples from the literature:
P1: longitude = -74.1 Deg, latitude = 40.7 Deg, bearing = 61.333333 Deg
P2: longitude = -9.2 Deg, latitude = 38.7 Deg, bearing = 311.81666 Deg
Result: longitude = -27.05 Deg, latitude = 48.08 Deg
(The program gives 152.95,-48.08 and 332.95,48.08 additionally)
P1: longitude = -80, latitude = 32.7, bearing = 112.78333
P2: longitude = -8.65, latitude = 41.166666, bearing = 250.86666
Result: longitude = -51.55, latitude = 18.09
(The program gives 128.45,-18.09 and 308.45,18.09 additionally)
Contact me via Packet or DK2ZA at DARC.DE if you have any questions
or if you have already found a better solution
73 de Helmut, DK2ZA
PRINT " This program calculates the coordinates of a point X"
PRINT " on the surface of the (assumed spherical) earth."
PRINT
PRINT " Given are the coordinates of two points P1 and P2"
PRINT " and the bearings from these points to point X."
PRINT
PRINT
PRINT " Input coordinates of P1 and P2:"
PRINT
PRINT " West of Greenwich: negative longitudes 0 ... -180"
PRINT " East of Greenwich: positive longitudes 0 ... 360"
PRINT
PRINT " North of the Equator: positive latitudes 0 ... 90"
PRINT " South of the Equator: negative latitudes 0 ... -90"
PRINT
INPUT " longitude of P1 (-180 ... 360) = ";lon1
lon1=PI*lon1/180
INPUT " latitude of P1 (-90 ... 90) = ";lat1
lat1=PI*lat1/180
PRINT
INPUT " longitude of P2 (-180 ... 360) = ";lon2
lon2=PI*lon2/180
INPUT " latitude of P2 (-90 ... 90) = ";lat2
lat2=PI*lat2/180
PRINT
PRINT
PRINT " Input bearings from P1 and P2 to Point X"
PRINT
PRINT " North = 0 or 360, East = 90, South = 180 or -180, West = 270 or -90"
PRINT
INPUT " bearing at P1 (-180 ... 360) = ";bear1
PRINT
INPUT " bearing at P2 (-180 ... 360) = ";bear2
'
' We now look at the Great Circle through P1 whose direction in P1
' is given by bear1 and at the Great Circle through P2 whose direction
' in P2 is bear2.
'
' These two Great Circles intersect each other at X1 and X2 which both
' can be the Point we are looking for. X1 and X2 lie on a straight line
' through the center of the earth.
'
' A bearing of 120 Deg gives the same Great Circle as a bearing of
' of -60 Deg so only bearings between -90 Deg and +90 Deg are necessary:
'
' TRUNC may be called INT in other languages
' TRUNC(3.2) = 3 TRUNC(-3.2) = -3 TRUNC(0.7) = 0 TRUNC(-0.7) = 0
'
bear1=bear1+180-180*TRUNC(bear1/180+1.5)
bear1=PI*bear1/180
'
bear2=bear2+180-180*TRUNC(bear2/180+1.5)
bear2=PI*bear2/180
'
'
' The Great Circle through P1 with direction bear1 in P1 intersects
' the Equator in S1 under the angle alpha1 (0 <= alpha1 <= 90):
'
alpha1=ACOS(SIN(ABS(bear1))*COS(lat1))
'
' c1 is the angle between straight lines coming from the center of
' the earth and going through S1 and through a point on the equator
' which has the same longitude as P1.
'
' c1 = ASIN(SIN(lat1)*SIN(ABS(bear1))/SIN(alpha1))
'
' First the argument of the ASIN (arcus sinus) function is calculated:
'
c1=SIN(lat1)*SIN(ABS(bear1))/SIN(alpha1)
'
' The ASIN function only likes arguments between -1 and +1 but
' c1 may be outside this interval due to rounding errors, so
'
IF c1<-1 THEN
c1=-1
ENDIF
IF c1>1 THEN
c1=1
ENDIF
'
' Now the ASIN function can be applied safely:
'
c1=ASIN(c1)
'
' Since alpha1 has to be between 0 and 90 Degrees we may have to take
' the other intersection of the Great Circle with the equator:
'
IF bear1>=0 THEN
lons1=lon1-c1
ELSE
lons1=lon1+c1+PI
ENDIF
'
'
' The Great Circle through P2 with direction bear2 in P2 intersects
' the Equator in S2 under the angle alpha2:
'
alpha2=ACOS(SIN(ABS(bear2))*COS(lat2))
'
c2=SIN(lat2)*SIN(ABS(bear2))/SIN(alpha2)
'
IF c2<-1 THEN
c2=-1
ENDIF
IF c2>1 THEN
c2=1
ENDIF
'
c2=ASIN(c2)
'
IF bear2>=0 THEN
lons2=lon2-c2
ELSE
lons2=lon2+c2+PI
ENDIF
'
'
' Now we have two Great Circles given by their intersection points
' with the equator (at lons1 and lons2) and the angles at which
' they intersect the equator (alpha1 and alpha2).
'
' These two Great Circles intersect each other two times.
' To find out where we use a formula that calculates the latitude
' of a point on a Great Circle given the longitude of that point.
'
' For the Great Circle through P1 the formula looks like this:
' (Given: longitude Result: latitude1)
'
' gamma1=ACOS(SIN(alpha1)*COS(longitude-lons1))
' latitude1=ASIN(SIN(longitude-lons1)*SIN(alpha1)/SIN(gamma1))
'
' The two circles intersect where latitude1 = latitude2 for the
' same longitude.
'
' Since the two intersections are 180 Deg apart there is exactly
' one intersection in the intervall 0 <= longitude < 180 Deg.
'
' So we are looking for a value of longitude (0 ... 180) for which
' latitude1 - latitude2 becomes 0.
'
' The simple method of repeatedly halving the intervall from
' longlo = 0 to longhi = 180 is used.
'
longlo=0
longhi=PI
'
REPEAT
'
longmid=(longlo+longhi)/2
'
'
longitude=longlo
'
' Calculate the latitude of a point on Great Circle 1 (latitude1)
' whose longitude is longlo
'
gamma1=ACOS(SIN(alpha1)*COS(longitude-lons1))
latitude1=SIN(longitude-lons1)*SIN(alpha1)/SIN(gamma1)
IF latitude1<-1 THEN
latitude1=-1
ENDIF
IF latitude1>1 THEN
latitude1=1
ENDIF
latitude1=ASIN(latitude1)
'
' Calculate the latitude of a point on Great Circle 2 (latitude2)
' whose longitude is longlo
'
gamma2=ACOS(SIN(alpha2)*COS(longitude-lons2))
latitude2=SIN(longitude-lons2)*SIN(alpha2)/SIN(gamma2)
IF latitude2<-1 THEN
latitude2=-1
ENDIF
IF latitude2>1 THEN
latitude2=1
ENDIF
latitude2=ASIN(latitude2)
'
' The difference of latitudes at longitude longlo is difflo
'
difflo=latitude1-latitude2
'
'
' Same procedure for longmid:
'
longitude=longmid
'
gamma1=ACOS(SIN(alpha1)*COS(longitude-lons1))
latitude1=SIN(longitude-lons1)*SIN(alpha1)/SIN(gamma1)
IF latitude1<-1 THEN
latitude1=-1
ENDIF
IF latitude1>1 THEN
latitude1=1
ENDIF
latitude1=ASIN(latitude1)
'
gamma2=ACOS(SIN(alpha2)*COS(longitude-lons2))
latitude2=SIN(longitude-lons2)*SIN(alpha2)/SIN(gamma2)
IF latitude2<-1 THEN
latitude2=-1
ENDIF
IF latitude2>1 THEN
latitude2=1
ENDIF
latitude2=ASIN(latitude2)
'
' The difference of latitudes at longitude longmid is diffmid
'
diffmid=latitude1-latitude2
'
' Now we choose that half which contains the zero of latitude1 - latitude2
' For this half the value of latitude1 - latitude2 has different signs
' at the borders of that half
' That means that the product of these numbers is negative
'
IF difflo*diffmid<0 THEN
longhi=longmid
ELSE
longlo=longmid
ENDIF
'
' The interval that contains the longitude of point X is now
' half as long as before
'
' The procedure is repeated until the difference between latitudes
' is small enough
'
UNTIL ABS(diffmid)<1.0E-10
'
'
longitude=180*longmid/PI
latitude=180*latitude1/PI
'
PRINT
PRINT
PRINT "Coordinates of point X:"
PRINT
PRINT "longitude: ";longitude,
PRINT "latitude: ";latitude
'
PRINT
PRINT "Second solution:"
PRINT
PRINT "longitude: ";longitude+180,
PRINT "latitude: ";-latitude
'
PRINT
PRINT "which is the same as:"
PRINT
PRINT "longitude: ";longitude-180,
PRINT "latitude: ";-latitude
Read previous mail | Read next mail
| |
