From: spamless@nil.nil Newsgroups: sci.math,sci.physics,sci.space.science,sci.space.tech Subject: Re: Angle of geosync orbit to surface Date: 27 Jan 98 12:27:52 GMT In sci.math Chris Wagner wrote: > Hi, I'm trying to find what east-west angle a satellite orbiting > geosynchronusly over a given longitute makes with the surface at another > longitude. I'm trying to figure out what angle a satellite dish has to > be pointed at to pick up the satellite. The satellite is orbiting at > 119 west longitude. I am at 84.6 west longitude. and latitude? (geostationary satellites are over the equator) > I've spent hours > searching the web for a formula I could use with no success. Here you go! Distance/direction on the earth (spherical model) ------------------------------------------------- We will take the earth as a sphere of radius=1 (distances measured in earth radii). Start with geocentric spherical coordinates (rho, theta, phi; rho the distance from the center; theta the longitude measured counter-clockwise, viz. east is positive; phi the colatitude, zero at the north pole, 180_degrees at the south pole, the complement of the latitude=90_degrees-latitude, where latitude is positive in the northern and negative in the southern hemisphere) and consider geocentric rectangular coordinates (positive z-axis being the axis of rotation, viz. goes through the north pole; x-axis from the earth center out through latitude zero [or colatitude=90_degrees], longitude zero; y-axis from the center out through latitude zero, longitude 90_degrees east). For points on the earth's surface, rho=1, however I will leave rho in the equations so we can use it to get distance/direction to, say, a satellite (in case you want to use this to determine the direction to point your satellite dish). Our equations to convert from spherical to rectangular coordinates (for two points; i=1 or 2) are: x(i)=rho(i)*cos(theta(i))*sin(phi(i)) y(i)=rho(i)*sin(theta(i))*sin(phi(i)) z(i)=rho(i)*cos(phi(i)) We are interested in the direction and distance from point 1 (i=1) to point 2 (i=2) in terms of local coordinates at point 1. We will need the local unit vectors giving 1) UP, 2) NORTH and 3) EAST at a point on the surface (i=1). The local UP vector (local unit vector in the direction of increasing rho) is nothing but a unit vector pointing in the same direction as the vector from the center to the local point, viz. is given by (lu1=local up at point 1) (in terms of rectangular coordinates): lu1 (local up vector) [cos(theta(1))*sin(phi(1)),sin(theta(1))*sin(phi(1)),cos(phi(1))] (for simplicity I am writing vectors as row vectors, instead of column vectors, which I would have preferred) The local north vector (local unit vector in direction of decreasing phi) is a linear combination of lu1 and the unit z vector ([0,0,1]) which is orthogonal to lu1 and has a positive z component. One can work this out (take a*lu1+[0,0,b] as the vector and take its dot product with lu1 which is a*1+b*cos(phi(1)), for this to be zero, we can take a=-1, b=sec(phi(1)) so that the z coordinate of a*lu1+[0,0,b] is sec(phi(1))-cos(phi(1))=sin(phi(1))*tan(phi(1)). The length of the resulting vector is SQR[sin^2(phi(1))+sin^2(phi(1))*tan^2(phi(1))] or ABS[sin(phi(1))]*SQR[1+tan^2(phi(1))] or ABS[sin(phi(1))]*ABS(sec(phi(1))] = ABS(tan(phi(1)). However, this vector has component sin(phi(1))*tan(phi(1)) as its z-component and this is negative in the southern (phi>90_degrees) hemisphere. We wish to normalize, viz. divide by its length, and keep a positive z-component. If we divide each coordinate of -lu1+[0,0,sec(phi(1))] by tan(phi(1)) (rather than the absolute value of tan(phi(1))), the z component is sin(phi(1)) which IS positive for phi from 0 to 180 degrees! Note that the division only changes the sign of the first two components of lu1 and changes the sin(phi) to a cos(phi) and for the last component changes the cos(phi) to sin(phi) so that ln1, the local north vector at point 1 is given by:) ln1 (local north vector) [-cos(theta(1))*cos(phi(1)),-sin(theta(1))*cos(phi(1)),sin(phi(1))] (exercise: check that this is indeed orthogonal to lu1) The local east vector (le1) (local unit vector in the direction of increasing theta) is orthogonal to lu1 and ln1 (the z-axis is a linear combination of these two, so this vector has no z-component). Forgetting the z components (as this vector has no z-component) it is orthogonal to (x,y components only) lu1 [cos(theta(1),sin(theta(1))]*sin(phi(1)) and so is (+/-) [-sin(theta(1)),cos(theta(1))] (this is normalized!) Note that choosing this vector (plus sign) does give a vector which is along the y-axis ([0,1]) at the point theta1=0 and is along the negative y-axis for theta1=180_degrees, i.e. this is in the east direction. Thus the local east vector is given by (with z-component zero): le1 (local east vector) [-sin(theta(1)),cos(theta(1)),0] -=-=-=-=-=- Now the calculations are easy: --- (see SIMPLIFICATION below) --- Given two points (not necessarily on the earth's surface, i.e. rho need not be one!), calculate their rectangular coordinates (x(i), y(i) and z(i)) and the displacement vector from point 1 to point 2, P(1)=[x(1),y(1),z(1)], P(2)=[x(2),y(2),z(2)] disp=P(2)-P(1)=[x(2)-x(1),y(2)-y(1),z(2)-z(1)] with dist = mag(disp) = magnitude of the displacement vector = distance. and let dir=disp/dist (a unit vector giving the direction from point 1 to point 2). Now use the local up, north and east vectors to get the spherical coordinates of point two relative to point one as follows: The altitude (direction up or down, alt) is given by: sin(alt) = dir.lu1 (dot product) We can get components of dir along the north and east vectors by taking: north=dir.ln1 and east=dir.le1 (dot products) to get the azimuth (azi)) (and I will calculate this as an angle measured clockwise from north) as cos(azi)=north/cos(alt), sin(azi)=east/cos(alt) (the two dimensional vector [north/cos(alt),east/cos(alt)] is a unit vector) (alt is between -90 and +90 degrees so cos(alt) is positive) (this is nothing but a conversion between rectangular and polar coordinates) (note that this does not take into account the singular case where cos(alt)=0, viz. the object is directly overhead where the altitude is 90_degrees and the northerly and easterly components are undefined... I will ignore those cases) We then convert the rectangular coordinates (north,east) to polar coordinates (radius,angle) with angle being the azimuth/bearing to get the azimuth. (note that I have not given the formula for the great circle distance from the point on the earth's surface directly above/below point_1 -- if point_1 is not on the earth's surface -- to the point directly above/below point_2, but that is easily calculated from the angle between the vertical vectors at points 1 and 2... and is handled in the simplifications below) -=-=-=-=-=- Simplification Note that the position vectors P(i)=[x(i),y(i),z(i)] are given by rho(i)*lu_i (lu2 being the local vertical vector at position 2) and the displacement vector is given by P(2)-P(1)=rho(2)*lu2-rho(1)*lu1 so that the up/north/east components of the displacement in terms of the local directions at point 1 are given by (since lu1 is orthogonal both to ln1 and le1 and lu1 has length one): NORTHerly component = rho(2)*lu2.ln1 EASTerly component = rho(2)*lu2.le1 VERTical component = rho(2)*lu2.lu1-rho(1) Let cos(X)=lu2.lu1 (dot product) (the angle from one location to the other as measured at the center of the earth) (angular distance between the locations) lu2 given (of course) by: [cos(theta(2))sin(phi(2)),sin(theta(2))sin(phi(2)),cos(phi(2))] (in particular, normalizing the two dimensional vector of the NORTHerly and EASTerly components, which eliminates the rho(2) scale factor, shows that only the altitude, not the azimuth, depends upon the altitudes/rhos) Thus, in the calculations above, we can replace disp (the displacement) vector with P(2)=[x(2),y(2),z(2)] and use the modified: cos(X) = lu2.lu1 (dot product) We can use a different formula for X, namely: 2sin(X/2) = ||lu2-lu1|| (where the norm of a vector, v, is given by ||v|| = sqrt(v.v)) which can be more accurate for calculating X than taking arccosine(lu2.lu1) when the corresponding points on the sphere are very close together, for in that case, lu2.lu1 is very close to 1. dist=SQR(rho(2)^2+rho(1)^2-2*rho(1)*rho(2)*cos(X)) (law of cosines) sin(alt) = [rho(2)*cos(X)-rho(1)]/dist north=lu2.ln1 (dot product) east =lu2.le1 (dot product) (as we are only using north and east to get the direction, we can scale them, i.e. drop the rho2 term in the NORTHerly and EASTerly components of the displacement and find the azimuth/bearing, the angle of the ray passing from the origin through the point (north,east)) We now find the bearing to pass through the point (north,east) (convert from the rectangular coordinates, (north,east), to the polar coordinates, (radius,angle) with angle being the azimuth). Finally, we can get the sub-location great circle distance (here I use the term subsatellite point to refer to the point directly below, for example, a satellite in geostationary orbit... here we calculate the great circle distance along the earth's surface from the point on the earth's surface directly above/below POINT_1, if that point is not on the earth's surface, to the point on the earth's surface directly above/below POINT_2 if that point is not on the earth's surface) in terms of the angle X, since this great circle distance is simply X (for X measured in radians) (units: earth radius = 1). The QBasic programme at the end uses these simplified calculations (using P(2) and the rhos instead of the displacement). Spherical Trigonometry One may notice that no "Spherical Trigonometry" has been used. One can recover the Spherical Trigonometry results by expanding the results, above, in terms of the trigonometric function representations for the rectangular coordinates (e.g. one might write lu1.lu2 as cos(theta1-theta2)*sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2)=cos(X)). One can, using trig identities derive spherical trig. from the rectangular coordinates and the angles/vectors above. However, it seems simplest to use the rectangular coordinates in the QBasic programme. --- See the QBasic programme at the end --- Examples: The following results were all calculated using the QBasic programme at the end. (NY City: 40.75 degrees North (+) latitude and 74 degrees West (-) longitude on the earth's surface) From NY City to London (51.5 degrees North, 0 degrees West): Drilling a tunnel at an angle of 25 degrees down (below the horizon) in direction 51.21 degrees East of North (azimuth=51.21_degrees) will get you to London in 3355 miles, while traveling overland (great circle) starting at an azimuth of 51.21_degrees will take 3465 miles. From NY City to a geostationary satellite over a little North of Quito, Ecuador (0 degrees North, 78.5 degrees West): (the programme inserts a value of zero for the latitude for geostationary satellites) Ignoring refraction, the satellite is 42.7 degrees above the horizon at 173 degrees West of North (azimuth 360-173=187 degrees, seven degrees West of South) The distance to the satellite is 23344 miles, but one would only have to travel 2830 miles overland (great circle route) to be under the satellite. From NY City to San Francisco (37.7775 degrees North, 122.4111 degrees West): A tunnel 18.56 degrees down at 78.4 degrees West of North (azimuth = 360-78.4 = 281.6 degrees: remember compass directions are measured east/clockwise from North, so an angle given as a certain number of degrees West of North is negative, so we add it/subtract its absolute value from 360 degrees if one wants the bearing given in terms of the angle measured clockwise from North) will reach San Francisco in 2521 miles (straight line) while being forced to follow the earth's surface increases the distance (great circle) to 2565 miles. From NY City to Sydney, Australia (33.88 degrees South, 151.2 degrees East): A tunnel 71.9 degrees down at 93.9 degrees West of North (azimuth = 266.1 degrees) will reach Sydney in 7526 miles, while going around the earth will take 9937 miles. From NY City to a point at the bottom of a deep well (very deep... 2000 miles deep!) under Sydney, Australia: In the programme, set to miles as the units, enter an altitude of -2000 for the destination point. The bearing (azimuth) remains the same, as does the distance along the great circle to the point at the top of the well (Sydney) but the location (bottom of the well) is closer to the center of the earth, so we dig down in a more vertical direction (78.21 degrees down) and reach the well bottom well before we come up at Sydney (if we travel by tunnel) in 5659 miles. -=-=-QBasic Programme-=-=- Do NOT get this with word wrap set! One of the lines in this programme is 80 characters long (namely, the test for a geostationary satellite to set the latitude to zero). To use it, make a copy of this file to have the text description of the calculations. From another copy delete all lines up through (and including) the CUT_HERE line. Save the remaining lines (the QBasic programme lines only). The programme can then be loaded/run from QBasic. ----- CUT_HERE ----- DECLARE FUNCTION lawofcos! (a AS SINGLE, b AS SINGLE, cosc AS SINGLE) DECLARE FUNCTION gkey$ () DECLARE FUNCTION dot! (a() AS SINGLE, b() AS SINGLE) DECLARE FUNCTION acs! (c AS SINGLE) DECLARE FUNCTION arg! (x AS SINGLE, y AS SINGLE) DECLARE FUNCTION asn! (s AS SINGLE) ' ------------------------------------------------------------------------ ' DIRection and DISTance (DIRDIST.BAS) between points on a spherical earth ' Programme: 13 April 1997 by John McGowan ' ------------------------------------------------------------------------ ' --------------------------------------------------------------------- ' SLIGHTLY MODIFIED: 21 April: to calculate x ' from cosine(x) = lu2().lu1() <-- dot product ' to 2*sine(x/2) = ||lu2()-lu1()|| ' --------------------------------------------------------------------- ' --------------------------------------------------------------------- ' User definable fields (units and base site) (and DIM vectors) ' --------------------------------------------------------------------- ' Latitude given in decimal degrees ' North = + (positive), South = - (negative) ' Longitude given in decimal degrees ' East = + (positive), West = - (negative) ' RHO, the distance from the center of the earth to the location ' is given in earth radii in the data for the base site ' set up arrays for vectors DIM lu1(3), ln1(3), le1(3), lu2(3), disp(3) ' ----- UNITS ----- ' UNITS includes the scale factor to change from a chosen set of units ' to/from earth radii. You can change this scale factor by replacing ' the number with the number of kilometers, centimeters, inches, ' furlongs, etc. which gives the radius of the earth and changing the ' string to label the units chosen. ' Earth radius units: ' RADIUS UNITS RADIUS is the radius of the earth in UNITS DATA 3959.0, "miles" ' ----- SITE at which values are desired ----- ' This data is for your home site, POINT 1 in the text description ' so that one does not have continually to enter this data. You ' can change the base point, but remember to change the remark ' indicating the site and the text in the Data statement! ' New York City: Latitude, Longitude, rho, name of site (decimal degrees) ' Remember: East=+, West=-; North=+, South=- basept: ' LAT LON RHO SITE (rho=1.0, earth surface) DATA 40.75, -74.00, 1.000000, "New York City" ' --------------------------------------------------------------------- ' --------------------------------------------------------------------- ' --------------------------------------------------------------------- ' programme begins here (above is user setable header information) ' --------------------------------------------------------------------- Pi = 3.141593 ' set units to be used for distances and altitudes RESTORE units READ radius, radius$ ' set base point -- POINT 1: the site at which we want the direction RESTORE basept READ lat1, theta1, rho1, site1$ ' convert to radians lat1 = lat1 * Pi / 180!: theta1 = theta1 * Pi / 180! phi1 = Pi / 2! - lat1 'colatitude of base point: phi1 ' --------------------------------------------------------------------- ' data input routine to get data for the second site ' --------------------------------------------------------------------- getdata: CLS : PRINT PRINT "Data input for destination site to which one desires the distance" PRINT "from "; site1$; ".": PRINT PRINT "In this data input section," PRINT "DEFAULTs are values which will be used if you do NOT use" PRINT "one of the defined keys. If the defined keys are [E/W] for" PRINT "east/west with a default of W, then pressing ANY key other" PRINT "than e or E (lower case is useable) will result in using" PRINT "the default value of W." ' Altitude -- RHO=earth_radius+altitude (in units of earth_radius=1) PRINT : PRINT "Is the destination on the earth's surface? [Y/N] - default=Y" sur$ = gkey$ IF ((sur$ = "n") OR (sur$ = "N")) THEN PRINT "Is this a geostationary satellite? [Y/N] - default=Y" geo$ = gkey$ IF (geo$ = "n" OR geo$ = "N") THEN PRINT "Input the altitude above the earth in "; radius$; ".": INPUT r rho2 = (r + radius) / radius ELSE rho2 = 6.61517 END IF ELSE rho2 = 1! END IF ' we keep the inputs in degrees in variables_dm2 for checking later ' Latitude label1: ' Geostationary satellites are over the equator IF (((sur$ = "n") OR (sur$ = "N")) AND ((geo$ <> "n") AND (geo$ <> "N"))) THEN latdm2 = 0!: hemn$ = "N": lat2 = 0!: phi2 = Pi / 2! ELSE PRINT : PRINT "What is the absolute value of the latitude in" PRINT "decimal degrees [NO SIGN!]": INPUT latdm2 IF latdm2 < 0! THEN CLS : PRINT : PRINT : GOTO label1 lat2 = latdm2 PRINT "Is this in the Northern or Southern Hemisphere? [N/S] - default=N" hemn$ = gkey$: IF ((hemn$ = "s") OR (hemn$ = "S")) THEN lat2 = -lat2 lat2 = lat2 * Pi / 180!: phi2 = Pi / 2! - lat2 END IF PRINT ' Longitude label2: PRINT "What is the absolute value of the longitude in" PRINT "decimal degrees? [NO SIGN!]": INPUT thetadm2 IF thetadm2 < 0! THEN CLS : PRINT : PRINT : GOTO label2 theta2 = thetadm2: theta2 = -theta2 'default is Western/negative PRINT "Is this in the Western or Eastern Hemisphere? [W/E] - default=W" hemw$ = gkey$ IF ((hemw$ = "e") OR (hemw$ = "E")) THEN theta2 = -theta2 theta2 = theta2 * Pi / 180! ' --------------------------------------------------------------------- ' Input data check ' --------------------------------------------------------------------- CLS : PRINT : PRINT PRINT "The Data you have entered is (in decimal degrees):": PRINT NS$ = "N": IF (hemn$ = "S") OR (hemn$ = "s") THEN NS$ = "S" EW$ = "W": IF (hemw$ = "e") OR (hemw$ = "E") THEN EW$ = "E" PRINT "Latitude: "; latdm2; NS$: PRINT PRINT "Longitude: "; thetadm2; EW$: PRINT IF (sur$ <> "n") AND (sur$ <> "N") THEN PRINT "On the earth's surface." ELSEIF (geo$ <> "n") AND (geo$ <> "N") THEN PRINT "Geostationary satellite orbit (over the equator)." ELSE PRINT "at altitude "; r; " "; radius$; "." END IF PRINT : PRINT "Is this correct? [Y/N] - default=Y": check$ = gkey$ IF ((check$ = "n") OR (check$ = "N")) THEN GOTO getdata ' --------------------------------------------------------------------- ' Calculate local vectors ' --------------------------------------------------------------------- ' local up (vertical) vectors lu1(1) = COS(theta1) * SIN(phi1): lu2(1) = COS(theta2) * SIN(phi2) lu1(2) = SIN(theta1) * SIN(phi1): lu2(2) = SIN(theta2) * SIN(phi2) lu1(3) = COS(phi1): lu2(3) = COS(phi2) 'local north vector ln1(1) = -COS(theta1) * COS(phi1) ln1(2) = -SIN(theta1) * COS(phi1) ln1(3) = SIN(phi1) 'local east vector le1(1) = -SIN(theta1) le1(2) = COS(theta1) le1(3) = 0! ' --------------------------------------------------------------------- ' Calculations ' --------------------------------------------------------------------- ' cosine of angle between points cosx = dot(lu2(), lu1()) 'dot product ' angle between points in RADIANS ' (this is the great circle distance between the points on the ' earth's surface above/under the sites -- if they are not on ' the earth's surface -- in units of earth_radius=1) ' x = acs(cosx) ' however, in hopes of getting a slightly more accurate value ' for x when the points are very close together (and cosx is very ' close to 1) we use 2*SIN(x/2)=||disp()|| where disp()=lu2()-lu1() ' is the displacement vector between the corresponding points on the ' surface of the earth. disp(1) = lu2(1) - lu1(1) disp(2) = lu2(2) - lu1(2) disp(3) = lu2(3) - lu1(3) x = 2! * asn(SQR(dot(disp(), disp())) / 2!) ' law of cosines to get distance (sides rho1, rho2 included angle x) ' between sites, whether or not on the earth's surface, in units ' of earth_radius=1 dist = lawofcos(rho1, rho2, cosx) sinalt = (rho2 * cosx - rho1) / dist ' scaled arcsine to get altitude in DEGREES alt = 180! * asn(sinalt) / Pi north = dot(lu2(), ln1()) 'dot product east = dot(lu2(), le1()) 'dot product ' scaled azimuth (bearing) in DEGREES azi = 180! * arg(north, east) / Pi ' --------------------------------------------------------------------- ' Output results ' --------------------------------------------------------------------- CLS : PRINT : PRINT PRINT "The distance from "; site1$; " to the site at" PRINT "Latitude "; latdm2; NS$; " and Longitude "; thetadm2; EW$ IF (sur$ <> "n") AND (sur$ <> "N") THEN PRINT "on the earth's surface is:" ELSEIF (geo$ <> "n") AND (geo$ <> "N") THEN PRINT "in geostationary satellite orbit (over the equator) is:" ELSE PRINT "at altitude "; r; " "; radius$; " is:" END IF PRINT : PRINT dist; " earth radii or"; dist * radius; " "; radius$: PRINT up$ = "above": IF alt < 0! THEN up$ = "below" PRINT "in a straight line at angle" PRINT ABS(alt); " (in decimal degrees) "; up$; " the horizon" PRINT : east$ = "East": IF azi < 0! THEN east$ = "West" PRINT "and at angle" PRINT ABS(azi); " (in decimal degrees) "; east$; " of North.": PRINT PRINT "The distance along a great circle route on the earth from "; site1$ PRINT "to the point directly below/above the chosen site (if it is not" PRINT "on the earth's surface) is:": PRINT PRINT x; " earth radii or"; x * radius; " "; radius$; ".": PRINT : PRINT PRINT "Do you want to find the distance to another site? [Y/N] - default=Y" more$ = gkey$ IF (more$ <> "n") AND (more$ <> "N") THEN GOTO getdata SYSTEM ' --------------------------------------------------------------------- ' END of main programme ' --------------------------------------------------------------------- FUNCTION acs (c AS SINGLE) ' ARCCOSINE of c in RADIANS ' Avoid floating point inaccuracies! (assume all errors are such). ' Use the half angle formula for the calculations. ' tan(pi/4-theta/2)=cos(theta)/(1+sin(theta)) where we calculate ' sin(theta) from c=cos(theta) as sin(theta)=sqrt(1-c^2) ' Use temp variable, t, and do NOT change the input value! t = c: IF ABS(c) >= 1! THEN t = SGN(c) * 1! acs = 3.141593 / 2! - 2! * ATN(t / (1 + SQR(1 - t * t))) END FUNCTION FUNCTION arg (x AS SINGLE, y AS SINGLE) ' Find the angle (in RADIANS) of the ray from the origin through the ' point x,y, i.e. the ARGument of the complex number x+iy. ' Use the half angle formula tan(theta/2)=sin(theta)/(1+cos(theta)) ' where cos(theta)=x/R and sin(theta)=y/R for R=sqrt(x^2+y^2). ' For accuracy, ensure that theta/2 is between either -Pi/4 and Pi/4 ' or (if x<0) rotate by 180_degrees to ensure this (then rotate back). ' arg(x,y) = atn2(y,x) IF (y = 0! AND x <= 0!) THEN arg = 3.141593 ELSEIF x < 0 THEN arg = 2! * ATN(-y / (ABS(x) + SQR(x * x + y * y))) + SGN(y) * 3.141593 ELSE arg = 2! * ATN(y / (x + SQR(x * x + y * y))) END IF END FUNCTION FUNCTION asn (s AS SINGLE) ' ARCSINE of s in RADIANS ' Avoid floating point inaccuracies! (assume all errors are such) ' Use a temporary variable, t, and do NOT change the input variable! ' Use half angle formula method ' tan(theta/2)=sin(theta)/(1+cos(theta)) where we calculate ' cos(theta) from s=sin(theta) as cos(theta)=sqrt(1-s^2). t = s: IF ABS(s) >= 1! THEN t = SGN(s) * 1! asn = 2! * ATN(t / (1! + SQR(1! - t * t))) END FUNCTION FUNCTION dot (a() AS SINGLE, b() AS SINGLE) ' dot product ' dot product of the vectors [a(1),a(2),a(3)] and [b(1),b(2),b(3)] dot = a(1) * b(1) + a(2) * b(2) + a(3) * b(3) END FUNCTION FUNCTION gkey$ ' routine to get one key stroke ' clear key buffer DO WHILE INKEY$ <> "": LOOP 'get key g$ = "": DO WHILE g$ = "": g$ = INKEY$: LOOP gkey$ = g$ END FUNCTION FUNCTION lawofcos (a AS SINGLE, b AS SINGLE, cosc AS SINGLE) ' Law of Cosines: c=sqrt(a^2+b^2-2*a*b*cos(C)) ' given a, b and the cosine of C ' Avoid floating point inaccuracies (assume errors in cos(C) are such). ' Use temporary variable t for cos(C) and do NOT change cos(C) globally t = cosc: IF ABS(cosc) >= 1! THEN t = SGN(cosc) * 1! lawofcos = SQR(a * a + b * b - 2! * a * b * t) END FUNCTION