# Astronomical Refraction for the HP-41

Overview

1°) Standard Atmosphere

a) Apparent Altitude >>> True Altitude
b) Other Simple Formulae
c) True Altitude >>> Apparent Altitude

2°) More general Programs

a) 2 short routines
b) Middle programs
c) A long program

-The light coming from a star is curved by the atmosphere,  so that its apparent altitude differs from its true altitude.
-Knowing the refraction R allows to convert the apparent altitude h0 and the true altitude h of a given star:   h = h0 - R
-The following programs use data from the Pulkovo Refraction Tables.

1°)  Standard Atmosphere(s)

a) Apparent Altitude >>> True Altitude

-The short routine hereafter gives relatively accurate results in the following atmospheric conditions:

Temperature:                                +15°Celsius
Pressure:                                    1013.25 mbar
Light wave-length:                         0.590 µm
Partial pressure of water vapor:          0   ( dry air )
Latitude:                                           45°
Observer's altitude:                            0   ( i-e at sea-level )

Formula:    R ~ (1°/62.8093)/Tan( h0 + 4.2206/( h0 + 15.1115/( h0 + 5.9431 ) ) )       where h0 is expressed in degrees

-Errors are smaller than  0"29  over the whole range [ 0° , 90° ]

 01 LBL "H0H"  02 DEG  03 HR  04 15.1115  05 RCL Y  06 5.9431  07 +  08 /  09 +  10 4.2206  11 X<>Y  12 /  13 +  14 TAN  15 1/X  16 62.8093  17 /  18 ST- Y  19 HMS  20 X<>Y  21 HMS  22 END

( 54 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y / R X h0 h

-All angles expressed in ° ' "

Examples:     With  h0 = 1°30'00"  ;   h0 = 27°  ;   h0 = 0

1.30  XEQ "H0H"  >>>>   h =   1°09'42"6    RDN    R = 0°20'17"4
27        R/S          >>>>   h = 26°58'08"3    RDN    R = 0°01'51"7
0         R/S          >>>>   h = -0°32'57"9     RDN   R =  0°32'57"9

-For the horizontal refraction, "H0H" yields 1977"880
-According to Pulkovo refraction tables, it should be  1977"971

-Laplace formula  R = 57".085/Tan h0 - 0".0666/Tan3 h0  gives even better results if h0 > 20° , more precisely:

errors are smaller than   0"02   for h0 > 20°
errors are smaller than   0"01   for h0 > 23°
errors are smaller than  0"002  for h0 > 30°

b) Other Simple Formulae

-We can build similar formulae for other atmospheric conditions.

R  ~   a/Tan( h0 + b/( h0 + c/( h0 + d ) ) )       where h0 is expressed in degrees

-For a dry air, light wave-length = 0.59 µm , latitude = 45° at sea-level, we find:

 t (°C) P (mbar) a b c d E error(0) -30 1013.25 67"8716 3.3124 13.2157 4.9474 0"20 -80"0 -10 1013.25 62"720 3.713 14.049 5.389 0"28 -33"6 +10 1013.25 58"292 4.093 14.561 5.730 0"30 -4"9 +30 1013.25 54"451 4.469 15.077 6.089 0"32 +14"4 15 500 28"265 4.438 15.125 5.892 0"17 0"17 15 700 39"576 4.350 15.088 5.904 0"23 0"09 15 900 50"885 4.253 14.937 5.881 0"29 0"23 15 1100 62"187 4.145 14.609 5.801 0"32 0"02

-  E is the maximum error ( in absolute value ) over the interval  [ 0°10' ; 90° ]
-  error(0) is the error at the horizon ( compared with the Pulkovo tables )

-As you can see, error(0) is much greater than E when  t < 15°C or t > 15°C , especially for very low temperatures,
and I don't really know the behavior of the function R(h0)  for  0° < h0 < 0°10'  if t # 15°C

c) True Altitude >>> Apparent Altitude

-We now assume t = 15°C & P = 1013.25 mbar again with the other standard parameters above.

Formula:    R ~ (1°/62.644)/Tan( h + 5.409/( h + 18.732/( h + 6.807 ) ) )       where h is expressed in degrees

