Air mass

In astronomy air mass is a measure of the amount of atmosphere between the observer and the object being observed. It is a function of the zenith angle (the angle between the line of sight an vertical) and the altitude of the observer. It is defined as the integral of the atmospheric density along the line of sight and is usually expressed relative to the air mass at zenith. Thus, looking straight up gives an air mass of one (regardless of observer's altitude) and viewing at any zenith angle greater than zero gives higher values.

You will need to integrate $\rho$ (h(a,z,x)) where $\rho$ (h) is the atmospheric density for a given height above sea level, and h(a,z,x) is the height above sea level for a point at distance x along the line of sight. Determining this last function requires some trigonometry.

For this task you can assume:

The density of Earth's atmosphere is proportional to exp(-a/8500 metres)
The Earth is a perfect sphere of radius 6731 km.

Task

Write a function that calculates the air mass for an observer at a given altitude a above sea level and zenith angle z.
Show the air mass for zenith angles 0 to 90 in steps of 5 degrees for an observer at sea level.
Do the same for the NASA SOFIA infrared telescope, which has a cruising altitude of 13,700 meters (about 8.3 miles),

it flies in a specially retrofitted Boeing 747 about four flights a week.

AWK

syntax: GAWK -f AIR_MASS.AWK
converted from FreeBASIC

BEGIN {

   dd = 0.001  # integrate in this fraction of the distance already covered
   DEG = 0.017453292519943295769236907684886127134 # degrees to radians
   RE = 6371000 # Earth radius in meters
   print("Angle          0 m      13700 m")
   for (z=0; z<=90; z+=5) {
     printf("%5d %12.8f %12.8f\n",z,am_airmass(0,z),am_airmass(13700,z))
   }
   exit(0)

} function am_airmass(a,z) {

   return am_column_density(a,z) / am_column_density(a,0)

} function am_column_density(a,z, d,delta,sum) { # integrates density along the line of sight

   while (d < 10000000) { # integrate only to a height of 10000km, effectively infinity
     delta = max(dd,(dd)*d) # adaptive step size to avoid it taking forever
     sum += am_rho(am_height(a,z,d+0.5*delta))*delta
     d += delta
   }
   return(sum)

} function am_height(a,z,d, aa,hh) {

a - altitude of observer
z - zenith angle in degrees
d - distance along line of sight

   aa = RE + a
   hh = sqrt(aa^2 + d^2 - 2*d*aa*cos((180-z)*DEG))
   return(hh-RE)

} function am_rho(a) { # density of air as a function of height above sea level

   return exp(-a/8500.0)

} function max(x,y) { return((x > y) ? x : y) } </lang>

Output:

Angle          0 m      13700 m
    0   1.00000000   1.00000000
    5   1.00380963   1.00380965
   10   1.01538466   1.01538475
   15   1.03517744   1.03517765
   20   1.06399053   1.06399093
   25   1.10305937   1.10306005
   30   1.15418974   1.15419083
   35   1.21998076   1.21998246
   40   1.30418931   1.30419190
   45   1.41234169   1.41234567
   50   1.55280404   1.55281025
   55   1.73875921   1.73876915
   60   1.99212000   1.99213665
   65   2.35199740   2.35202722
   70   2.89531368   2.89537287
   75   3.79582352   3.79596149
   80   5.53885809   5.53928113
   85  10.07896219  10.08115981
   90  34.32981136  34.36666557

Factor

Translation of: FreeBASIC

Works with: Factor version 0.99 2021-02-05

<lang factor>USING: formatting io kernel math math.functions math.order math.ranges math.trig sequences ;

CONSTANT: RE 6,371,000 ! Earth's radius in meters CONSTANT: dd 0.001 ! integrate in this fraction of the distance already covered CONSTANT: FIN 10,000,000 ! integrate to a height of 10000km

! the density of air as a function of height above sea level

rho ( a -- x ) neg 8500 / e^ ;

! z = zenith angle (in degrees) ! d = distance along line of sight ! a = altitude of observer

height ( a z d -- x )

   RE a + :> AA
   AA sq d sq + 180 z - deg>rad cos AA * d * 2 * - sqrt RE - ;

column-density ( a z -- x )

   ! integrates along the line of sight
   0 0 :> ( s! d! )
   [ d FIN < ] [
       dd dd d * max :> delta   ! adaptive step size to avoid taking it forever
       s a z d 0.5 delta * + height rho delta * + s!
       d delta + d!
   ] while s ;

airmass ( a z -- x )

   [ column-density ] [ drop 0 column-density ] 2bi / ;

"Angle 0 m 13700 m" print "------------------------------------" print 0 90 5 <range> [

   dup [ 0 swap airmass ] [ 13700 swap airmass ] bi
   "%2d %15.8f %17.8f\n" printf

] each</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0      1.00000000        1.00000000
 5      1.00380963        1.00380965
10      1.01538466        1.01538475
15      1.03517744        1.03517765
20      1.06399053        1.06399093
25      1.10305937        1.10306005
30      1.15418974        1.15419083
35      1.21998076        1.21998246
40      1.30418931        1.30419190
45      1.41234169        1.41234567
50      1.55280404        1.55281025
55      1.73875921        1.73876915
60      1.99212000        1.99213665
65      2.35199740        2.35202722
70      2.89531368        2.89537287
75      3.79582352        3.79596149
80      5.53885809        5.53928113
85     10.07896219       10.08115981
90     34.32981136       34.36666557

FreeBASIC

define DEG 0.017453292519943295769236907684886127134 'degrees to radians
define RE 6371000 'Earth radius in meters
define dd 0.001 'integrate in this fraction of the distance already covered
define FIN 10000000 'integrate only to a height of 10000km, effectively infinity
define max(a, b) iif(a>b,a,b)

function rho(a as double) as double

   'the density of air as a function of height above sea level
   return exp(-a/8500.0)

end function

function height( a as double, z as double, d as double ) as double

   'a = altitude of observer
   'z = zenith angle (in degrees)
   'd = distance along line of sight
   dim as double AA = RE + a, HH
   HH = sqr( AA^2 + d^2 - 2*d*AA*cos((180-z)*DEG) )
   return HH - RE

end function

function column_density( a as double, z as double ) as double

   'integrates density along the line of sight
   dim as double sum = 0.0, d = 0.0, delta
   while d<FIN
       delta = max(dd, (dd)*d)  'adaptive step size to avoid it taking forever:
       sum += rho(height(a, z, d+0.5*delta))*delta
       d += delta
   wend
   return sum

end function

function airmass( a as double, z as double ) as double

   return column_density( a, z ) / column_density( a, 0 )

end function

print "Angle 0 m 13700 m" print "------------------------------------" for z as double = 0 to 90 step 5.0

   print using "##      ##.########      ##.########";z;airmass(0, z);airmass(13700, z)

next z </lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

Julia

Translation of: FreeBASIC

<lang julia>using Printf

const DEG = 0.017453292519943295769236907684886127134 # degrees to radians const RE = 6371000 # Earth radius in meters const dd = 0.001 # integrate in this fraction of the distance already covered const FIN = 10000000 # integrate only to a height of 10000km, effectively infinity

""" the density of air as a function of height above sea level """ rho(a::Float64)::Float64 = exp(-a / 8500.0)

""" a = altitude of observer

   z = zenith angle (in degrees)
   d = distance along line of sight """

height(a, z, d) = sqrt((RE + a)^2 + d^2 - 2 * d * (RE + a) * cosd(180 - z)) - RE

""" integrates density along the line of sight """ function column_density(a, z)

   dsum, d = 0.0, 0.0
   while d < FIN
       delta = max(dd, (dd)*d)  # adaptive step size to avoid it taking forever:
       dsum += rho(height(a, z, d + 0.5 * delta)) * delta
       d += delta
   end
   return dsum

end

airmass(a, z) = column_density(a, z) / column_density(a, 0)

println("Angle 0 m 13700 m\n", "-"^36) for z in 0:5:90

   @printf("%2d      %11.8f      %11.8f\n", z, airmass(0, z), airmass(13700, z))

end

</lang>

Output:

Angle           0 m          13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

Nim

Translation of: Wren

<lang Nim>import math, strformat

const

 Re = 6371000  # Radius of earth in meters.
 Dd= 0.001     # Integrate in this fraction of the distance already covered.
 Fin = 1e7     # Integrate only to a height of 10000km, effectively infinity.

func rho(a: float): float =

 ## The density of air as a function of height above sea level.
 exp(-a / 8500)

func height(a, z, d: float): float =

 ## Height as a function of altitude (a), zenith angle (z)
 ## in degrees and distance along line of sight (d).
 let aa = Re + a
 let hh = sqrt(aa * aa + d * d - 2 * d * aa * cos(degToRad(180-z)))
 result = hh - Re

func columnDensity(a, z: float): float =

 ## Integrates density along the line of sight.
 var d = 0.0
 while d < Fin:
   let delta = max(Dd, Dd * d)   # Adaptive step size to avoid it taking forever.
   result += rho(height(a, z, d + 0.5 * delta)) * delta
   d += delta

func airmass(a, z: float): float =

 columnDensity(a, z) / columnDensity(a, 0)

echo "Angle 0 m 13700 m" echo "------------------------------------" var z = 0.0 while z <= 90:

 echo &"{z:2}      {airmass(0, z):11.8f}      {airmass(13700, z):11.8f}"
 z += 5</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

Perl

Translation of: Raku

<lang perl>use strict; use warnings; use feature <say signatures>; no warnings 'experimental::signatures'; use List::Util 'max';

use constant PI => 2*atan2(1,0); # π use constant DEG => PI/180; # degrees to radians use constant RE => 6371000; # Earth radius in meters use constant dd => 0.001; # integrate in this fraction of the distance already covered use constant FIN => 10000000; # integrate only to a height of 10000km, effectively infinity

Density of air as a function of height above sea level

sub rho ( $a ) {

   exp( -$a / 8500 );

}

sub height ( $a, $z, $d ) {

   # a = altitude of observer
   # z = zenith angle (in degrees)
   # d = distance along line of sight
   my $AA = RE + $a;
   my $HH = sqrt $AA**2 + $d**2 - 2 * $d * $AA * cos( (180-$z)*DEG );
   $HH - RE;

}

Integrates density along the line of sight

sub column_density ( $a, $z ) {

   my $sum = 0;
   my $d   = 0;
   while ($d < FIN) {
       my $delta = max(dd, dd * $d);  # Adaptive step size to avoid it taking forever
       $sum += rho(height($a, $z, $d + $delta/2))*$delta;
       $d   += $delta;
   }
   $sum;

}

sub airmass ( $a, $z ) {

   column_density($a, $z) / column_density($a, 0);

}

say 'Angle 0 m 13700 m'; say '------------------------------------'; for my $z (map{ 5*$_ } 0..18) {

   printf "%2d      %11.8f      %11.8f\n", $z, airmass(0, $z), airmass(13700, $z);

}</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

Phix

constant RE = 6371000,  // radius of earth in meters
         DD = 0.001,    // integrate in this fraction of the distance already covered
         FIN = 1e7      // integrate only to a height of 10000km, effectively infinity
 
// The density of air as a function of height above sea level.
function rho(atom a) return exp(-a/8500) end function
 
// a = altitude of observer
// z = zenith angle (in degrees)
// d = distance along line of sight
function height(atom a, z, d)
    atom aa = RE + a,
         hh = sqrt(aa*aa + d*d - 2*d*aa*cos((180-z)*PI/180))
    return hh - RE
end function
 
// Integrates density along the line of sight.
function columnDensity(atom a, z)
    atom res = 0,
         d = 0
    while d<FIN do
        atom delta = max(DD, DD*d) // adaptive step size to avoid it taking forever
        res += rho(height(a, z, d + 0.5*delta))*delta
        d += delta
    end while
    return res
end function
 
function airmass(atom a, z) return columnDensity(a,z)/columnDensity(a,0) end function
 
printf(1,"Angle     0 m              13700 m\n")
printf(1,"------------------------------------\n")
for z=0 to 90 by 5 do
    printf(1,"%2d      %11.8f      %11.8f\n", {z, airmass(0,z), airmass(13700,z)})
end for

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

Raku

Translation of: FreeBASIC

<lang perl6>constant DEG = pi/180; # degrees to radians constant RE = 6371000; # Earth radius in meters constant dd = 0.001; # integrate in this fraction of the distance already covered constant FIN = 10000000; # integrate only to a height of 10000km, effectively infinity

| Density of air as a function of height above sea level

sub rho ( \a ) { exp( -a / 8500 ) }

sub height ( \a, \z, \d ) {

   # a = altitude of observer
   # z = zenith angle (in degrees)
   # d = distance along line of sight
   my \AA = RE + a;
   sqrt( AA² + d² - 2*d*AA*cos((180-z)*DEG) ) - AA;

}

| Integrates density along the line of sight

sub column_density ( \a, \z ) {

   my $sum = 0;
   my $d   = 0;
   while $d < FIN {
       my \delta = max(dd, (dd)*$d);  # Adaptive step size to avoid it taking forever
       $sum += rho(height(a, z, $d + delta/2))*delta;
       $d   += delta;
   }
   $sum;

}

sub airmass ( \a, \z ) {

   column_density( a, z ) /
   column_density( a, 0 )

}

say 'Angle 0 m 13700 m'; say '------------------------------------'; say join "\n", (0, 5 ... 90).hyper(:3batch).map: -> \z {

   sprintf "%2d      %11.8f      %11.8f", z, airmass(    0, z), airmass(13700, z);

}</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

REXX

Translation of: FreeBASIC

<lang rexx>/*REXX pgm calculates the air mass above an observer and an object for various angles.*/ numeric digits (length(pi()) - length(.)) % 4 /*calculate the number of digits to use*/ parse arg aLO aHI aBY oHT . /*obtain optional arguments from the CL*/ if aLO== | aLO=="," then aLO= 0 /*not specified? Then use the default.*/ if aHI== | aHI=="," then aHI= 90 /* " " " " " " */ if aBY== | aBY=="," then aBY= 5 /* " " " " " " */ if oHT== | oHT=="," then oHT= 13700 /* " " " " " " */ w= 30; @ama= 'air mass at' /*column width for the two air_masses. */ say 'angle|'center(@ama "sea level", w) center(@ama commas(oHT) 'meters', w) /*title*/ say "─────┼"copies(center(, w, "─"), 2)'─' /*display the title sep for the output.*/ y= left(, w-20) /*Y: for alignment of the output cols.*/

     do j=aLO  to aHI  by aBY;        am0= airM(0, j);                 amht= airM(oHT, j)
     say center(j, 5)'│'right( format(am0, , 8), w-10)y  right( format(amht, , 8), w-10)y
     end   /*j*/

say "─────┴"copies(center(, w, "─"), 2)'─' /*display the foot separator for output*/ exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ airM: procedure; parse arg a,z; if z==0 then return 1; return colD(a, z) / colD(a, 0) d2r: return r2r( arg(1) * pi() / 180) /*convert degrees ──► radians. */ pi: pi= 3.1415926535897932384626433832795028841971693993751058209749445923078; return pi rho: procedure; parse arg a; return exp(-a / 8500) r2r: return arg(1) // (pi() * 2) /*normalize radians ──► a unit circle. */ e: e= 2.718281828459045235360287471352662497757247093699959574966967627724; return e commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? /*──────────────────────────────────────────────────────────────────────────────────────*/ cos: procedure; parse arg x; x= r2r(x); a= abs(x); numeric fuzz min(6, digits() - 3)

     hpi= pi*.5;  if a=pi    then return -1;   if a=hpi | a=hpi*3  then return 0;    z= 1
                  if a=pi/3  then return .5;   if a=pi*2/3         then return -.5;  _= 1
     x= x*x;  p= z;      do k=2  by 2;  _= -_ * x / (k*(k-1));     z= z + _
                         if z=p  then leave;   p= z;   end;                    return z

/*──────────────────────────────────────────────────────────────────────────────────────*/ exp: procedure; parse arg x; ix= x%1; if abs(x-ix)>.5 then ix= ix + sign(x); x= x-ix

     z=1;  _=1;   w=z;     do j=1; _= _*x/j;  z=(z+_)/1;  if z==w  then leave;  w=z;  end
     if z\==0  then z= z * e() ** ix;                                          return z/1

/*──────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d= digits(); numeric digits; h= d+6

     numeric form; parse value format(x,2,1,,0) 'E0'  with  g 'E' _ .;  g= g * .5'e'_ % 2
     m.=9;     do j=0  while h>9;       m.j= h;               h= h%2 + 1;      end  /*j*/
               do k=j+5  to 0  by -1;   numeric digits m.k;   g= (g+x/g)*.5;   end  /*k*/
     numeric digits d;                  return g/1

/*──────────────────────────────────────────────────────────────────────────────────────*/ elev: procedure; parse arg a,z,d; earthRad= 6371000 /*earth radius in meters.*/

     aa= earthRad + a;  return sqrt(aa**2 + d**2 - 2*d*aa*cos( d2r(180-z) ) )  - earthRad

/*──────────────────────────────────────────────────────────────────────────────────────*/ colD: procedure; parse arg a,z; sum= 0; d= 0; dd= .001; infinity= 10000000

                  do while d<infinity;  delta= max(dd, dd*d)
                  sum= sum  +  rho( elev(a, z, d + 0.5*delta) ) * delta;     d= d + delta
                  end   /*while*/
     return sum</lang>

output when using the default inputs:

angle|    air mass at sea level        air mass at 13,700 meters
─────┼─────────────────────────────────────────────────────────────
  0  │          1.00000000                     1.00000000
  5  │          1.00380963                     1.00380965
 10  │          1.01538466                     1.01538475
 15  │          1.03517744                     1.03517765
 20  │          1.06399053                     1.06399093
 25  │          1.10305937                     1.10306005
 30  │          1.15418974                     1.15419083
 35  │          1.21998076                     1.21998246
 40  │          1.30418931                     1.30419190
 45  │          1.41234169                     1.41234567
 50  │          1.55280404                     1.55281025
 55  │          1.73875921                     1.73876915
 60  │          1.99212000                     1.99213665
 65  │          2.35199740                     2.35202722
 70  │          2.89531368                     2.89537287
 75  │          3.79582352                     3.79596149
 80  │          5.53885809                     5.53928113
 85  │         10.07896219                    10.08115981
 90  │         34.32981136                    34.36666557
─────┴─────────────────────────────────────────────────────────────

Wren

Translation of: FreeBASIC

Library: Wren-math

Library: Wren-fmt

<lang ecmascript>import "/math" for Math import "/fmt" for Fmt

// constants var RE = 6371000 // radius of earth in meters var DD = 0.001 // integrate in this fraction of the distance already covered var FIN = 1e7 // integrate only to a height of 10000km, effectively infinity

// The density of air as a function of height above sea level. var rho = Fn.new { |a| Math.exp(-a/8500) }

// a = altitude of observer // z = zenith angle (in degrees) // d = distance along line of sight var height = Fn.new { |a, z, d|

   var aa = RE + a
   var hh = (aa * aa + d * d - 2 * d * aa * (Math.radians(180-z).cos)).sqrt
   return hh - RE

}

