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.

Write a function that calculates the air mass for an observer at a given altitude a above sea level and zenith angle z.

For this exercise 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.

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 SOFIA infrared telescope, which has a cruising altitude of 13,700 m.

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

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

REXX

Translation of: FreeBASIC

<lang rexx></lang> output

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