-Errors are smaller than  0"62  over the whole range [ -0°32'58" , 90° ]

 01 LBL "HH0"  02 DEG  03 HR  04 18.732  05 RCL Y  06 6.807  07 +  08 /  09 +  10 5.409  11 X<>Y  12 /  13 +  14 TAN  15 1/X  16 62.644  17 /  18 ST+ Y  19 HMS  20 X<>Y  21 HMS  22 END

( 50 bytes / SIZE 000 )

 STACK INPUTS OUTPUTS Y / R X h h0

-All angles expressed in ° ' "

Example:        With  h = 1°30'00"

1.30  XEQ "H-H0"  >>>>  h0 = 1°48'38"8    RDN    R = 0°18'38"8

2°)  More General Programs

a) 2 Short Routines

-We still assume a dry air, light wave-length = 0.59 µm , latitude = 45° , altitude = 0
but now, the temperature t may be different from 15°C and the pressure P may be different from 1013.25 mbar.

Formula:    R ~ (P/1013.25).(288.15/(t+273.15)).(1°/62.8093)/Tan( h0 + 4.2206/( h0 + 15.1115/( h0 + 5.9431 ) ) )       where h0 is expressed in degrees

Data Registers:    R00 = unused                              ( Registers R01 & R02 are to be initialized before executing "REF" )

•  R01 = t     ( in °C )
•  R02 = P    ( in mbar )
Flags: /
Subroutines: /

 01  LBL "REF" 02  DEG 03  HR 04  15.1115 05  RCL Y       06  5.9431 07  + 08  / 09  + 10  4.2206 11  X<>Y       12  / 13  + 14  TAN 15  1/X 16  220.8625 17  / 18  RCL 01 19  273.15 20  + 21  / 22  RCL 02       23  * 24  X<0? 25  CLX 26  ST- Y 27  HMS 28  X<>Y         29  HMS 30  END

( 68 bytes / SIZE 003 )

 STACK INPUTS OUTPUTS Y / R X h0 h

-All angles expressed in ° ' "

Example:       t = -10°C , P = 1100 mbar      -10  STO 01   1100   STO 02

-If   h0 =  12°34'56"    12.3456  XEQ "REF"  >>>>   h = 12°29'58"4   X<>Y   R = 0°04'57"6
-If   h0 =  10°12'34"    10.1234       R/S          >>>>   h = 10°06'29"5   X<>Y   R = 0°06'04"5

-The results are acceptable for  h0 > 10°  ( errors are of the order of 1 or 2 arcseconds for h0 = 10° )
-For heights of a few degrees, the outputs are much more doubtful.

-The following program takes the true altitude ( in decimal degrees ) and returns the apparent altitude ( in decimal degrees )

Formula:

h0 - h  ~  (P/1013.25).(288.15/(t+273.15)).(1°/62.644)/Tan( h + 5.409/( h + 18.732/( h + 6.807 ) ) )

Data Registers:    R00 = unused                                                            ( Registers R01 & R02 are to be initialized before executing "HH0" )

•  R01 = t     ( in °C )
•  R02 = P    ( in mbar )
Flags: /
Subroutines: /

 01 LBL "HH0"  02 DEG  03 18.732  04 RCL Y            05 6.807 06 +  07 /  08 +  09 5.409  10 X<>Y 11 /  12 +  13 TAN  14 1/X  15 220.281 16 /  17 RCL 01          18 273.15  19 +  20 / 21 RCL 02          22 *  23 +  24 END

( 57 bytes / SIZE 003 )

 STACK INPUT OUTPUT X h h0

-All angles expressed in decimal degrees.

Example:       t = 0°C , P = 1100 mbar        0  STO 01   1100   STO 02

if h = 20°     20  XEQ "HH0"  >>>     h0 = 20°0495
if h = 10°     10      R/S           >>>     h0 = 10°0987

b) Middle Programs

-We suppose a light wave-length = 0.59 µm , latitude = 45° , altitude = 0
but the temperature t , the pressure P and the humidity f may have non-standard values.

Formulae:   The standard refraction is computed by  R0 ~ (1°/62.8093)/Tan( h0 + 4.2206/( h0 + 15.1115/( h0 + 5.9431 ) ) )

R = R0 ( P/1013.25 ) [ 286.68 / ( t + 271.68 ) ] [ 1 - f / 6579 - ( f/426 )2 ] [ 1 + { (15-t) + (15-t)2 / 196 } exp(-h0/7) / (1+h0) / 157 ]
x [ 1 + (P-1013.25) exp(-0.4h0)/12181 ] [ 1 - ( f + f2/16 ) exp(-0.18h0) / (1+h0) / 5198 ]

