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 (h(a,z,x)) where (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 SOFIA infrared telescope, which has a cruising altitude of 13,700 meters (about 8.3 miles),
- flying in a specially retrofitted Boeing 747 about four flights a week).
Factor
<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
<lang 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
<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
Perl
<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
<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 ) {
return 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² + d² - 2*d*AA*cos((180-z)*DEG) ); return 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;
}
return $sum;
}
sub airmass ( \a, \z ) {
return column_density( a, z )
/ column_density( a, 0 );
}
say 'Angle 0 m 13700 m'; say '------------------------------------'; for 0, 5 ... 90 -> \z {
say 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
<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 comma(oHT) "meters", w) say '─────┼'copies(center("", w, '─'), 2)'─' y= left(, w-20) /*pad 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*/
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 comma: 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
<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
<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