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 h₀ and the true altitude h of a given star:   h = h₀ - 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( h₀ + 4.2206/( h₀ + 15.1115/( h₀ + 5.9431 ) ) )       where h₀ 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  h₀ = 1°30'00"  ;   h₀ = 27°  ;   h₀ = 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 h₀ - 0".0666/Tan³ h₀  gives even better results if h₀ > 20° , more precisely:

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

     b) Other Simple Formulae
 

-We can build similar formulae for other atmospheric conditions.

         R  ~   a/Tan( h₀ + b/( h₀ + c/( h₀ + d ) ) )       where h₀ 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(h₀)  for  0° < h₀ < 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"  >>>>  h₀ = 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( h₀ + 4.2206/( h₀ + 15.1115/( h₀ + 5.9431 ) ) )       where h₀ 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   h₀ =  12°34'56"    12.3456  XEQ "REF"  >>>>   h = 12°29'58"4   X<>Y   R = 0°04'57"6
-If   h₀ =  10°12'34"    10.1234       R/S          >>>>   h = 10°06'29"5   X<>Y   R = 0°06'04"5

-The results are acceptable for  h₀ > 10°  ( errors are of the order of 1 or 2 arcseconds for h₀ = 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:

      h₀ - 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"  >>>     h₀ = 20°0495
  if h = 10°     10      R/S           >>>     h₀ = 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  R₀ ~ (1°/62.8093)/Tan( h₀ + 4.2206/( h₀ + 15.1115/( h₀ + 5.9431 ) ) )

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

Data Registers:    R00 = h₀                                                                 ( 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   h₀ =  0°                     0      XEQ "REF2"  >>>>   h = -0°33'32"       X<>Y    R = 0°33'32"
-If   h₀ =  10°23'45"    10.2345       R/S          >>>>   h = 10°19'02"8    X<>Y    R = 0°04'42"2
-If   h₀ =  49°12'34"    49.1234       R/S          >>>>   h = 49°11'47"9    X<>Y    R = 0°00'46"1

-If  h₀ > 20° more accurate results will be obtained if  R₀ is computed by Laplace's formula.
-The accuracy is of the order of 10" near the horizon, but errors rapidly decrease as h₀ increases.
-However, disturbances of the atmosphere make accurate results more theoretical than real for very small h₀-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"  >>>>    h₀ = 1°20'18"7
  h = 10°23'45"     10.2345        R/S          >>>>    h₀ = 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:

   R₀ ~ (1°/63.05561)/Tan( h₀ + 3.81451/( h₀ + 6.04529/( h₀ + 8.42681/( h₀ + 23.82074/( h₀ + 7.40780 ) ) ) ) )

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

  R = R₀ (1.0552126/(1+0.00368084.t)).(1+A).(P/1013.25).(1+B).(0.98282+0.005981/Lwl²).(1+C).(1-0.152 10 ^-3 f -0.55 10 ^-5 f ² ).(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+h₀)

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

-The exponential terms on the right are quite arbitrary and the results are uncertain for 0° < h₀ < 0°10' ( though correct for h₀ = 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(h₀) in this interval.
-Though slightly less accurate, the following formulae may also be used:

-If t = -30°C    10⁵A = Max ( -15 - 387 x + 82156 x² + 142877 x³ - 963927 x⁴ + 1845007 x⁵ - 1614635 x⁶ + 545974 x⁷ ; 0 )
-If t = -10°C    10⁵A = Max (  -5 - 466 x + 49149 x² + 49994 x³ - 454890 x⁴ + 891332 x⁵ - 779644 x⁶ + 262223 x⁷ ; 0 )
-If t = +10°C    10⁵A = Max ( -1 - 113 x + 10184 x² + 7212 x³ - 84702 x⁴ + 168894 x⁵ - 147638 x⁶ + 49394 x⁷ ; 0 )
-If t = +15°C          A = 0
-If t = +30°C    10⁵A = Min ( +1 + 469 x - 32520 x² - 6816 x³ + 214150 x⁴ - 444530 x⁵ + 390988 x⁶ - 130572 x⁷ ; 0 )
 

Coefficient B:

-If P = 500 mbar          10⁵B =  -27 + 909 x - 42020 x² + 102902 x³ - 101640 x⁴ + 16348 x⁵ + 39269 x⁶ - 19816 x⁷
-If P = 700 mbar          10⁵B =  -16 + 506 x - 24962 x² + 58265 x³ - 49889 x⁴ - 6869 x⁵ + 35957 x⁶ - 15541 x⁷
-If P = 900 mbar          10⁵B =   -7 + 229 x - 9556 x² + 23689 x³ - 25749 x⁴ + 9819 x⁵ + 3176 x⁶ - 2541 x⁷
-If P = 1013.25 mbar        B = 0
-If P = 1100 mbar        10⁵B =   4 - 153 x + 7206 x² - 18115 x³ + 21595 x⁴ - 12458 x⁵ + 2134 x⁶ + 572 x⁷

Coefficient C:

  10⁵C = [ 473 ( 0.59-Lwl ) + 1570 ( 0.59-Lwl )² + 2911 ( 0.59-Lwl )³ ] exp ( -0.472 h₀^0.866 )

Coefficient D:

  10⁵D =  ( -14.6 f + 2.556 f ² - 0.12445 f ³ + f ⁴/214 - f ⁵/16540 )/( 1 + 1.057 h₀ + 0.29 h₀² + h₀³/80 )

Coefficient E:

       E = ( -1/260 ) Cos ( 2.Lat ) exp ( -0.467 h₀^0.8215 )

Coefficient F:

       F = [ exp ( - alt/18031 ) - 1 ] exp ( -1.106 h₀^0.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 = h₀                                                               ( 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 )

    h₀ = 0               XEQ "REFR"  >>>>   h =   -0°30'03"88    X<>Y    R =  0°30'03"88  
    h₀ = 1°                    R/S           >>>>   h =    0°37'43"50    X<>Y    R =  0°22'16"50
    h₀ = 12°34'56"        R/S           >>>>   h =  12°30'52"86    X<>Y    R =  0°04'03"14
    h₀ = 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° < h₀ < 0°10' ). Interpolation may also decrease the precision.
- errors ~  0"5 for h₀ = 5°
- errors ~  0"2 for h₀ = 10°
-For h₀ > 20° the accuracy is determined by the formula which computes R₀ , therefore use Laplace's formula
 as explained on the right of the beginning of this listing if you need more accurate results for these h₀-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