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 (i1) 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_{i1}CA_{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_{i1} and
r_{i} = r_{i1}+dy_{i1}
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}+fz_{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.
PcCalculator 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+zarcsin(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(yY) P=0.0396(T/@)^(K/X)
Y=51 X=2.8 T=270.65>@+X(yY) P=0.6694(T/@)^(K/X)
Y=47 X=0 T=270.65 P=1.1091e(K(yY)/T)
Y=32 X=2.8 T=228.65>@+X(yY) P=8.6802(T/@)^(K/X)
Y=20 X=1 T=216.65>@+X(yY) P=54.7489(T/@)^(K/X)
Y=11 X=0 T=216.65 P=226.3206e(K(yY)/T)
Y=0 X=6.5 T=288.15>@+X(yY) 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((T273)^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=1z<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+zasx 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