Data Registers:    R00 = h0                                                                 ( Registers R01 thru R03 are to be initialized before executing "H0H" )

•  R01 = t     ( in °C )           ( -10 <= t <= +30 )
•  R02 = P    ( in mbar )      ( 700 <= P <= 1100 )
•  R03 = f     ( in mbar )          ( 0 <= f <= 20 )
Flags: /
Subroutines: /

 01  LBL "H0H" 02  DEG 03  HR 04  STO 00 05  15.1115 06  RCL 00         07  5.9431 08  + 09  / 10  + 11  4.2206 12  X<>Y 13  / 14  + 15  TAN 16  RCL 02 17  X<>Y 18  / 19  221.995 20  / 21  RCL 01         22  271.68 23  + 24  / 25  RCL 03 26  6579 27  / 28  RCL 03 29  426 30  / 31  X^2 32  + 33  * 34  - 35  15 36  RCL 01 37  - 38  ENTER^ 39  X^2 40  196 41  / 42  + 43  157 44  / 45  RCL 00         46  7 47  / 48  E^X 49  / 50  RCL 00 51  1 52  + 53  / 54  * 55  + 56  RCL 02 57  1013.25 58  - 59  12181 60  / 61  RCL 00 62  .4 63  * 64  E^X 65  / 66  * 67  + 68  RCL 03        69  RCL 03 70  4 71  / 72  X^2 73  + 74  5198 75  / 76  RCL 00 77  .18 78  * 79  E^X 80  / 81  RCL 00 82  1 83  + 84  / 85  * 86  - 87  RCL 00        88  X<>Y 89  ST- Y 90  HMS 91  X<>Y 92  HMS 93  END

( 155 bytes / SIZE 004 )

 STACK INPUTS OUTPUTS Y / R X h0 h

