|
Hi all, Back on trig functions again – trying to write a buffer function
that doesn’t create hundreds of points for each node. Overall, what I am trying to get is a calculator that will calculate
the endpoint coordinates from a calculated bearing and distance from an origin
point in lat/ long. I am using the function in cogoline.mb, or rather a
function that originated from the same source (same comments). To get a
line parallel from a line, I am first calculating a set of polygons offset to
the left and right from each segment. The trouble is when I calculate a
point a distance at 90 degrees from the calculated bearing of the segment, it
isn’t quite 90 degrees. My segment bearing calculation is not a spherical
calculation and doesn’t work correctly (I degree in latitude does not
equal 1 degree in longitude). TO solve it, I think that I need to put in
a bearing calculation that takes this into account. I have put in a great circle bearing calculation function below that
was adapted from a published source. The issue that I have is when I run
it on a set of lines that smallish, I get a domain error on the sqr function on
the line below. It is trying to find the square root of a negative number. It is adapted from a program in QBASIC, which does not have an ACos
function. I get the same error as when I run the program using ACos to
get the azimuth from the angle. I have actually used this function before
to calculate bearings between airports and it worked OK. I suspect rounding
errors. The values it fails on is when the CosAzimuth = -1.00059 or 1.00059 CosAzimuth = Sin(sglLat2 * Deg2Rad) - Cos(GreatCircleAngle) *
Sin(sgllat1 * Deg2Rad) CosAzimuth = CosAzimuth / Sin(GreatCircleAngle) / Then FAILS HERE: Azimuth = Atn(Sqr(1 - CosAzimuth * CosAzimuth) / CosAzimuth) Has anyone got a calculator that gives a bearing that accounts for the
difference in distance for units of latitude? I actually want to use it
in a VB program as well, so would prefer to use algebraic calculations than
Mapbasic if possible. I thought that a simple If Abs(CosAzimuth) > 1 Then make it 1 or -1
would work, but it seemed to spike it off wildly. Perhaps I did something
wrong. On further reading, the great circle bearing needs to be adjusted all
the way along its trajectory to get to the right destination, so perhaps this
is not the way to go at all. What function do surveyors use to do
this? I have a line and I want to calculate a parallel line 1.5 m to the
left? Any suggestions? R Function fGC_Bearing(ByVal sgllat1 As Float,
ByVal sglLong1 As Float,
ByVal sglLat2 As Float,
ByVal sglLong2 As Float) As Float 'This program is used to calculate the great circle bearings to and ' from any two sites on the Earth and the distances between the ' sites. Both the short and long paths along the great circle ' encompassing the two sites are considered. Standard spherical ' trignometrical formulae are used to perform these calculations. 'DIM LON(2) As Float
'array to hold the longitudes of the two sites 'DIM LAT(2) As Float
'array to hold the latitude of the two sites Dim TmpVal As Float Dim pi, Deg2Rad, RadiusEarth As Float Dim DiffLong As Float Dim CosCGA, GreatCircleAngle As Float Dim CosAzimuth, Azimuth As Float Dim Bearing As Float Dim Dist_Km As Float OnError GoTo ErrorHandler If sgllat1 = sglLat2 And sglLong1 = sglLong2 Then fGC_Bearing = 0 Exit Function End If pi = 3.14159 'the ratio of
circular circumference to diameter Deg2Rad = pi /
180 'constant to convert
degrees into radians RadiusEarth = 6372 'average radius of
Earth in kilometres 'Lat(1) = sgllat1 'sglLat2= sglLat2 'sglLong1 = sglLong1 'sglLong2 = sglLong2 DiffLong = sglLong2 - sglLong1 'difference in longitude between the two
sites 'now we calculate the COSINE of the Great Circle Angle
(GreatCircleAngle) between ' the two sites - this takes the next two lines - using spherical trig ' [ note that CosCGA is COS(GreatCircleAngle) ] CosCGA = Cos(sgllat1 * Deg2Rad) * Cos(sglLat2 * Deg2Rad) * Cos(DiffLong
* Deg2Rad) CosCGA = CosCGA + Sin(sgllat1 * Deg2Rad) * Sin(sglLat2 * Deg2Rad) 'we now take the inverse cosine of CosCGA to give GreatCircleAngle ' because QBASIC does not have an inverse cosine function we must ' use the inverse tangent function along with pythagorus ' this gives us the GreatCircleAngle in radians GreatCircleAngle = Atn(Sqr(1 - CosCGA * CosCGA) / CosCGA) 'we must now ensure that the GreatCircleAngle lies between 0 and 2*pi
radians If GreatCircleAngle <= 0 Then GreatCircleAngle = GreatCircleAngle + pi End If 'we now calculate the distance between the sites along the great circle ' this distance will be in the same units as RadiusEarth (ie
kilometres) Dist_Km = RadiusEarth * GreatCircleAngle 'we now calculate the COSINE of the bearing or azimuth ' this will take the next two lines - using spherical trigonometry ' [ note that CosAzimuth is COS(Azimuth) ] CosAzimuth = Sin(sglLat2 * Deg2Rad) - Cos(GreatCircleAngle) *
Sin(sgllat1 * Deg2Rad) CosAzimuth = CosAzimuth / Sin(GreatCircleAngle) / 'we now calculate the inverse cosine of CosAzimuth to give Azimuth the
azimuth ' of second site from the first Print "CosAzimuth " + CosAzimuth Azimuth = Atn(Sqr(1 - CosAzimuth * CosAzimuth) / CosAzimuth) 'and again we want the azimuth or bearing to lie from 0 to 2* pi (360
deg) If Azimuth < 0 Then Azimuth = pi + Azimuth End If 'this not only requires checking Azimuth itself but also the longitude ' difference If Sin(DiffLong * Deg2Rad) < 0 Then Azimuth = 2 * pi - Azimuth End If 'we now print bearing, after changing to degrees and rounding Bearing = Azimuth / Deg2Rad fGC_Bearing = Bearing Print "GC ANGLE" + GreatCircleAngle/ Deg2Rad + "
Bearing: " + Bearing 'Next line is used for debugging 'PRINT "Airport: " + AirportName + " Bearing is "+
Bearing+ " deg Distance is "+ Dist_Km+ " km" Exit Function ErrorHandler: Call
DISPLAYMESSAGE("ERROR","fGC_Bearing") fGC_Bearing = -999 End Function ------------------------------------------- Robert
Crossley Agtrix P/L
Far P: 61 2 6680 1309 F: 61 2 6680 5214 Cooparoo 4151 P: 61 7 3843 3363 |
_______________________________________________ MapInfo-L mailing list MapInfo-L@lists.directionsmag.com http://www.directionsmag.com/mailman/listinfo/mapinfo-l
