Calculates the latitude and longitude of a specified number of points on a great circle between a given origin and destination. GREAT CIRCLE Program Author: David L. Deever DDEEVER@BCLCL1.IM.BATTELLE.ORG Otterbein College Westerville, Ohio 43081 (614) 891-2041 Date: April, 1993 This program asks the user to enter the latitude and longitude of an origin and destination. These latitudes and longitudes are entered as pairs of numbers in parentheses. They may also be stored as variables to be used as input. Easterly longitudes and southerly latitudes are designated by negative numbers. (This scheme follows that of J. Power's program BEAMHEAD.) The user is then asked for the desired number of segments. The program displays the length (in statute miles) of each of the segments. (They are all the same length.) It then proceeds to display the latitude and longitude of each of the points where the segments join and also gives the heading crossing the point. On a true great circle, of course, the heading is continuously varying. The program pauses after the display of each point. The command (40+7/60,82+56/60)->Wstrvll stores the values corresponding to 40o7' North latitude and 82o56' West longitude into the variable . Significant portions of a sample run could then be: Origin: Wstrvlle --could come from VARS menu Dest: (33,97) --approximately DALLAS # of segments: 5 ---------------------------------- Segment Length = 184.06 ---------------------------------- Point # = 0 Lat. = 40o7'0" Long. = 82o56'0" Hdg. = 242o9'11" ---------------------------------- Point # = 1 Lat. = 38o49'54" Long. = 85o57'34" Hdg. = 240o13'44" ---------------------------------- etc. Each dashed line represents a pause. There are two special cases to consider. The first is simple; if the origin and destination are within 0.05 miles of each other (about 2.6" of latitude) they are considered the same point and no route is calculated. The second case is more involved; if the origin and destination are within 0.05 miles of being antipodal points (on direct opposite sides of the earth) then any great circle route from the origin will arrive at the destination. In this case the user is asked to supply an original heading and that route is used. Segments are restricted to central angles of less than 90 degrees. This program may be used to plot a complete orbit of the earth by choosing antipodal points and giving an initial heading to get the first half and then using the same points and giving the opposite heading (+ or - 180o) to get the second half in reverse order. The points so calculated will not, of course, take into account the rotation of the earth under an orbiting sattelite. OUTLINE OF CALCULATION METHOD The calculations in this program make extensive use of the vector operations supplied with the TI-85. The calculations are simplified by placing the origin point at longtide 0 for internal calculation and taking the radius of the earth to be 1 unit. VA is a unit vector to the origin point. VB is a unit vector to the destination. AN is a unit vector tangent to the surface at the origin and pointing to the north polar axis. DD is the central angle between the origin and destination. DC = dot(VA,VB) DD = \cos^-1\(DC) DD/NN is the central angle of each segment AP is a unit vector to the current point DP = tan(DD/NN) is the length of a vector tangent to the surface from AP to the extended AP vector of the next point. (dot(VB,AP)/norm (AP)*AP is the projection of VB on AP BT = VB-(dot(VB,AP)/norm (AP)*AP is a tangent vector from AP toward VB Special case for antipodal points: The first BT = (cos AZ)*AN+(sin AZ)*cross(VA,AN) where AZ is the user entered initial heading. AP (new value ) = unitV (AP+(DP/norm (BT))*BT) The latitude of AP is given by \sin^-1\(AP(3)) The longitude of AP (from VA) is given by angle(AP(2),AP(1)) ( angle is one of the functions.) The formula cos(HDG) = (cos(A)sin(B)-sin(A)cos(B)cos(LoA-LoB)/sin(D) where A = start latitude B = ending latitude LoA-LoB = longitude diff from end to start D = central angle between start and end If LoA-LoB then HDG = 360-HDG is used to calculate each heading. Several minor features of the program are worth noting. The Degree/Radian state of the computer is saved and restored. Also, by clearing the screen, first displaying a value (with Disp) and then placing the annotation (with Outpt) we can get annotation and value on the same line. Outpt cannot display a value in degrees, minutes and seconds. (Note to TI: an improvement could be made here.) Care must be taken with inverse trig functions \cos^-1\ and \sin^-1\. If, even due to roundoff error, the argument to either of these functions is outside the range [-1,1] the TI-85 produces a type result (correctly) and goes merrily on its way. (Another note to TI: this fact should probably be documented in the instruction manual.)