( All angles are expressed in ° ' " )   Execution time = 6 seconds.

Example:         t = 0°C ;  P = 900 mbar  ;  f = 12 mbar

0  STO 01    900  STO 02    12  STO 03

-If   h0 =  0°                     0      XEQ "REF2"  >>>>   h = -0°33'32"       X<>Y    R = 0°33'32"
-If   h0 =  10°23'45"    10.2345       R/S          >>>>   h = 10°19'02"8    X<>Y    R = 0°04'42"2
-If   h0 =  49°12'34"    49.1234       R/S          >>>>   h = 49°11'47"9    X<>Y    R = 0°00'46"1

-If  h0 > 20° more accurate results will be obtained if  R0 is computed by Laplace's formula.
-The accuracy is of the order of 10" near the horizon, but errors rapidly decrease as h0 increases.
-However, disturbances of the atmosphere make accurate results more theoretical than real for very small h0-values:
near the horizon, the refraction may fluctuate by several arcminutes !

Data Registers:    R00 = h                                                                 ( Registers R01 thru R03 are to be initialized before executing "HH0" )

•  R01 = t     ( in °C )           ( -10 <= t <= +30 )
•  R02 = P    ( in mbar )      ( 700 <= P <= 1100 )
•  R03 = f     ( in mbar )          ( 0 <= f <= 20 )
Flags: /
Subroutines: /

 01 LBL "HH0"  02 DEG  03 HR  04 STO 00  05 18.732  06 RCL 00            07 6.807  08 +  09 /  10 +  11 5.409  12 X<>Y  13 /  14 +  15 TAN  16 1/X  17 62.644  18 /  19 STO 06  20 +  21 STO 04 22 7  23 /  24 CHS  25 E^X  26 15  27 RCL 01            28 -  29 STO 05  30 196  31 /  32 RCL 05  33 +  34 *  35 1  36 RCL 04  37 +  38 STO 07  39 /  40 157  41 /  42 1 43 +  44 RCL 04  45 .4  46 *  47 CHS  48 E^X  49 12181  50 /  51 RCL 02            52 1013.25  53 STO 05  54 -  55 *  56 1  57 +  58 *  59 RCL 03  60 6579  61 /  62 RCL 03  63 426 64 /  65 X^2  66 +  67 CHS  68 1  69 +  70 *  71 RCL 03            72 RCL 03  73 4  74 /  75 X^2  76 +  77 5198  78 /  79 RCL 04  80 .18  81 *  82 CHS  83 E^X  84 * 85 RCL 07  86 /  87 CHS  88 1  89 +  90 *  91 286.677  92 *  93 RCL 02            94 *  95 RCL 05  96 /  97 RCL 01  98 271.677  99 + 100 / 101 RCL 06 102 * 103 RCL 00 104 + 105 HMS 106 END

( 170 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS X h h0

( All angles are expressed in ° ' " )   Execution time = 7 seconds.

Example:         t = 0°C ;  P = 900 mbar  ;  f = 12 mbar

0  STO 01    900  STO 02    12  STO 03

h = 1°                       1        XEQ "HH0"  >>>>    h0 = 1°20'18"7
h = 10°23'45"     10.2345        R/S          >>>>    h0 = 10°28'25"3

-The following program tries to produce the accurate values of the Pulkovo refraction tables, and the other parameters are now taken into account:

c) A Long Program

-This program uses the following formula to compute the atmospheric refraction in a standard atmosphere:

R0 ~ (1°/63.05561)/Tan( h0 + 3.81451/( h0 + 6.04529/( h0 + 8.42681/( h0 + 23.82074/( h0 + 7.40780 ) ) ) ) )

-Then, refraction R in more general conditions is calculated by:

R = R0 (1.0552126/(1+0.00368084.t)).(1+A).(P/1013.25).(1+B).(0.98282+0.005981/Lwl2).(1+C).(1-0.152 10 -3 f -0.55 10 -5 f 2 ).(1+D).(1+E).(1+F)

where

t = temperature (°C) ; P = pressure (mbar) ; f = partial pressure of water vapor (mbar)
Lwl = Light wave-length ( µm ) ; lat = latitude
alt = observer's altitude ( over sea-level ) in meters.

- A , B , C , D , E , F are coefficients which depend on the temperature, pressure, Light wave-length, humidity, latitude and observer's altitude respectively.

-They are very small near the zenith, but they can't be neglected near the horizon!
-Reference [1] provides tables for these coefficients, but "REFR" uses approximate formulae instead:

-Lagange's interpolation formula is used to obtain A and B for other t- and P-values.

Coefficient A:         With  x = 1/(1+h0)

-If t = -30°C    105A = Max ( -2 - 1411 x + 100967 x2 + 3583 x3 - 465432 x4 + 928890 x5 - 783471 x6 + 251549 x7 + 2377 e -43.h0 ; 0 )
-If t = -10°C    105A = Max (  -880 x + 57082 x2 - 6928 x3 - 250807 x4 + 515833 x5 - 438687 x6 + 141374 x7 + 976 e -41.h0 ; 0 )
-If t = +10°C    105A = Max ( -175 x + 11332 x2 - 1318 x3 - 54120 x4 + 112625 x5 - 96545 x6 + 31284 x7 + 147 e -30.h0 ; 0 )
-If t = +15°C          A = 0
-If t = +30°C    105A = Min ( -1 + 589 x - 34750 x2 + 9753 x3 + 154745 x4 - 335229 x5 + 291742 x6 - 95395 x7 - 284 e -37.h0 ; 0 )

-The exponential terms on the right are quite arbitrary and the results are uncertain for 0° < h0 < 0°10' ( though correct for h0 = 0° and 0°10' )
-The coefficient A increases more rapidly ( for t < 15°C ) or less rapidly ( for t > 15°C ) near the horizon,
and I don't really know the behavior of the function A(h0) in this interval.
-Though slightly less accurate, the following formulae may also be used:

-If t = -30°C    105A = Max ( -15 - 387 x + 82156 x2 + 142877 x3 - 963927 x4 + 1845007 x5 - 1614635 x6 + 545974 x7 ; 0 )
-If t = -10°C    105A = Max (  -5 - 466 x + 49149 x2 + 49994 x3 - 454890 x4 + 891332 x5 - 779644 x6 + 262223 x7 ; 0 )
-If t = +10°C    105A = Max ( -1 - 113 x + 10184 x2 + 7212 x3 - 84702 x4 + 168894 x5 - 147638 x6 + 49394 x7 ; 0 )
-If t = +15°C          A = 0
-If t = +30°C    105A = Min ( +1 + 469 x - 32520 x2 - 6816 x3 + 214150 x4 - 444530 x5 + 390988 x6 - 130572 x7 ; 0 )

Coefficient B:

-If P = 500 mbar          105B =  -27 + 909 x - 42020 x2 + 102902 x3 - 101640 x4 + 16348 x5 + 39269 x6 - 19816 x7
-If P = 700 mbar          105B =  -16 + 506 x - 24962 x2 + 58265 x3 - 49889 x4 - 6869 x5 + 35957 x6 - 15541 x7
-If P = 900 mbar          105B =   -7 + 229 x - 9556 x2 + 23689 x3 - 25749 x4 + 9819 x5 + 3176 x6 - 2541 x7
-If P = 1013.25 mbar        B = 0
-If P = 1100 mbar        105B =   4 - 153 x + 7206 x2 - 18115 x3 + 21595 x4 - 12458 x5 + 2134 x6 + 572 x7

Coefficient C:

105C = [ 473 ( 0.59-Lwl ) + 1570 ( 0.59-Lwl )2 + 2911 ( 0.59-Lwl )3 ] exp ( -0.472 h00.866 )

Coefficient D:

105D =  ( -14.6 f + 2.556 f 2 - 0.12445 f 3 + f 4/214 - f 5/16540 )/( 1 + 1.057 h0 + 0.29 h02 + h03/80 )

Coefficient E:

E = ( -1/260 ) Cos ( 2.Lat ) exp ( -0.467 h00.8215 )

Coefficient F:

F = [ exp ( - alt/18031 ) - 1 ] exp ( -1.106 h00.805 )

-All these formulae may certainly be improved.
-They are empirical and I 've found many of them with my HP-48 and Sune Bredahl's excellent "SOLVESYS" library ( cf  http://www.hpcalc.org )

Data Registers:    R00 = h0                                                               ( Registers R01 thru R06 are to be initialized before executing "REFR" )

•  R01 = t     ( in °C )         ( -30 <= t <= +30 )
•  R02 = P    ( in mbar )     ( 500 <= P <= 1100 )            These limits are not absolute...
•  R03 = f     ( in mbar )      ( 0 <= f <= 30 )
•  R04 = Lwl  ( in µm )       ( 0.4 <= lwl <= 0.7 )                     R08 thru R24: temp
•  R05 = lat   ( in °. ' " )
•  R06 = alt   ( in meters )   ( 0 <= alt <= 1000 )

When the program stops,  R07 = Refraction ( in degrees and decimals )

t = temperature ; P = pressure ; f = partial pressure of water vapor
Lwl = Light wave-length ; lat = latitude
alt = observer's altitude ( over sea-level )

Flags: /
Subroutine:  "LAGR"  ( cf "Lagrange Interpolation formula for the HP-41" )

-If you want to take advantage of Laplace's formula,   replace lines 27-28 by   LBL 02
and add the following instructions after line 04:

21                 X<>Y         ENTER^        *                  1/X           GTO 02
X>Y?           TAN           X^2                CHS            +               LBL 01
GTO 01        1/X             185 E-7          63.064         *               CLX

-If you want to use the second formulae to calculate A, replace lines 77 to 102 by

15                   82156             CHS                CHS
STO 08          STO 10          STO 12            545974
387                 142877          1845007           XEQ 01
CHS                STO 11         STO 13
STO 09           963927         1614635

-Replace lines 106 to 132 by

5                     49419           454890              779644
STO 08           STO 10         CHS                  CHS
466                  575               STO 12             262223
CHS                +                   891332              XEQ 01
STO 09           STO 11         STO 13

-Replace lines 136 to 160 by

1                    10184             CHS                 CHS
STO 08          STO 10          STO 12            49394
113                7212               168894             XEQ 01
CHS              STO 11          STO 13
STO 09          84702            147638

-And replace lines 164 to 189 by

1                       32520                STO 11          STO 13
CHS                  CHS                  214150          390988
STO 08             STO 10             STO 12          130572
469                    6816                 444530           CHS
STO 09              CHS                 CHS               XEQ 01

 01 LBL "REFR"  02 DEG  03 HR  04 STO 00  05 23.82074  06 RCL Y  07 7.4078  08 +  09 /  10 +  11 8.42681  12 X<>Y  13 /  14 +  15 6.04529  16 X<>Y  17 /  18 +  19 3.81451  20 X<>Y  21 /  22 +  23 TAN  24 1/X  25 63.05561  26 /  27 X<0?  28 CLST  29 RCL 02  30 *  31 960.233  32 /  33 RCL 01  34 271.677  35 /  36 1  37 +  38 /  39 1  40 RCL 03            41 18 E4  42 /  43 6579  44 1/X  45 +  46 RCL 03  47 *  48 -  49 *  50 5  51 RCL 04  52 X^2  53 836  54 *  55 /  56 .98282  57 +  58 *  59 STO 07  60 10  61 STO 19  62 CHS  63 STO 17  64 15  65 STO 21  66 ST+ X  67 STO 23  68 CHS  69 STO 15 70  CLX   71  STO 22   72  SIGN   73  RCL 00   74  +   75  1/X   76  STO 14   77  2               78  STO 08   79  1411         80  CHS        81  STO 09   82  100967    83  STO 10   84  3583   85  STO 11   86  465432   87  CHS   88  STO 12   89  928890   90  STO 13   91  783471   92  CHS   93  251549   94  XEQ 01   95  2377   96  RCL 00          97  43   98  *   99  CHS 100  E^X 101  * 102  + 103  X<0? 104  CLX 105  STO 16 106  CLX       107  STO 08 108  880        109  CHS      110  STO 09         111  57082   112  STO 10 113  6928 114  CHS 115  STO 11 116  250807 117  CHS 118  STO 12 119  515833 120  STO 13 121  438687 122  CHS 123  141374 124  XEQ 01 125  976 126  RCL 00 127  41 128  * 129  CHS 130  E^X 131  * 132  + 133  X<0? 134  CLX 135  STO 18 136  175       137  CHS 138  STO 09 139  11332   140  STO 10 141  1318    142  CHS     143  STO 11 144  54120 145  CHS 146  STO 12 147  112625 148  STO 13 149  96545 150  CHS 151  31284 152  XEQ 01 153  147 154  RCL 00 155  30 156  * 157  CHS 158  E^X 159  * 160  + 161  X<0? 162  CLX 163  STO 20        164  1        165  STO 08 166  589          167  STO 09    168  34750  169  CHS    170  STO 10  171  9753 172  STO 11 173  154745 174  STO 12 175  335229 176  CHS 177  STO 13 178  291742 179  95395 180  CHS 181  XEQ 01 182  284 183  RCL 00         184  37 185  * 186  CHS 187  E^X 188  * 189  - 190  X>0? 191  CLX 192  STO 24 193  RCL 01 194  XEQ 02 195  500 196  STO 15 197  700 198  STO 17 199  90 200  ST* 19 201  67.55 202  ST* 21 203  1100 204  STO 23 205  27 206  STO 08 207  909 208  STO 09 209  42020 210  CHS 211  STO 10 212  102902 213  STO 11 214  101640 215  CHS 216  STO 12 217  16348 218  STO 13 219  39269 220  19816 221  CHS 222  XEQ 01 223  STO 16 224  16 225  STO 08 226  506 227  STO 09 228  24962 229  CHS 230  STO 10 231  58265 232  STO 11        233  49889 234  CHS 235  STO 12 236  6869 237  CHS 238  STO 13 239  35957 240  15541 241  CHS 242  XEQ 01 243  STO 18 244  7 245  STO 08 246  229 247  STO 09 248  9556 249  CHS 250  STO 10 251  23689 252  STO 11         253  25749 254  CHS 255  STO 12 256  9819 257  STO 13 258  3176 259  2541 260  CHS 261  XEQ 01 262  STO 20 263  4 264  CHS 265  STO 08 266  153 267  CHS 268  STO 09 269  7206 270  STO 10 271  18115 272  CHS 273  STO 11 274  21595 275  STO 12 276  12458 277  CHS 278  STO 13 279  2134 280  572 281  XEQ 01 282  STO 24 283  RCL 02 284  XEQ 02 285  RCL 03 286  RCL 03 287  RCL 03 288  16540 289  / 290  214 291  1/X 292  - 293  * 294  .12445 295  + 296  * 297  2.556 298  - 299  * 300  14.6 301  - 302  * 303  RCL 00        304  80 305  / 306  .29 307  + 308  RCL 00 309  * 310  1.057 311  + 312  RCL 00 313  * 314  1 315  + 316  / 317  XEQ 03 318  .59 319  RCL 04         320  - 321  1570 322  RCL Y 323  2911 324  * 325  + 326  * 327  473 328  + 329  * 330  RCL 00 331  .866 332  Y^X 333  .472 334  * 335  E^X 336  / 337  XEQ 03 338  GTO 04 339  LBL 01 340  RCL 14 341  STO T 342  * 343  + 344  * 345  RCL 13 346  + 347  * 348  RCL 12 349  + 350  * 351  RCL 11 352  + 353  * 354  RCL 10 355  + 356  * 357  RCL 09 358  + 359  * 360  RCL 08 361  - 362  RTN 363  LBL 02 364  15.024 365  X<>Y 366  XEQ "LAGR" 367  LBL 03 368   E5 369  ST+ Y 370  / 371  ST* 07 372  RTN 373  LBL 04 374  1 375  RCL 05 376  HR 377  ST+ X 378  COS 379  260 380  / 381  RCL 00         382  .8215 383  Y^X 384  .467 385  * 386  E^X 387  / 388  - 389  ST* 07 390  RCL 06 391  18031 392  / 393  CHS 394  E^X-1 395  RCL 00 396  .805 397  Y^X 398  1.106 399  * 400  E^X 401  / 402  1 403  + 404  ST* 07 405  RCL 00 406  RCL 07 407  ST- Y 408  HMS 409  X<>Y 410  HMS 411  END

( 852 bytes / SIZE 025 )

 STACK INPUTS OUTPUTS Y / R X h0 h

( All angles expressed in ° ' " )   Execution time = 64 seconds.

Example:   With  t = 20°C , P = 1000 mbar , f = 12 mbar , Lwl = 0.5 µm , lat = 30° , alt = 500 m           ( store these 6 numbers into R01 thru R06 )

h0 = 0               XEQ "REFR"  >>>>   h =   -0°30'03"88    X<>Y    R =  0°30'03"88
h0 = 1°                    R/S           >>>>   h =    0°37'43"50    X<>Y    R =  0°22'16"50
h0 = 12°34'56"        R/S           >>>>   h =  12°30'52"86    X<>Y    R =  0°04'03"14
h0 = 41°16'24"        R/S           >>>>   h =  41°15'20"85    X<>Y    R =  0°01'03"15

-The accuracy is of the order of 1 or 2 arcseconds near the horizon ( with a greater uncertainty for 0° < h0 < 0°10' ). Interpolation may also decrease the precision.
- errors ~  0"5 for h0 = 5°
- errors ~  0"2 for h0 = 10°
-For h0 > 20° the accuracy is determined by the formula which computes R0 , therefore use Laplace's formula
as explained on the right of the beginning of this listing if you need more accurate results for these h0-values.

-As mentioned above, the accurate results near the horizon are often unrealistic without knowing the detailed structure of low atmosphere.
( cf references [3] & [9] for a detailed analysis on low-altitude refraction )

References:

[1]  Jean Kovalevsky et al. - "Introduction aux Ephemerides Astronomiques" - EDP Sciences - ISBN 2-86883-298-9  ( in French )
[2]  Abalakin 1985 - "Refraction Tables of Pulkovo Observatory" 5th edition - Nauka, Leningrad
[3]  Andrew T. Young - "Sunset Science IV. Low-Altitude Refraction" - 2004, the Astronomical Journal, 127 , 3622
[4]  Lawrence H. Auer & E. Myles Standish - "Astronomical Refraction, Computational method for all zenith angles" 2000 AJ, 119 , 2472
[5]  Krystyna Kurzynska - "Precision in determination of astronomical refraction from aerological data" - 1987, Astron. Nachr. 308, 323
[6]  Krystyna Kurzunska - "Local effects in pure astronomical refraction" - 1988, Astron. Nachr. 309, 57
[7]  Krystyna Kurzynska - "On the accuracy of the Refraction Tables of Pulkovo Observatory, 5th edition" 1988, Astron. Nachr. 309, 213
[8]  Minodora Lipcanu - "A direct method for the calculation of astronomical refraction"
[9]  Andrew T. Young - "Understanding Astronomical Refraction" - 2006, The Observatory, Vol. 126, N° 1191, pp. 82-115