Angle difference between two bearings: Difference between revisions

Rather than calculate angle differences and mess about with folding the results into ±180 and getting the sign right, why not use some mathematics? These days, trigonometrical functions are calculated swiftly by specialised hardware (well, microcode), and with the availability of functions working in degrees, matters are eased further nor is precision lost in converting from degrees to radians. So, the first step is to convert a bearing into an (x,y) unit vector via function CIS(t) = cos(t) + i.sin(t), which will handle all the annoyance of bearings specified in values above 360. Then, using the dot product of the two vectors allows the cosine of the angle to be known, and the cross product determines the sign.

However, this relies on the unit vectors being accurately so, and their subsequent dot product not exceeding one in size: given the rounding of results with the limited precision actual floating-point arithmetic, there may be problems. Proving that a calculation will not suffer these on a specific computer is difficult, especially as the desire for such a result may mean that any apparent pretext leading to that belief will be seized upon. Because calculations on the IBM pc and similar computers are conducted with 80-bit floating-point arithmetic, rounding errors for 64-bit results are likely to be small, but past experience leads to a "fog of fear" about the precise behaviour of floating-point arithmetic.

As it happens, the test data did not provoke any objections from the ACOSD function, but even so, a conversion to using ''arctan'' instead of ''arccos'' to recover angles would be safer. By using the four-quadrant ''arctan(x,y)'' function, the sign of the angle difference is also delivered and although that result could be in 0°-360° it turns out to be in ±180° as desired. On the other hand, the library of available functions did not include an arctan for complex parameters, so the complex number Z had to be split into its real and imaginary parts, thus requiring two appearances and to avoid repeated calculation, a temporary variable Z is needed. Otherwise, the statement could have been just <code>T = ATAN2D(Z1*CONJG(Z2))</code> and the whole calculation could be effected in one statement, <code>T = ATAN2D(CIS(90 - B1)*CONJG(CIS(90 - B2)))</code> And, since ''cis(t) = exp(i.t)'', <code>T = ATAN2D(EXP(CMPLX(0,90 - B1))*CONJG(EXP(CMPLX(0,90 - B2))))</code> - although using the arithmetic statement function does seem less intimidating.

COMPLEX*16 CIS,Z1,Z2,Z !Scratchpads.

Z = Z1*CONJG(Z2) !(Z1x,Z1y)(Z2x,-Z2y) = (Z1x.Z2x + Z1y.Z2y, Z1y.Z2x - Z1x.Z2y)

T = ATAN2D(AIMAG(Z),REAL(Z)) !Madly, arctan(x,y) is ATAN(Y,X)!

1 -45D0,125D0, -45D0,145D0, 29.4803D0,-88.6381D0,

3 -70099.74233810938D0, 29840.67437876723D0,

5 1174.8380510598456D0, -154146.66490124757D0,

WRITE (6,1) ("B",I,"x","y", I = 1,2) !Or, one could just list them twice.

* 2(A13,I1,"(",A9,",",A9,")"),A9) !Compare format 10, above.

B1( x, y) B2( x, y)

29.4803( 0.492124, 0.870525) -88.6381(-0.999718, 0.023767) -118.118