Refraction integral

On a ray tracing method to calculate atmospheric refraction

Aimo Niemi, 2009

Refraction is the bending of a light ray, when it travels in the earth atmosphere. If the ray is coming from the space, the bending is called astronomical refraction. If both ends of the ray are inside the atmosphere, the phenomenon is called terrestrial refraction. The reason for the bending is the increasing density of the air and its increasing refractive index, when the ray penetrates deeper and deeper into the atmosphere. If now the refractive index n is known throughout the path of the ray, that bending is easily calculated. See the picture below, click to enlarge it.

Fig. 1

Referring to fig.1 the observer at point A_o, on height y above the ground level, sees a ray coming from an apparent zenith distance z_o. In the picture the atmosphere is thought to be divided into concentric layers with the heights dy_i, where the refractive index is supposed to have a constant value n_i. At the point A_i on the border of the layers (i-1) and i light is refracted according to Snell's law so that

If further r_i is the distance between point A_i and the earth center C, then, using the sine rule in the triangle A_i-1CA_i, we get

from which

Starting from a value f=0 and using some values of our choice for dy we can calculate f=f+f_i-1 and r_i = r_i-1+dy_i-1 which enables us to move from one layer to the next and continue the calculation there the same way. If terrestrial refraction is calculated, we can stop, when f or r_i has reached a predefined value. If we calculate astronomical refraction, we may stop, when the height of the light ray is something over 80 km. True zenith distance of the ray is then z_i+f, where i refers to the last layer, which was calculated and the refraction, of course, is then z_i+f-z_o (If all this is too simple, try this ;)

Refractive index can be calculated using formula

where p is pressure in kPa, t is temperature in Celsius, and RH is relative humidity in percent. The formula is given at NIST where is also a more precise online tool, which calculates the index of refraction of air and wavelength of light in air as a function of various input parameters. Usually, however, p and t are known only at the ground level and we must use some standard theory to calculate how they change in the course of the light ray. For this purpose The US Standard Atmosphere 1976 is a good choice. In it the atmosphere is divided into several regions where inside the temperature and pressure are defined with simple formulas. For instance, its lowest part, Troposphere, extends from the ground level to the height of 11 km at which the temperature changes linearly from +15°C to -56.5°C and the pressure exponentially from 1013.25 hPa to 223.32 hPa.

Pc-Calculator program
After the previous short introduction we can now create programs to calculate the refraction. The example below was created with Pc Calculator, a clever calculator program combined with a smart text and formula editor. With Pc Calculator anybody can make simple calculations but here we use also its advanced properties, which are shortly discussed below.

The working memory of the calculator is organised as lines, which the user can execute one by one as she or he wills.

The cursor can be moved onto any line and when ENTER is pressed, the result of the line is put to the adjacent answer window.
Each line has also a unique line number, which calculators line function Ł (pound sign) can use as an argument to execute that line. Except numbers, the argument can also be an expression or a register name, where some constant values are stored.
Each line may also have a name, which the quote function " can use as an argument to execute that line.

Look now the example below and find there the line 14, which ends as :tp.y That is how line's name is defined. If you now want to execute that line, you can do so by executing "tp.y. This would work fine but, if speed is important that is not the best way because each time the formula interpreter finds such a line reference, it must search the whole memory to find the right line. A faster way is used in the chapter "Initiate coefficients" and there in the expression A="tp.y#. Here A gets the real line number of line tp.y, which can now be executed faster with ŁA. After the execution the line (or quote) function returns also the value of the line, which it just executed. If there are many expressions on the line, the value of the last one is returned.

Now look again at the line tp.y and find there the function #. It returns here the real line number, where it is situated and, if we now execute Ł(#-1), the line immediately before is executed. For instance, if y on the line tp.y has a value 33, then the logical expressions 11<y, 20<y, 32<y all have the value -1 while the others are zeroes and the line function there becomes equivalent to Ł(#-4) which will execute the line, where Y gets the value 32. Note that on that line there is also an expression 228.65>@. In Pc Calculator this is not a logical operation. Instead here the value 228.65 is stored to the temporary register @ where it is used later when pressure P is calculated. Note also that here the capital letter Y refers to a different memory place (or a register) than the small letter y.

Main part of the calculation is performed on the first line in the chapter "Calculate refraction". Readers who are not used Pc Calculator may find it difficult to understand so, it is rewritten below and each expression on the line is put there on a line of its own.
    G=D/(P+30)
    S=R+G
    y=y+G
    x=R*sin(z)/S
    f=f+z-arcsin(x)
    R=S
    z=arcsin[N*x/(ŁB)]
This line is then repeated in the loop, which the next line, "trace.refraction", creates. The loop ends, when y has become greater than 85 km. If you calculate terrestrial refraction or if you, for instance, want to know how far away a refraction raised horizon is, you can change the loop end test here. To test the program load Pc Calculator and use copy and paste.


'PC CALCULATOR, Aimo Niemi 2009
'**********************************************
'Standard Atmosphere T(y) P(y)
'**********************************************
K=-9806.65/287.05307
Y=84.852 X=0 T=1 P=0
Y=71 X=-2 T=214.65>@+X(y-Y) P=0.0396(T/@)^(K/X)
Y=51 X=-2.8 T=270.65>@+X(y-Y) P=0.6694(T/@)^(K/X)
Y=47 X=0 T=270.65 P=1.1091e(K(y-Y)/T)
Y=32 X=2.8 T=228.65>@+X(y-Y) P=8.6802(T/@)^(K/X)
Y=20 X=1 T=216.65>@+X(y-Y) P=54.7489(T/@)^(K/X)
Y=11 X=0 T=216.65 P=226.3206e(K(y-Y)/T)
Y=0 X=-6.5 T=288.15>@+X(y-Y) P=1013.25(T/@)^(K/X)
Ł(#-1+11<y+20<y+32<y+47<y+51<y+71<y+84<y):tp.y

'**********************************************
'Refractive index N(T,P,Q=RH)
'**********************************************
M= 7.839/10^5 :700nm 
M= 7.860/10^5 :633nm 
M= 7.897/10^5 :550nm 
M= 7.930/10^5 :500nm
ŁA N=1+MP/T+(y<11)15((T-273)^2+160)Q/10^12 :n.qth

'**********************************************
'Initiate coefficients, adjust R to your local latitude
'**********************************************
R=6371+y K=-9806.65/287.05307 M="550 f=0 D=1-z<80
Ł(#-1) A="tp.y# B="n.qth# @=ŁA N=ŁB @="trace :startref

'**********************************************
'Calculate refraction G=dy z=z(i) R=r(i) S=r(i+1)
'**********************************************
G=D/(P+30) S=R+G y=y+G x=Rsz/S f=f+z-asx R=S z=as(Nx/ŁB)
Ł(#-1)Ł(#-85<y) :trace.refraction
f+z+( ?g' beep and print the result

'f+z+( ?g'= 79°41'01.6502" M="550nm Q=0%
'f+z+( ?g'= 79°41'01.0704" M="550nm Q=100%

'**********************************************
'give Q=RH, y=ground height, z=apparent zenith distance and
'execute "startref to get the true zenith distance = (f+z)
'**********************************************
Q=0 y=0 z=79.6 Z=z @="startref