Aviation formulary: by Ed Williams V1.12 ----------------------------------- Great Circle Navigation Formulae This introduction is written for pilots (and others) who are interested in great circle navigation and would like to know how to compute courses, headings and other quantities of interest. These formulae can be programmed into your calculator or spreadsheet. I'll attempt to include enough information that those familiar with plane trigonometry can derive additional results if required. It is a well known that the shortest distance between two points is a straight line. However anyone attempting to fly from Los Angeles to New York on the straight line connecting them would have to dig a very substantial tunnel first. The shortest distance, *following the earth's surface* lies vertically above the aforementioned straight line route. This route can be constructed by slicing the earth in half with an imaginary plane through LAX and JFK. This plane cuts the (assumed spherical) earth in a circular arc connecting the two points, called a *great* circle. Only planes through the center of the earth give rise to great circles. An arbitrary plane will cut a sphere in a circle, but the resulting little circle is not the shortest distance between the points it connects. A little thought will show that lines of longitude (meridians) are great circles, but lines of latitude, with the exception of the equator, are not. I will assume the reader is familiar with latitude and longitude as a means of designating locations on the earth's surface. For the convenience of North Americans I will take North latitudes and West longitudes as positive and South and East negative. The longitude is the opposite of the usual mathematical convention. True course is defined as usual, as the angle between the course line and the local meridian measured clockwise. The first important fact to realise is that in general a great circle route has a true course that varies from point to point. For instance the great circle route between two points of equal (non-zero) latitude does not follow the line of latitude in an E-W direction, but arcs towards the pole. It is possible to fly between two points using an unvarying true course, but the resulting route differs from the great circle route and is called a *rhumb* line or loxodrome. Unlike a great circle which is a closed curve encircling the earth, a rhumb line spirals indefinitely poleward. Natural questions are to seek the great circle distance between two specified points and the true course at points along the route. The required spherical trigonometric formulae are greatly simplified if angles and distances are measured in the appropriate natural units, which are both radians! A radian, by definition, is the angle subtended by a circular arc of unit length and unit radius. Since the length of a complete circular arc of unit radius is 2*pi, the conversion is 360 degrees equals 2*pi radians, or: angle_radians=(pi/180)*angle_degrees ( and angle_degrees=(180/pi)*angle_radians ) Great circle distance can be likewise be expressed in radians by defining the distance to be the angle subtended by the arc at the center of the earth. Since by definition, one nautical mile subtends one minute (=1/60 degree) of arc, we have: distance_radians=(pi/(180*60))*distance_nm (and distance_nm=((180*60)/pi)*distance_radians ) In all subsequent formulae all distances and angles, such as latitudes, longitudes and true courses will be assumed to be given in radians, greatly simplifying them, and in applications the above formulae and their inverses are necessary to convert back and forth between natural and practical units. Examples of this process are given later. Some great circle formulae: The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by: d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2)) a mathematically equivalent formula, which is less subject to rounding error for short distances is: d=2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2)) Note: ^ denotes the exponentiation operator, sqrt is the square root function, acos the arc-cosine (or inverse cosine) function and asin is the arc-sine function. If asin or acos are unavailable they can be implemented using the atan2 function: acos(x)=atan2(sqrt(1-x^2),x) asin(x)=atan2(x,sqrt(1-x^2))} [Note: Here atan2 has the conventional (C) ordering of arguments, namely atan2(y,x). This is not universal, Excel for instance uses atan2(x,y), but it has asin and acos anyway. Be warned.] Further note: if your calculator is so impoverished that only atan is available then use: atan2(y,x)=atan(y/x) x>0 atan2(y,x)=atan(y/x)+pi x<0, y>=0 atan2(y,x)=pi/2 x=0, y>0 atan2(y,x)=atan(y/x)-pi x<0, y<0 atan2(y,x)=-pi/2 x=0, y<0 atan2(0,0) is undefined and should give an error. Note on the mod function. This appears to be implemented differently in different languages. Mod(y,x) is the remainder on dividing y by x and always lies in the range 0 <=mod =0 mod=y- x*int(y/x) else mod=y+ x*(int(-y/x)+1) endif ------------- The initial course, tc1, (at point 1) from point 1 to point 2 is given by: if sin(lon2-lon1)<0 tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1))) else tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1))) endif An alternative formula, not requiring the pre-computation of d, the distance between the points, is: tc1=mod(atan2(sin(lon1-lon2)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon1-lon2)), 2*pi) Intermediate points {lat,lon} lie on the great circle connecting points 1 and 2 when: lat=atan((sin(lat1)*cos(lat2)*sin(lon-lon2) -sin(lat2)*cos(lat1)*sin(lon-lon1))/(cos(lat1)*cos(lat2)*sin(lon1-lon2))) (not applicable for meridians. i.e if sin(lon1-lon2)=0) A point {lat,lon} is a distance d out on the tc radial from point 1 if: lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc)) if (cos(lat)=0) lon=lon1 // endpoint a pole else lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi endif This algorithm is limited to distances such that dlonsqrt(A^2 + B^2)) "no crossing" else dlon = acos(C/sqrt(A^2+B^2)) lon3_1=mod(lon1+dlon+lon+pi, 2*pi)-pi lon3_2=mod(lon1-dlon+lon+pi, 2*pi)-pi endif _________________________ Cross track error: Suppose you are proceeding on a great circle route from A to B (course =crs_AB) and end up at D, perhaps off course. You can calculate the course from A to D (crs_AD) and the distance from A to D (dist_AD) using the formulae above. In terms of these the cross track error, XTD, (distance off course) is given by XTD =asin(sin(dist_AD)*sin(crs_AD-crs_AB)) (positive XTD means right of course, negative means left) Examples: --------- Suppose point 1 is LAX: (33deg 57min N, 118deg 24min W) Suppose point 2 is JFK: (40deg 38min N, 73deg 47min W) In radians LAX is (33+57/60)*pi/180=0.592539, (118+24/60)*pi/180=2.066470 and JFK is (0.709186,1.287762) The distance from LAX to JFK is d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2)) =acos(sin(0.592539)*sin(0.709186)+ cos(0.592539)*cos(0.709186)*cos(0.778708)) =acos(0.811790) =0.623585 radians =0.623585*180*60/pi=2144nm The initial true course out of LAX is: sin(-0.778708)=-0.702<0 so tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1))) =acos((sin(0.709186)-sin(0.592539)*cos(0.623585))/ (sin(0.623585)*cos(0.592535)) =acos(0.408455) =1.150035 radians =66 degrees An enroute waypoint 100nm from LAX on the 66 degree radial (100nm along the GC to JFK) has lat and long given by: 100nm = 100*pi/(180*60)=0.0290888 radians lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc)) =asin(sin(0.592539)*cos(0.020888) +cos(0.592539)*sin(0.020888)*cos(1.150035)) =asin(0.568087) =0.604180 radians =34 degrees 37 min N lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi =mod(2.066470- asin(sin(1.150035)*sin(0.020888)/cos(0.604180))+pi,2*pi)-pi =mod(2.034206+pi,2*pi)-pi =2.034206 radians =116 degrees 33 min W The great circle route from LAX to JFK crosses the 111degree W meridian at a latitude of: (111degrees=1.937315 radians) lat=atan((sin(lat1)*cos(lat2)*sin(lon-lon2) -sin(lat2)*cos(lat1)*sin(lon-lon1))/(cos(lat1)*cos(lat2)*sin(lon1-lon2))) =atan((sin(0.592539)*cos(0.709186)*sin(0.649553) -sin(0.709186)*cos(0.592539)*sin(-0.129154))/(cos(0.592539)*cos(0.709186) *sin(0.778708))) =atan(0.737110) =0.635200radians =36 degrees 24min ____________________ Cross track error: Suppose enroute from JFK to LAX you find yourself at (D) N34:30 W116:30, which in radians is (0.6021386,2.033309) (See earlier for LAX, JFK coordinates and course) From LAX to D the distance is: dist_AD=acos(sin(0.592539)*sin(0.6021386)+ cos(0.592539)*cos(0.6021386)*cos(2.066470-2.033309)) =0.02905 radians (99.8665 nm) From LAX to D the course is: crs_AD=acos((sin(0.6021386)-sin(0.592539)*cos(0.02905))/ (sin(0.02905)*cos(0.592539))) =1.22473 radians (70.17 degrees) At point D the cross track error is: xtk= asin(sin(0.02905)*sin(1.22473-1.15003)) =0.00216747 radians = 0.00216747*180*60/pi =7.4512 nm right of course ___________________- Example of the intersection calc (briefly): Let point 1 be REO (42.60N,117.866W)=(0.74351,2.05715)rad Let point 2 be BKE (44.84N,117.806W)=(0.782606,2.056103)rad The 51 degree (=0.890118rad) bearing from REO intersects with 137 degree (=2.391101rad) from BKE at (lat3,lon3): Then: dst12=0.039103 crs12=0.018996 crs21=3.161312 ang1=0.871122 ang2=0.770211 ang3=1.500667 dst13=0.02729 dst23=0.029986 lat3=0.760473 =43.5N lon3=2.027876 =116.2W at BOI! ------------ Some general spherical triangle formulae. A spherical triangle is on whose sides are all great circular arcs. Let the sides have lengths a,b and c radians, and the opposite angles be A, B and C radians. c A -------B \ | \ | \b |a \ | \ | \ | \C| \| (The angle at B is not necessarily a right angle) sin(a) sin(b) sin(c) ----- = ------ = ------ sin(A) sin(B) sin(C) cos(a)=cos(b)*cos(c)+sin(b)*sin(c)*cos(A) cos(b)=cos(c)*cos(a)+sin(c)*sin(a)*cos(B) cos(c)=cos(a)*cos(b)+sin(a)*sin(b)*cos(C) cos(A)=-cos(B)*cos(C)+sin(B)*sin(C)*cos(a) cos(B)=-cos(C)*cos(A)+sin(C)*sin(A)*cos(b) cos(C)=-cos(A)*cos(B)+sin(A)*sin(B)*cos(c) Some useful consequences of these are: tan(A)=sin(B)*sin(a)/(sin(c)*cos(a)-cos(B)*cos(c)*sin(a)) tan(B)=sin(C)*sin(b)/(sin(a)*cos(b)-cos(C)*cos(a)*sin(b)) tan(C)=sin(A)*sin(c)/(sin(b)*cos(c)-cos(A)*cos(b)*sin(c)) tan(a)=sin(b)*sin(A)/(sin(C)*cos(A)+cos(b)*cos(C)*sin(A)) tan(b)=sin(c)*sin(B)/(sin(A)*cos(B)+cos(c)*cos(A)*sin(B)) tan(c)=sin(a)*sin(C)/(sin(B)*cos(C)+cos(a)*cos(B)*sin(C)) given *any* three of {a,b,c,A,B,C} the remaining sides and angles can be found from the above formulae. Note that for a spherical triangle A+B+C is not pi (180 degrees) but greater. _____________________________________________________________________ Rhumb lines (loxodromes) Rhumb lines or loxodromes are tracks of constant true course. Except for N/S courses, they are not the same as great circles. When two points (lat1,lon1), (lat2,lon2) are connected by a rhumb line with true course tc : lon2-lon1=-tan(tc)*(log((1+sin(lat2))/cos(lat2))+log((1+sin(lat1))/cos(lat1))) (logs are "natural" logarithms to the base e) The dist, d between the points is given by: d= abs((lat2-lat1)/cos(tc)) (cos(tc) > 0) = mod(lon2-lon1,2*pi)*cos(lat1) (tc=3*pi/2, lat1=lat2) = mod(lon1-lon2.2*pi)*cos(lat1) (tc=pi/2, lat1=lat2) To figure the rhumb line course,tc, and distance, d, from (lat1,lon1) to (lat2,lon2) allowing for the possibility of date-line crossings and E-W courses (neither point a pole!): dlon_W=mod(lon2-lon1,2*pi); dlon_E=mod(lon1-lon2,2*pi); dphi=log((1+sin(lat2))/cos(lat2))-log((1+sin(lat1))/cos(lat1)); if (dlon_Wsqrt(tol)){ dlon=dphi*tan(tc); } else{ // along parallel dlon=sin(tc)*d/cos(lat1); } lon=mod(lon1-dlon+pi,2*pi)-pi; } Example: rhumb line course from LAX to JFK: LAX (0.592539,2.066470) and JFK is (0.709186,1.287762) dlon_W=mod(1.287762-2.066470,2*pi)=5.50448 dlon_E=mod(2.066470-1.287762,2*pi)=0.778708 dphi=log((1+sin(0.709186))/cos(0.709186))-log((1+sin(0.592539))/cos(0.592539)) =0.146802 tc=mod(atan2(0.778708,0.146802),2*pi)=1.38446 = 79.3 degrees d=abs((0.709186-0.592539)/cos(1.38446))=0.62965 = 2164.6nm Compare this with the great circle course of 66 degrees and distance of 2144 nm. ---------------------------------------------------------------------------- Wind Triangles In all formulae, all angles are in radians. Convert back and forth as in the Great Circle section. [This is unnecessary on calculators which have a "degree mode" for trig functions. Most programming languages provide only "radian mode".] angle_radians=(pi/180)*angle_degrees ( and angle_degrees=(180/pi)*angle_radians ) A further conversion is required if using degrees/minutes/seconds: angle_degrees=degrees+(minutes/60.)+(seconds/3600.) degrees=int(angle_degrees) minutes=int(60*(angle_degrees-degrees)) seconds=int(3600*(angle_degrees-degrees-60*minutes)) [ You may have a built-in HH <-> HH:MM:SS conversion to do this efficiently] Let CRS=course, HD=heading, WD=wind direction (from), TAS=True airpeed, GS=groundspeed, WS=windspeed. Units of the speeds do not matter as long as they are all the same. (1) Unknown Wind: WS=sqrt( (TAS-GS)^2+ 4*TAS*GS*(sin((HD-CRS)/2))^2 ) WD=CRS + atan2(TAS*sin(HD-CRS), TAS*cos(HD-CRS)-GS) (**) if (WD<0) then WD=WD+2*pi if (WD>2*pi) then WD=WD-2*pi ( (**) assumes atan2(y,x), reverse arguments if your implementation has atan2(x,y) ) (2) Find HD, GS SWC=(WS/TAS)*sin(WD-CRS) if (abs(SWC)>1) then "course cannot be flown-- wind too strong" else HD=CRS+asin(SWC) if (HD<0) HD=HD+2*pi if (HD>2*pi) HD=HD-2*pi GS=TAS*sqrt(1-SWC^2)-WS*cos(WD-CRS) endif Note: The purpose of the "if (HD<0) HD=HD+2*pi; if (HD>2*pi) HD=HD-2*pi" is to ensure the final heading ends up in the range (0, 2*pi). Another way to do this, with the MOD function available is: HD=MOD(HD,2*pi) (3) Find CRS, GS GS=sqrt(WS^2 + TAS^2 - 2*WS*TAS*cos(HD-WD)) WCA=atan2(WS*sin(HD-WD),TAS-WS*cos(HD-WD)) (*) CRS=MOD(HD+WCA,2*pi) (*) WCA=asin((WS/GS)*sin(HD-WD)) works if the wind correction angle is less than 90 degrees, which will always be the case if WS < TAS. The listed formula works in the general case. _______________________________________________________________________ Variation formula for US. I did a least squares polynomial fit to the NFDC airport database. x=latitude (N degrees) y=longitude (W degrees) var= variation (degrees) var= -65.6811 + 0.99*x + 0.0128899*x^2 - 0.0000905928*x^3 + 2.87622*y - 0.0116268*x*y - 0.00000603925*x^2*y - 0.0389806*y^2 - 0.0000403488*x*y^2 + 0.000168556*y^3 Continental US only, 3771 points, RMS error 1 degree All within 2 degrees except for the following airports: MO49 MO86 MO50 3K6 02K and KOOA (2454 130 h_Tr) = 15-.0019812*h(ft) C (h < 36089.24ft) Variation of pressure with altitude: p= p_0*(1-6.87535*10^-6*h)^5.2561 h<36,089.24ft p_Tr= 0.22336*p_0 p=p_Tr*exp(-4.80634*10^-5*(h-36089.24)) h>36,089.24ft Variation of density with altitude: rho=rho_0*(1.- 6.87535*10^-6* h)^4.2561 h<36,089.24ft rho_Tr=0.29707*rho_0 rho=rho_Tr*exp(-4.80634*10^-5*(h-36089.24)) h>36,089.24ft _____________ Relationship of pressure and indicated altitude: alt_set in inches, heights in feet P_alt_corr= 145447*(1- (alt_set/29.92)^.19026) or P_alt_corr= (29.92-alt_set)*1000 (simple approximation) P_alt= Ind_Alt + P_alt_corr ______________ Relationship of pressure and density altitude: D_Alt=P_alt+(T_s/T_r)*(1.-(T_s/T)^.23498) (Standard temp T_s and actual temp T in Kelvin) An approximate, but fairly accurate formula is: D_Alt=P_Alt+118.6*(T-T_s) where T and T_s may (both) be either Celsius or Kelvin [[ Density altitude example: Let pressure altitude (P_alt) be 8000 ft, temperature 18C. Standard temp (T_s) is given by T_s=15-.0019812*8000=-0.85C = (273.15-0.85)K=272.3K Actual temperature (T) is 18C=(273.15+18)K=291.2K Density altitude (D_Alt) = 8000+ (272.3/.0019812)*(1-(272.3/291.2)^.23498) =8000+ 2150=10150ft or approximately: Density Altitude=8000 +118.6*(18+0.85)=10236ft ]] _______________ Relationship of true and calibrated (indicated) altitude: TA= CA + (CA-FE)*(ISADEV)/(273+OAT) where TA= True Altitude above sea-level FE= Field Elevation of station providing the altimeter setting CA= Calibrated altitude= Altitude indicated by altimeter when set to the altimeter setting, corrected for calibration error. ISADEV= Average deviation from standard temperature from standard in the air column between the station and the aircraft (in C) OAT= Outside air temperature (at altitude) The above is more precise than provided by the E6B or similar. ____________________________________________- Mach numbers, true vs calibrated airspeeds etc. Mach Number (M) = TAS/CS CS = sound speed= 38.9793*sqrt(T+273.15) where T is the OAT in celsius. TAS is true airspeed in knots. Because of compressibility, the measured IAT (indicated air temperature) is higher than the actual true OAT. Approximately: IAT=OAT+K*TAS^2/7592 The recovery factor K, depends on installation, and is usually in the range 0.95 to 1.0, but can be as low as 0.7. Temperatures are Celsius, TAS in knots. Also: OAT = (IAT + 273.15) / (1 + 0.2*K*M^2) - 273.15 The airspeed indicator measures the differential pressure, DP, between the pitot tube and the static port, the resulting indicated airspeed (IAS), when corrected for calibration and installation error is called "calibrated airspeed" (CAS). For low-speed (M<0.3) airplanes the true airspeed can be obtained from CAS and the density altitude, DA. TAS = CAS*(rho_0/rho)^0.5=CAS/(1.- 6.87535*10^-6 DA)^2.1280 (DA<36,089.24ft) Roughly, TAS increases by 1.5% per 1000ft. When compressibility is taken into account, the calculation of the TAS is more elaborate: DP=P_0*((1+0.2*(IAS/CS_0)^2)^3.5 -1) M=(5*( (DP/P+1)^(2/7) -1) )^0.5 TAS= M*CS P_0 is is (standard) sea-level pressure, CS_0 is the speed of sound at sea-level, CS is the speed of sound at altitude, and P is the pressure at altitude. These are given by earlier formulae: P_0= 29.92"Hg =1013.25 mB = 2116.2lbs/ft^2 P= P_0*(1-6.87535*10^-6 PA)^5.2561, pressure altitude, PA<36,089.24ft CS= 38.9793*sqrt(T+273.15) where T is the (static/true) OAT in Celsius. CS_0=38.9793*sqrt(15+273.15)=661.67 knots [Example: CAS=250 knots, PA=10000ft, IAT=2C, recovery factor=0.8 DP=29.92*((1+0.2*(250/661.67)^2)^3.5 -1)= 3.098" P=29.92*(1-6.87535*10^-6 *10000)^5.2561= 20.576" M= (5*( (3.098/20.576 +1)^(2/7) -1) )^0.5= 0.452 Mach OAT=(2+273.15)/(1 + 0.2*0.8*0.452^2) - 273.15= -6.7C CS= 38.9793*sqrt(-6.7+273.15)=636.27 knots TAS=636.27*0.452=287.6 knots] ------------------ { Some notes on the origins of some of the "magic number constants in the preceeding section: 6.87535*10^-6 = T'/T_0 where T' is the standard temperature lapse rate and T_0 is the standard sea-level temperature. 5.2561 = Mg/RT_0 M is the (average) molecular weight of air, g is the acceleration of gravity and R is the gas constant. 0.22336 is the ratio of the pressure at the tropopause to sea-level pressure. 4.80634*10^-5 is Mg/RT_tr where T_tr is the temperature at the tropopause. 4.2561 = Mg/RT_0 -1 0.29707= ratio of the density at the tropopause to the density at SL (rho_0) 145447= T_0/T' .19026=RT_0/Mg 145447= T_0/T' 38.9793 =sqrt(gamma R T_0/M) } -------- HUMIDITY, WET & DRY BULBS etc: The relative humidity, f (as a fraction) is related to the temperature, T and dewpoint Td by: f= exp(17.27(Td/(Td+237.3)-T/(T+237.3))) and to the frostpoint temperature Tf by: f= exp(21.87(Tf/(Tf+265.5)-T/(T+265.5))) Temperatures are in Celsius. Multiply f by 100 if you want a percentage. The above are based on an empirical fit to the saturation vapor pressure of water due to O. Tetens in Zeitschrift fur Geophysik, Vol VI (1930), quoted in "Principles of Meteorological Analysis" by W. J. Saucier (Dover NY 1983). This fit is: e_s=6.11 * exp(bT/(T+a)) for the saturation vapor pressure e_s in mbar over water a=237.3, b=17.27 over ice a=265.5, b=21.87 An alternative slightly more accurate fit (over water) is: e_s = 6.10779 + T * (4.43652e-1 + T * (1.42894e-2 + T * (2.65064e-4 + T * (3.03124e-6 + T * (2.03408e-8 + (6.13682e-11 * T)))))) (from Lowe, JAM (1977), 103) Tables of Relative Humidity and Dewpoint vs Temperature and Wet Bulb Temperature can be found in "Introduction to Meteorology" by Franklyn Cole (Wiley NY 1975). Inverting this to find dewpoint in terms of temp and RH: Dewpoint Td=237.3/(1/(ln(f)/17.27+T/(T+237.3))-1) Frostpoint Tf=265.5/(1/(ln(f)/21.87+T/(T+265.5))-1) Given the wet bulb temperature Tw (C), the dry bulb temperature T (C), and the pressure, p in mbar one gets the (approximate) relative humidity and dewpoint by the following: ed= 6.11*exp(17.27*T/(T+237.3)) /* SVP at dry-bulb temp ew= 6.11*exp(17.27*Tw/(Tw+237.3)) /* SVP at wet-bulb temp wd=0.62197*ed/(p-ed) /* saturation mixing ratio at T ww=0.62197*ew/(p-ew) /* saturation mixing ratio at Tw w=(2500.0*ww-1.0046*(T-Tw))/(2500.0+1.81*(T-Tw)) /* mixing ratio f= w/wd /* relative humidity as a fraction e= p*w/(0.62197+w) /* vapor pressure (mb) Td=(237.3*log10(e)-186.527)/(8.286-log10(e)) /* the dewpoint (C) This uses the Tetens fit for the saturated vapor pressure and treat water vapor as an ideal gas, both of which are pretty good approximations. If you want better refer to the Smithsonian Meteorological Tables (Smithsonian Institute 1963) _________________________ A related formula gives the increase in effective density altitude due to humidity. It only addresses the reduction of air density, and not the effect on engine power output: Increase(ft)=0.267*RH*(T+273)*exp(17.3*T/(T+237))*(1-0.00000688*H)^(-5.26) RH (f above) is the relative humidity expressed as a fraction, T is the temperature in Celsius and H is the pressure altitude in feet. Examples are: SL/30C/100% -> 565' increase in DA 10000/5C/80% -> 124' increase in DA 5000/40C/80% -> 1004' increase in DA. In terms of the dewpoint, Td the formula is: Increase(ft)=0.267*(T+273)*exp(17.3*Td/(Td+237))*(1-0.00000688*H)^(-5.26) which clearly agrees with the above when T=Td and RH=1. ------------------------------------------------------------------------------ Bellamy's formula for the wind drift and (single) wind correction angle is as follows: Drift (nm) = 21500*(p2-p1)/(sin(latitude)*TAS) (p2-p1 in inches) = 635 *(p2-p1)/(sin(latitude)*TAS) (p2-p1 in mB) Wind Correction Angle= 1230000*(p2-p1)/(sin(latitude)*TAS*Dist) (inches) = 36300* (p2-p1)/(sin(latitude)*TAS*Dist) (mB) p2-p1 is the difference between the destination and departure pressures. latitude is the average latitude on the route. TAS is the true airspeed in knots. Dist is the distance in nm. If the destination pressure is higher, the drift is to the left, and the required WCA is to the right (and vice-versa). Example: SFO -> LAX 300nm at 100 knots, latitude 36 degrees. Suppose the LAX altimeter setting is 0.2" higher (better the actual pressure difference at cruise altitude if you can get it). Drift = 21500*0.2/(sin(36)*100)= 73nm left WCA=1230000*0.2/(sin(36)*100*300)= 14 degrees right A discussion of this is in Barry Schiff's "Proficient Pilot I". ________________________________________________________ 1 knot = 1.15077 mph 1 mph = 0.86898 knot 1 knot = 1.85200 km/hr 1 km/hr= 0.53996 knot 1 mph = 1.60936 km/hr 1 km/hr= 0.62137 mph --------------------------------------------------------- Ellipsoidal parameters: Name Major axis (km) Ellipticity (f) WGS84 6378.13700 1/298.257223563 GRS80/NAD83 6378.13700 1/298.257222101 WGS66 6378.145 1/298.25 GRS67/IAU68 6378.16000 1/298.2472 WGS72 6378.135 1/298.26 Krasovsky 6378.245 1/298.3 Clarke66/NAD27 6378.2064 1/298.97866982 Reference: Coordinate Systems and Map Projections, D. H. Maling (Pergamon 1992) To convert between geocentric (radius r, geocentric latitude u) and geodetic coordinates (geodetic latitude v, height above the ellipsoid h): tan(u) = tan(v)*(h*sqrt((a*cos(v))^2+(b*sin(v))^2) +b^2)/ (h*sqrt((a*cos(v))^2+(b*sin(v))^2) +a^2) r^2 = h^2 + 2*h*sqrt((a*cos(v))^2+(b*sin(v))^2)+ (a^4-(a^4-b^4)*(sin(v))^2)/(a^2-(a^2-b^2)*(sin(v))^2) a and b are the semi-major axes of the ellipsoid b=a*(1-f) where f is the flattening Note that geocentric and geodetic longitudes are equal. ________________________________________________ Rate and radius of turn and pivotal altitude. In a steady turn, in no wind, with bank angle, b at an airspeed v tan(b)= v^2/(R g) and v= w R where g is the acceleration due to gravity, R is the radius of turn and w is the rate of turn. Pivotal altitude h_p is given by h = v^2/g With R in feet, v in knots, b in degrees and w in degrees/sec (inconsistent units!), numerical constants are introduced: R =v^2/(11.23*tan(0.01745*b)) At 100 knots, with a 45 degree bank, the radius of turn is 100^2/(11.23*tan(0.01745*45))= 891 feet. The rate of turn w is given by: w = 96.7*v/R (Example = 96.7*100/891= 10.9 degs/sec ) The bank angle b_s for a standard rate turn is given by: b_s = 57.3*atan(v/362.1) (Example) for 100 knots, b_s = 57.3*atan(100/362.1) = 15.4 degrees The pivotal altitude is given by: h_p = v^2/11.23 (Example) At 100 knots groundspeed the pivotal altitude is 100^2/11.23 = 890 feet. ---------------------------- Comments, corrections, suggestions to: Ed Williams CIS 72347,1516 72347.1516@compuserve.com ________________________________________________ CHANGES 1.04 (11/11/95) Add formula for computing lat/long of a point a given radial and distance valid when the distance can exceed one quarter of the earth's circumference. Note that atan2(0,0) should return an error. Add rhumb line formulae and example. Change intersection calculation to only provide result when intersection of radials exists. 1.05 (12/17/96) Correct test for pole in formula for computing lat/long of a point a given radial and distance: lat=0 => cos(lat)=0 1.06 (3/4/97) Definitions of a and b reversed in Tejens' fit to SVP. 1.07 (4/1/97) Add additional spherical triangle formulae. Correct the condition (dlon