refrac.c 1.52 KB
 Jeroen Vreeken committed Oct 08, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 ``````/* Atmospheric refraction * Returns correction in degrees to be added to true altitude * to obtain apparent altitude. * * -- S. L. Moshier */ extern double atpress; /* millibars */ extern double attemp; /* degrees C */ extern double DTR; /* pi/180 */ #if __STDC__ double tan (double); #else double tan(); #endif double refrac(alt) double alt; /* altitude in degrees */ { double y, y0, D0, N, D, P, Q; int i; if( (alt < -2.0) || (alt >= 90.0) ) return(0.0); /* For high altitude angle, AA page B61 * Accuracy "usually about 0.1' ". */ if( alt > 15.0 ) { D = 0.00452*atpress/((273.0+attemp)*tan( DTR*alt )); return(D); } /* Formula for low altitude is from the Almanac for Computers. * It gives the correction for observed altitude, so has * to be inverted numerically to get the observed from the true. * Accuracy about 0.2' for -20C < T < +40C and 970mb < P < 1050mb. */ /* Start iteration assuming correction = 0 */ y = alt; D = 0.0; /* Invert Almanac for Computers formula numerically */ P = (atpress - 80.0)/930.0; Q = 4.8e-3 * (attemp - 10.0); y0 = y; D0 = D; for( i=0; i<4; i++ ) { N = y + (7.31/(y+4.4)); N = 1.0/tan(DTR*N); D = N*P/(60.0 + Q * (N + 39.0)); N = y - y0; y0 = D - D0 - N; /* denominator of derivative */ if( (N != 0.0) && (y0 != 0.0) ) /* Newton iteration with numerically estimated derivative */ N = y - N*(alt + D - y)/y0; else /* Can't do it on first pass */ N = alt + D; y0 = y; D0 = D; y = N; } return( D ); }``````