[futurebasic] Re: [FB] Great Circle Cal

Message: < previous - next > : Reply : Subscribe : Cleanse
Home   : January 1999 : Group Archive : Group : All Groups

From: Al Boyd <aboyd@...>
Date: Fri, 8 Jan 99 15:24:04 -0600
>The reason this appears here is that
>the formula that calculates "a" is never _supposed_ to exceed 1, but
>depending on the rounding errors of the particular computer, you could
>possibly get a number slightly greater than 1.  the "min" is just a way
>to limit the number to no greater than 1.
>
>- Rick

Thanks Rick,

Actually the information I had obtained about great circle calculations 
explained this, but I was so confused over the whole thing... as the 
saying goes "Couldn't see the trees because of the forest" (or something 
like that). Anyway, here's the formula done in FB, if anyone sees a 
problem, please let me know. I don't know how accurate the routine is, 
but there are a couple of websites I know of where you can enter 
coordinates and get distances, and the routine is giving me same or near 
same results. One site gives me small differences, but it probably has to 
do with what value is used for the radius of the earth.

Al

'---------------------

WINDOW 1

'================ FUNCTIONS ===================

'===============================================================
'Information for the great circle routine was obtained from:
'http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1
'The great circle formula in the information was stated to have
'come from:
'Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
'Sky and Telescope, vol. 68, no.2, 1984, p. 159)
'===============================================================

LOCAL FN arcsin#(value#)
  'FutureBASIC does not have a ARCSIN function, so the
  'below formula is used. Obtained from the Math
  'Appendix in the FB Reference Manual
END FN = ATN(value#/SQR(1-value# * value#))

LOCAL FN greatCircleDis#(R#,lat1#,lat2#,dlat#,dlon#)
  a# = (SIN(dlat#/2))^2 + COS(lat1#) * COS(lat2#) * (SIN(dlon#/2))^2
  aa# = SQR(a#)
  'If the SQR(a#) is greater then 1 the the inverse of sine (or arcsin) 
wouldn't work
  IF aa# > 1 THEN aa# = 1
  b# = FN arcsin#(aa#)
  c# = 2 * b#
  d# = R# * c#
END FN = d#

'================= ROUTINE ====================
'R# is the radius of the earth. The earth is not a
'perfect sphere, but these calculations are accurate
'enough for my needs.
'Statute Miles  = 3963.1
'Nautical Miles = 3443.9
'Kilometers     = 6378.0
'
'If coordinates are in degrees/minutes/seconds you must convert to
'decimal. Then for the calculations to work properly 
'the degrees must be converted to radians by
'multipling degrees by pi/180 (or .017453293)
'If coordinates cross prime meridian and/or equator, one
'degree must be positive, and the other negitive.

'You'll need Latitude/Longitude coordinates from two different locations.
'Example coordinates: 23.7S 18.3E to 14.4N 9.5W

lati1# = -23.7 * .017453293
long1# = 18.3  * .017453293

lati2# = 14.4 * .017453293
long2# = -9.5 * .017453293

dlong# = long2# - long1#
dlati# = lati2# - lati1#

'''''''''''''
Rad# = 3963.1
dis# = FN greatCircleDis#(Rad#,lati1#,lati2#,dlati#,dlong#)
TEXT ,,_boldBit%
PRINT
PRINT dis#;" Statute Miles"
'''''''''''''
Rad# = 3443.9
dis# = FN greatCircleDis#(Rad#,lati1#,lati2#,dlati#,dlong#)
PRINT
PRINT dis#;" Nautical Miles"
'''''''''''''
Rad# = 6378.0
dis# = FN greatCircleDis#(Rad#,lati1#,lati2#,dlati#,dlong#)
PRINT
PRINT dis#;" Kilometers"
'''''''''''''
STOP