// Integrates density along the line of sight. var columnDensity = Fn.new { |a, z|

   var sum = 0
   var d = 0
   while (d < FIN) {
       var delta = Math.max(DD, DD * d) // adaptive step size to avoid it taking forever
       sum = sum + rho.call(height.call(a, z, d + 0.5 * delta)) * delta
       d = d + delta
   }
   return sum

}

var airmass = Fn.new { |a, z| columnDensity.call(a, z) / columnDensity.call(a, 0) }

System.print("Angle 0 m 13700 m") System.print("------------------------------------") var z = 0 while (z <= 90) {

   Fmt.print("$2d      $11.8f      $11.8f", z, airmass.call(0, z), airmass.call(13700, z))
   z = z + 5

}</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557

XPL0

Translation of: FreeBASIC

<lang XPL0>define DEG = 0.017453292519943295769236907684886127134; \degrees to radians define RE = 6371000.; \Earth radius in meters define DD = 0.001; \integrate in this fraction of the distance already covered define FIN = 10000000.; \integrate only to a height of 10000km, effectively infinity

function real Max(A, B); real A, B;

   return (if A>B then A else B);

function real Rho(A); real A; [ \the density of air as a function of height above sea level

   return Exp(-A/8500.0)

end; \function

function real Height( A, Z, D ); real A, \= altitude of observer

    Z, \= zenith angle (in degrees)
    D; \= distance along line of sight

real AA, HH; [ AA:= RE + A;

   HH:= sqrt( AA*AA + D*D - 2.*D*AA*Cos((180.-Z)*DEG) );
   return HH - RE;

end; \function

function real Column_density( A, Z ); real A, Z; \integrates density along the line of sight real Sum, D, Delta; [ Sum:= 0.0; D:= 0.0;

   while D<FIN do
       [Delta:= Max(DD, (DD)*D); \adaptive step size to avoid it taking forever:
       Sum:= Sum + Rho(Height(A, Z, D+0.5*Delta))*Delta;
       D:= D + Delta;
       ];
   return Sum;

end; \function

function real Airmass( A, Z ); real A, Z; [ return Column_density( A, Z ) / Column_density( A, 0. ); end; \function

real Z; [Text(0, "Angle 0 m 13700 m^M^J");

Text(0, "------------------------------------^M^J");

Z:= 0.; while Z<=90. do

   [Format(2, 0);  RlOut(0, Z);
   Format(8, 8);   RlOut(0, Airmass(0., Z));
   RlOut(0, Airmass(13700., Z));  CrLf(0);
   Z:= Z + 5.;
   ]

]</lang>

Output:

Angle     0 m              13700 m
------------------------------------
 0       1.00000000       1.00000000
 5       1.00380963       1.00380965
10       1.01538466       1.01538475
15       1.03517744       1.03517765
20       1.06399053       1.06399093
25       1.10305937       1.10306005
30       1.15418974       1.15419083
35       1.21998076       1.21998246
40       1.30418931       1.30419190
45       1.41234169       1.41234567
50       1.55280404       1.55281025
55       1.73875921       1.73876915
60       1.99212000       1.99213665
65       2.35199740       2.35202722
70       2.89531368       2.89537287
75       3.79582352       3.79596149
80       5.53885809       5.53928113
85      10.07896219      10.08115981
90      34.32981136      34.36666557