Air mass: Difference between revisions
m (added some facts about the SOFIA flying telescope.) |
m (→{{header|Wren}}: Minor tidy) |
||
(35 intermediate revisions by 19 users not shown) | |||
Line 6: | Line 6: | ||
For this task you can assume: |
For this task you can assume: |
||
* The density of Earth's atmosphere is proportional to exp(-a/8500 metres) |
:* The density of Earth's atmosphere is proportional to exp(-a/8500 metres) |
||
* The Earth is a perfect sphere of radius 6731 km. |
:* The Earth is a perfect sphere of radius 6731 km. |
||
;Task: |
;Task: |
||
:* Write a function that calculates the air mass for an observer at a given altitude ''' ''a'' ''' above sea level and zenith angle ''' ''z.'' ''' |
:* 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. |
:* 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), |
:* Do the same for the [https://www.nasa.gov/mission_pages/SOFIA/overview/index.html 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. |
||
<br><br> |
<br><br> |
||
=={{header|11l}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="11l">V DEG = 0.017453292519943295769236907684886127134 |
|||
V RE = 6371000 |
|||
V dd = 0.001 |
|||
V FIN = 10000000 |
|||
F rho(a) |
|||
‘ the density of air as a function of height above sea level ’ |
|||
R exp(-a / 8500.0) |
|||
F height(Float a; z, d) |
|||
‘ |
|||
a = altitude of observer |
|||
z = zenith angle (in degrees) |
|||
d = distance along line of sight |
|||
’ |
|||
R sqrt((:RE + a) ^ 2 + d ^ 2 - 2 * d * (:RE + a) * cos((180 - z) * :DEG)) - :RE |
|||
F column_density(a, z) |
|||
‘ integrates density along the line of sight ’ |
|||
V (dsum, d) = (0.0, 0.0) |
|||
L d < :FIN |
|||
V delta = max(:dd, (:dd) * d) |
|||
dsum += rho(height(a, z, d + 0.5 * delta)) * delta |
|||
d += delta |
|||
R dsum |
|||
F airmass(a, z) |
|||
R column_density(a, z) / column_density(a, 0) |
|||
print("Angle 0 m 13700 m\n "(‘-’ * 36)) |
|||
L(z) (0.<91).step(5) |
|||
print(f:‘{z:3} {airmass(0, z):12.7} {airmass(13700, z):12.7}’)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Angle 0 m 13700 m |
|||
------------------------------------ |
|||
0 1.0000000 1.0000000 |
|||
5 1.0038096 1.0038096 |
|||
10 1.0153847 1.0153848 |
|||
15 1.0351774 1.0351776 |
|||
20 1.0639905 1.0639909 |
|||
25 1.1030594 1.1030601 |
|||
30 1.1541897 1.1541908 |
|||
35 1.2199808 1.2199825 |
|||
40 1.3041893 1.3041919 |
|||
45 1.4123417 1.4123457 |
|||
50 1.5528040 1.5528102 |
|||
55 1.7387592 1.7387692 |
|||
60 1.9921200 1.9921366 |
|||
65 2.3519974 2.3520272 |
|||
70 2.8953137 2.8953729 |
|||
75 3.7958235 3.7959615 |
|||
80 5.5388581 5.5392811 |
|||
85 10.0789622 10.0811598 |
|||
90 34.3298114 34.3666656 |
|||
</pre> |
|||
=={{header|Ada}}== |
|||
{{trans|C}} |
|||
{{works with|Ada 2012}} |
|||
<syntaxhighlight lang="ada"> |
|||
with Ada.Text_IO; use Ada.Text_IO; |
|||
with Ada.Numerics.Generic_Elementary_Functions; |
|||
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO; |
|||
procedure Main is |
|||
subtype double is Long_Float; |
|||
package double_io is new Ada.Text_IO.Float_IO (double); |
|||
use double_io; |
|||
package Elementary_Double is new Ada.Numerics.Generic_Elementary_Functions |
|||
(Float_Type => double); |
|||
use Elementary_Double; |
|||
-- degrees to radians |
|||
Deg : constant := 0.017_453_292_519_943_295_769_236_907_684_886_127_134; |
|||
-- Earth radius in meters |
|||
Re : constant := 6_371_000.0; |
|||
-- integrate in this fraction of the distance already covered |
|||
Dd : constant := 0.001; |
|||
-- integrate only to a height of 10000km. efectively infinity |
|||
Fin : constant := 10_000_000.0; |
|||
function rho (a : double) return double is (Exp (-a / 8_500.0)); |
|||
function height (a : double; z : double; d : double) return double is |
|||
aa : double := Re + a; |
|||
hh : double := |
|||
Sqrt (aa * aa + d * d - 2.0 * d * aa * Cos ((180.0 - z) * Deg)); |
|||
begin |
|||
return hh - Re; |
|||
end height; |
|||
function column_density (a : double; z : double) return double is |
|||
sum : double := 0.0; |
|||
d : double := 0.0; |
|||
d_delta : double; |
|||
begin |
|||
while d < Fin loop |
|||
-- adaptive step size to avoid it taking forever |
|||
d_delta := Dd * d; |
|||
if d_delta < Dd then |
|||
d_delta := Dd; |
|||
end if; |
|||
sum := sum + rho (height (a, z, d + 0.5 * d_delta)) * d_delta; |
|||
d := d + d_delta; |
|||
end loop; |
|||
return sum; |
|||
end column_density; |
|||
function air_mass (a : double; z : double) return double is |
|||
(column_density (a, z) / column_density (a, 0.0)); |
|||
z : double := 0.0; |
|||
begin |
|||
Put_Line ("Angle 0 m 13700 m"); |
|||
Put_Line ("------------------------------------"); |
|||
while z <= 90.0 loop |
|||
Put(Item => Integer(z), Width => 2); |
|||
Put (Item => air_mass (0.0, z), Fore => 8, Aft => 8, Exp => 0); |
|||
Put (Item => air_mass (13_700.0, z), Fore => 8, Aft => 8, Exp => 0); |
|||
New_Line; |
|||
z := z + 5.0; |
|||
end loop; |
|||
end Main; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|AWK}}== |
|||
<syntaxhighlight lang="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) } |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|BASIC}}== |
|||
==={{header|BASIC256}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="freebasic">global RE, dd, LIM |
|||
RE = 6371000 #Earth radius in meters |
|||
dd = 0.001 #integrate in this fraction of the distance already covered |
|||
LIM = 10000000 #integrate only to a height of 10000km, effectively inLIMity |
|||
print "Angle 0 m 13700 m" |
|||
print "------------------------------------" |
|||
for z = 0 to 90 step 5 |
|||
print rjust(z,2); " "; ljust(airmass(0, z),13,"0"); " "; ljust(airmass(13700, z),13,"0") |
|||
next z |
|||
end |
|||
function max(a, b) |
|||
if a > b then return a else return b |
|||
end function |
|||
function rho(a) |
|||
#the density of air as a function of height above sea level |
|||
return exp(-a/8500.0) |
|||
end function |
|||
function height(a, z, d) |
|||
#a = altitude of observer |
|||
#z = zenith angle (in degrees) |
|||
#d = distance along line of sight |
|||
AA = RE + a |
|||
HH = sqr(AA^2 + d^2 - 2*d*AA*cos(radians(180-z))) |
|||
return HH - RE |
|||
end function |
|||
function column_density(a, z) |
|||
#integrates density along the line of sight |
|||
sum = 0.0 |
|||
d = 0.0 |
|||
while d < LIM |
|||
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 |
|||
end while |
|||
return sum |
|||
end function |
|||
function airmass(a, z) |
|||
return column_density(a, z) / column_density(a, 0) |
|||
end function</syntaxhighlight> |
|||
==={{header|FreeBASIC}}=== |
|||
<syntaxhighlight 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 |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
==={{header|True BASIC}}=== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="qbasic">FUNCTION max(a, b) |
|||
IF a > b then LET max = a else LET max = b |
|||
END FUNCTION |
|||
FUNCTION rho(a) |
|||
!the density of air as a function of height above sea level |
|||
LET rho = exp(-a/8500) |
|||
END FUNCTION |
|||
FUNCTION height(a, z, d) |
|||
!a = altitude of observer |
|||
!z = zenith angle (in degrees) |
|||
!d = distance along line of sight |
|||
LET aa = re+a |
|||
LET hh = sqr(aa^2+d^2-2*d*aa*cos((180-z)*deg)) |
|||
LET height = hh-re |
|||
END FUNCTION |
|||
FUNCTION columndensity(a, z) |
|||
!integrates density along the line of sight |
|||
LET sum = 0 |
|||
LET d = 0 |
|||
DO while d < lim |
|||
LET delta = max(dd, (dd)*d) !adaptive step size to avoid it taking forever: |
|||
LET sum = sum+rho(height(a, z, d+.5*delta))*delta |
|||
LET d = d+delta |
|||
LOOP |
|||
LET columndensity = sum |
|||
END FUNCTION |
|||
FUNCTION airmass(a, z) |
|||
LET airmass = columndensity(a, z)/columndensity(a, 0) |
|||
END FUNCTION |
|||
LET deg = .0174532925199433 !degrees to radians |
|||
LET re = 6371000 !Earth radius in meters |
|||
LET dd = .001 !integrate in this fraction of the distance already covered |
|||
LET lim = 10000000 !integrate only to a height of 10000km, effectively infinity |
|||
PRINT "Angle 0 m 13700 m" |
|||
PRINT "------------------------------------" |
|||
FOR z = 0 to 90 step 5 |
|||
PRINT using "## ##.######## ##.########": z, airmass(0, z), airmass(13700, z) |
|||
NEXT z |
|||
END</syntaxhighlight> |
|||
=={{header|C}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="c">#include <math.h> |
|||
#include <stdio.h> |
|||
#define DEG 0.017453292519943295769236907684886127134 // degrees to radians |
|||
#define RE 6371000.0 // Earth radius in meters |
|||
#define DD 0.001 // integrate in this fraction of the distance already covered |
|||
#define FIN 10000000.0 // integrate only to a height of 10000km, effectively infinity |
|||
static double rho(double a) { |
|||
// the density of air as a function of height above sea level |
|||
return exp(-a / 8500.0); |
|||
} |
|||
static double height(double a, double z, double d) { |
|||
// a = altitude of observer |
|||
// z = zenith angle (in degrees) |
|||
// d = distance along line of sight |
|||
double aa = RE + a; |
|||
double hh = sqrt(aa * aa + d * d - 2.0 * d * aa * cos((180 - z) * DEG)); |
|||
return hh - RE; |
|||
} |
|||
static double column_density(double a, double z) { |
|||
// integrates density along the line of sight |
|||
double sum = 0.0, d = 0.0; |
|||
while (d < FIN) { |
|||
// adaptive step size to avoid it taking forever |
|||
double delta = DD * d; |
|||
if (delta < DD) |
|||
delta = DD; |
|||
sum += rho(height(a, z, d + 0.5 * delta)) * delta; |
|||
d += delta; |
|||
} |
|||
return sum; |
|||
} |
|||
static double airmass(double a, double z) { |
|||
return column_density(a, z) / column_density(a, 0.0); |
|||
} |
|||
int main() { |
|||
puts("Angle 0 m 13700 m"); |
|||
puts("------------------------------------"); |
|||
for (double z = 0; z <= 90; z+= 5) { |
|||
printf("%2.0f %11.8f %11.8f\n", |
|||
z, airmass(0.0, z), airmass(13700.0, z)); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Delphi}}== |
|||
{{works with|Delphi|6.0}} |
|||
{{trans|GO}} |
|||
{{libheader|SysUtils,StdCtrls}} |
|||
<syntaxhighlight lang="Delphi"> |
|||
const RE = 6371000; { radius of earth in meters} |
|||
const DD = 0.001; { integrate in this fraction of the distance already covered} |
|||
const FIN = 1e7; { integrate only to a height of 10000km, effectively infinity} |
|||
function rho(a: double): double; |
|||
{ The density of air as a function of height above sea level.} |
|||
begin |
|||
Result:=Exp(-a / 8500); |
|||
end; |
|||
function Radians(degrees: double): double; |
|||
{ Converts degrees to radians} |
|||
begin |
|||
Result:= degrees * Pi / 180 |
|||
end; |
|||
function Height(A, Z, D: double): double; |
|||
{ a = altitude of observer} |
|||
{ z = zenith angle (in degrees)} |
|||
{ d = distance along line of sight} |
|||
var AA,HH: double; |
|||
begin |
|||
AA := RE + A; |
|||
HH := Sqrt(AA*AA + D*D - 2*D*AA*Cos(Radians(180-z))); |
|||
Result:= HH - RE; |
|||
end; |
|||
function ColumnDensity(A, Z: double): double; |
|||
{ Integrates density along the line of sight.} |
|||
var Sum,D,Delta: double; |
|||
begin |
|||
Sum := 0.0; |
|||
D := 0.0; |
|||
while D < FIN do |
|||
begin |
|||
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; |
|||
end; |
|||
Result:= Sum; |
|||
end; |
|||
function AirMass(A, Z: double): double; |
|||
begin |
|||
Result:= ColumnDensity(A, Z) / ColumnDensity(a, 0); |
|||
end; |
|||
procedure ShowAirMass(Memo: TMemo); |
|||
var Z: integer; |
|||
begin |
|||
Memo.Lines.Add('Angle 0 m 13700 m'); |
|||
Memo.Lines.Add('------------------------------------'); |
|||
Z:=0; |
|||
while Z<=90 do |
|||
begin |
|||
Memo.Lines.Add(Format('%2d %11.8f %11.8f', [z, airmass(0, Z), airmass(13700, Z)])); |
|||
Z:=Z+5; |
|||
end; |
|||
end; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
Elapsed Time: 189.304 ms. |
|||
</pre> |
|||
=={{header|EasyLang}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang=easylang> |
|||
func rho a . |
|||
return pow 2.718281828459 (-a / 8500) |
|||
. |
|||
func height a z d . |
|||
AA = 6371000 + a |
|||
HH = sqrt (AA * AA + d * d - 2 * d * AA * cos (180 - z)) |
|||
return HH - 6371000 |
|||
. |
|||
func density a z . |
|||
while d < 10000000 |
|||
delta = higher 0.001 (0.001 * d) |
|||
sum += rho height a z (d + 0.5 * delta) * delta |
|||
d += delta |
|||
. |
|||
return sum |
|||
. |
|||
func airmass a z . |
|||
return density a z / density a 0 |
|||
. |
|||
numfmt 8 2 |
|||
print "Angle 0 m 13700 m" |
|||
print "------------------------" |
|||
for z = 0 step 5 to 90 |
|||
print z & " " & airmass 0 z & " " & airmass 13700 z |
|||
. |
|||
</syntaxhighlight> |
|||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
{{works with|Factor|0.99 2021-02-05}} |
{{works with|Factor|0.99 2021-02-05}} |
||
< |
<syntaxhighlight lang="factor">USING: formatting io kernel math math.functions math.order |
||
math.ranges math.trig sequences ; |
math.ranges math.trig sequences ; |
||
Line 53: | Line 646: | ||
dup [ 0 swap airmass ] [ 13700 swap airmass ] bi |
dup [ 0 swap airmass ] [ 13700 swap airmass ] bi |
||
"%2d %15.8f %17.8f\n" printf |
"%2d %15.8f %17.8f\n" printf |
||
] each</ |
] each</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 79: | Line 672: | ||
</pre> |
</pre> |
||
=={{header| |
=={{header|Go}}== |
||
{{trans|FreeBASIC}} |
|||
<lang freebasic> |
|||
<syntaxhighlight lang="go">package main |
|||
#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) |
|||
import ( |
|||
function rho(a as double) as double |
|||
"fmt" |
|||
'the density of air as a function of height above sea level |
|||
"math" |
|||
return exp(-a/8500.0) |
|||
) |
|||
end function |
|||
const ( |
|||
function height( a as double, z as double, d as double ) as double |
|||
RE = 6371000 // radius of earth in meters |
|||
DD = 0.001 // integrate in this fraction of the distance already covered |
|||
'z = zenith angle (in degrees) |
|||
FIN = 1e7 // integrate only to a height of 10000km, effectively infinity |
|||
'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 |
|||
// The density of air as a function of height above sea level. |
|||
function column_density( a as double, z as double ) as double |
|||
func rho(a float64) float64 { return math.Exp(-a / 8500) } |
|||
'integrates density along the line of sight |
|||
dim as double sum = 0.0, d = 0.0, delta |
|||
// Converts degrees to radians |
|||
while d<FIN |
|||
func radians(degrees float64) float64 { return degrees * math.Pi / 180 } |
|||
delta = max(dd, (dd)*d) 'adaptive step size to avoid it taking forever: |
|||
sum += rho(height(a, z, d+0.5*delta))*delta |
|||
// a = altitude of observer |
|||
// z = zenith angle (in degrees) |
|||
// d = distance along line of sight |
|||
func height(a, z, d float64) float64 { |
|||
aa := RE + a |
|||
hh := math.Sqrt(aa*aa + d*d - 2*d*aa*math.Cos(radians(180-z))) |
|||
return hh - RE |
|||
} |
|||
// Integrates density along the line of sight. |
|||
func columnDensity(a, z float64) float64 { |
|||
sum := 0.0 |
|||
d := 0.0 |
|||
for d < FIN { |
|||
delta := math.Max(DD, DD*d) // adaptive step size to avoid it taking forever |
|||
sum += rho(height(a, z, d+0.5*delta)) * delta |
|||
d += delta |
d += delta |
||
} |
|||
return sum |
return sum |
||
} |
|||
end function |
|||
func airmass(a, z float64) float64 { |
|||
return |
return columnDensity(a, z) / columnDensity(a, 0) |
||
} |
|||
end function |
|||
func main() { |
|||
print "Angle 0 m 13700 m" |
|||
fmt.Println("Angle 0 m 13700 m") |
|||
print "------------------------------------" |
|||
fmt.Println("------------------------------------") |
|||
for z as double = 0 to 90 step 5.0 |
|||
for z := 0; z <= 90; z += 5 { |
|||
print using "## ##.######## ##.########";z;airmass(0, z);airmass(13700, z) |
|||
fz := float64(z) |
|||
next z |
|||
fmt.Printf("%2d %11.8f %11.8f\n", z, airmass(0, fz), airmass(13700, fz)) |
|||
</lang> |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Java}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="java">public class AirMass { |
|||
public static void main(String[] args) { |
|||
System.out.println("Angle 0 m 13700 m"); |
|||
System.out.println("------------------------------------"); |
|||
for (double z = 0; z <= 90; z+= 5) { |
|||
System.out.printf("%2.0f %11.8f %11.8f\n", |
|||
z, airmass(0.0, z), airmass(13700.0, z)); |
|||
} |
|||
} |
|||
private static double rho(double a) { |
|||
// the density of air as a function of height above sea level |
|||
return Math.exp(-a / 8500.0); |
|||
} |
|||
private static double height(double a, double z, double d) { |
|||
// a = altitude of observer |
|||
// z = zenith angle (in degrees) |
|||
// d = distance along line of sight |
|||
double aa = RE + a; |
|||
double hh = Math.sqrt(aa * aa + d * d - 2.0 * d * aa * Math.cos(Math.toRadians(180 - z))); |
|||
return hh - RE; |
|||
} |
|||
private static double columnDensity(double a, double z) { |
|||
// integrates density along the line of sight |
|||
double sum = 0.0, d = 0.0; |
|||
while (d < FIN) { |
|||
// adaptive step size to avoid it taking forever |
|||
double delta = Math.max(DD * d, DD); |
|||
sum += rho(height(a, z, d + 0.5 * delta)) * delta; |
|||
d += delta; |
|||
} |
|||
return sum; |
|||
} |
|||
private static double airmass(double a, double z) { |
|||
return columnDensity(a, z) / columnDensity(a, 0.0); |
|||
} |
|||
private static final double RE = 6371000.0; // Earth radius in meters |
|||
private static final double DD = 0.001; // integrate in this fraction of the distance already covered |
|||
private static final double FIN = 10000000.0; // integrate only to a height of 10000km, effectively infinity |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|jq}}== |
|||
'''Adapted from [[#Wren|Wren]]''' |
|||
{{works with|jq}} |
|||
'''Works with gojq, the Go implementation of jq''' |
|||
'''Preliminaries''' |
|||
<syntaxhighlight lang="jq">def pi: 4 * (1|atan); |
|||
def radians: . * pi / 180; |
|||
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .; |
|||
# Input: a number |
|||
# Output: a string with $digits fractional decimal digits, with proper rounding |
|||
def fmt($width; $digits): |
|||
. as $in |
|||
| tostring |
|||
| index(".") as $ix |
|||
| if test("[eE]") then . |
|||
elif $ix |
|||
then pow(10; $digits) as $p |
|||
| ($in * $p | round | tostring) as $s |
|||
| if test("[eE]") then $s |
|||
else ($s | index(".")) as $ix |
|||
| if $ix then $s[:$ix + 1] + $s[$ix+1: $ix+1+$digits] |
|||
else $s[:-$digits] + "." + $s[-$digits:] |
|||
end |
|||
end |
|||
else . + "." + "0" * digits |
|||
end |
|||
| lpad($width);</syntaxhighlight> |
|||
'''Physics''' |
|||
<syntaxhighlight lang="jq"># constants |
|||
def RE: 6371000; # radius of earth in meters |
|||
def DD: 0.001; # integrate in this fraction of the distance already covered |
|||
def FIN: 1e7; # integrate only to a height of 10000km, effectively infinity |
|||
# The density of air as a function of height above sea level. |
|||
def rho: (-./8500) | exp; |
|||
# a = altitude of observer (in m) |
|||
# z = zenith angle (in degrees) |
|||
# d = distance along line of sight (in m) |
|||
def height($a; $z; $d): |
|||
(RE + $a) as $aa |
|||
| (($aa * $aa + $d * $d - 2 * $d * $aa * ((180-$z)|radians|cos) )|sqrt ) - RE; |
|||
# Integrates density along the line of sight. |
|||
def columnDensity($a; $z): |
|||
{ sum: 0, d: 0 } |
|||
| until (.d >= FIN; |
|||
([DD, DD * .d] | max) as $delta # adaptive step size to avoid it taking forever |
|||
| .sum = .sum + ((height($a; $z; .d + 0.5 * $delta))|rho) * $delta |
|||
| .d += $delta ) |
|||
| .sum ; |
|||
def airmass(a; z): columnDensity(a; z) / columnDensity(a; 0); |
|||
"Angle 0 m 13700 m", |
|||
"------------------------------------", |
|||
( range(0; 91; 5) |
|||
| "\(lpad(2)) \(airmass(0; .)|fmt(11;8)) \(airmass(13700; .)|fmt(11;8))" )</syntaxhighlight> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 151: | Line 914: | ||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
< |
<syntaxhighlight lang="julia">using Printf |
||
const DEG = 0.017453292519943295769236907684886127134 # degrees to radians |
const DEG = 0.017453292519943295769236907684886127134 # degrees to radians |
||
Line 183: | Line 946: | ||
@printf("%2d %11.8f %11.8f\n", z, airmass(0, z), airmass(13700, z)) |
@printf("%2d %11.8f %11.8f\n", z, airmass(0, z), airmass(13700, z)) |
||
end |
end |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
Angle 0 m 13700 m |
Angle 0 m 13700 m |
||
Line 207: | Line 970: | ||
90 34.32981136 34.36666557 |
90 34.32981136 34.36666557 |
||
</pre> |
</pre> |
||
=={{header|Nim}}== |
|||
{{trans|Wren}} |
|||
<syntaxhighlight 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</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|Perl}}== |
|||
{{trans|Raku}} |
|||
<syntaxhighlight 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); |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">constant</span> <span style="color: #000000;">RE</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">6371000</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">// radius of earth in meters</span> |
<span style="color: #008080;">constant</span> <span style="color: #000000;">RE</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">6371000</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">// radius of earth in meters</span> |
||
<span style="color: #000000;">DD</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.001</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">// integrate in this fraction of the distance already covered</span> |
<span style="color: #000000;">DD</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0.001</span><span style="color: #0000FF;">,</span> <span style="color: #000080;font-style:italic;">// integrate in this fraction of the distance already covered</span> |
||
Line 245: | Line 1,146: | ||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d %11.8f %11.8f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">airmass</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">airmass</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13700</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)})</span> |
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d %11.8f %11.8f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">z</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">airmass</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">),</span> <span style="color: #000000;">airmass</span><span style="color: #0000FF;">(</span><span style="color: #000000;">13700</span><span style="color: #0000FF;">,</span><span style="color: #000000;">z</span><span style="color: #0000FF;">)})</span> |
||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 270: | Line 1,171: | ||
90 34.32981136 34.36666557 |
90 34.32981136 34.36666557 |
||
</pre> |
</pre> |
||
=={{header|Python}}== |
|||
<syntaxhighlight lang="python">""" Rosetta Code task: Air_mass """ |
|||
from math import sqrt, cos, exp |
|||
DEG = 0.017453292519943295769236907684886127134 # degrees to radians |
|||
RE = 6371000 # Earth radius in meters |
|||
dd = 0.001 # integrate in this fraction of the distance already covered |
|||
FIN = 10000000 # integrate only to a height of 10000km, effectively infinity |
|||
def rho(a): |
|||
""" the density of air as a function of height above sea level """ |
|||
return exp(-a / 8500.0) |
|||
def height(a, z, d): |
|||
""" |
|||
a = altitude of observer |
|||
z = zenith angle (in degrees) |
|||
d = distance along line of sight |
|||
""" |
|||
return sqrt((RE + a)**2 + d**2 - 2 * d * (RE + a) * cos((180 - z) * DEG)) - RE |
|||
def column_density(a, z): |
|||
""" integrates density along the line of sight """ |
|||
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 |
|||
return dsum |
|||
def airmass(a, z): |
|||
return column_density(a, z) / column_density(a, 0) |
|||
print('Angle 0 m 13700 m\n', '-' * 36) |
|||
for z in range(0, 91, 5): |
|||
print(f"{z: 3d} {airmass(0, z): 12.7f} {airmass(13700, z): 12.7f}") |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
Angle 0 m 13700 m |
|||
------------------------------------ |
|||
0 1.0000000 1.0000000 |
|||
5 1.0038096 1.0038096 |
|||
10 1.0153847 1.0153848 |
|||
15 1.0351774 1.0351776 |
|||
20 1.0639905 1.0639909 |
|||
25 1.1030594 1.1030601 |
|||
30 1.1541897 1.1541908 |
|||
35 1.2199808 1.2199825 |
|||
40 1.3041893 1.3041919 |
|||
45 1.4123417 1.4123457 |
|||
50 1.5528040 1.5528102 |
|||
55 1.7387592 1.7387692 |
|||
60 1.9921200 1.9921366 |
|||
65 2.3519974 2.3520272 |
|||
70 2.8953137 2.8953729 |
|||
75 3.7958235 3.7959615 |
|||
80 5.5388581 5.5392811 |
|||
85 10.0789622 10.0811598 |
|||
90 34.3298114 34.3666656 |
|||
</pre> |
|||
=={{header|Raku}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="raku" line>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); |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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</pre> |
|||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
{{trans|FreeBASIC}} |
{{trans|FreeBASIC}} |
||
< |
<syntaxhighlight 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*/ |
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*/ |
parse arg aLO aHI aBY oHT . /*obtain optional arguments from the CL*/ |
||
Line 281: | Line 1,308: | ||
if oHT=='' | oHT=="," then oHT= 13700 /* " " " " " " */ |
if oHT=='' | oHT=="," then oHT= 13700 /* " " " " " " */ |
||
w= 30; @ama= 'air mass at' /*column width for the two air_masses. */ |
w= 30; @ama= 'air mass at' /*column width for the two air_masses. */ |
||
say 'angle|'center(@ama "sea level", w) |
say 'angle|'center(@ama "sea level", w) center(@ama commas(oHT) 'meters', w) /*title*/ |
||
say |
say "─────┼"copies(center('', w, "─"), 2)'─' /*display the title sep for the output.*/ |
||
y= left('', w-20) /* |
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) |
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 |
say center(j, 5)'│'right( format(am0, , 8), w-10)y right( format(amht, , 8), w-10)y |
||
end /*j*/ |
end /*j*/ |
||
say "─────┴"copies(center('', w, "─"), 2)'─' /*display the foot separator for output*/ |
|||
exit 0 /*stick a fork in it, we're all done. */ |
exit 0 /*stick a fork in it, we're all done. */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
Line 296: | Line 1,325: | ||
r2r: return arg(1) // (pi() * 2) /*normalize radians ──► a unit circle. */ |
r2r: return arg(1) // (pi() * 2) /*normalize radians ──► a unit circle. */ |
||
e: e= 2.718281828459045235360287471352662497757247093699959574966967627724; return e |
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) |
cos: procedure; parse arg x; x= r2r(x); a= abs(x); numeric fuzz min(6, digits() - 3) |
||
Line 317: | Line 1,346: | ||
aa= earthRad + a; return sqrt(aa**2 + d**2 - 2*d*aa*cos( d2r(180-z) ) ) - earthRad |
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 |
colD: procedure; parse arg a,z; sum= 0; d= 0; dd= .001; infinity= 10000000 |
||
do while d<infinity; delta= max(dd, dd*d) |
do while d<infinity; delta= max(dd, dd*d) |
||
sum= sum + rho( elev(a, z, d + 0.5*delta) ) * delta; d= d + delta |
sum= sum + rho( elev(a, z, d + 0.5*delta) ) * delta; d= d + delta |
||
end /*while*/ |
end /*while*/ |
||
return sum</ |
return sum</syntaxhighlight> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |
||
Line 345: | Line 1,374: | ||
85 │ 10.07896219 10.08115981 |
85 │ 10.07896219 10.08115981 |
||
90 │ 34.32981136 34.36666557 |
90 │ 34.32981136 34.36666557 |
||
─────┴───────────────────────────────────────────────────────────── |
|||
</pre> |
</pre> |
||
=={{header|RPL}}== |
|||
{{trans|FreeBASIC}} |
|||
≪ → a ≪ a NEG 8500 / EXP ≫ ≫ ‘'''RHO'''’ STO |
|||
≪ ''''RHO'''(√(aa^2+D^2-2*aa*D*COS(180-Z))-re)' EVAL |
|||
≫ ‘'''COLD'''’ STO |
|||
≪ |
|||
6371000 3 PICK OVER + → a z re aa |
|||
≪ DEG |
|||
z 'Z' STO ''''COLD'''' { D 0 1E7 } 1E-7 ∫ DROP |
|||
0 'Z' STO ''''COLD'''' { D 0 1E7 } 1E-7 ∫ DROP / |
|||
≫ ≫ ‘'''AM'''’ STO |
|||
≪ { } 0 90 '''FOR''' z z + z '''AM''' + 5 '''STEP''' ≫ |
|||
{{out}} |
|||
<pre> |
|||
1: { 0 1 |
|||
5 1.00380686363 |
|||
10 1.01537368745 |
|||
15 1.03515302646 |
|||
20 1.06394782383 |
|||
25 1.10299392042 |
|||
30 1.15409753978 |
|||
35 1.21985818096 |
|||
40 1.30403285254 |
|||
45 1.41214767977 |
|||
50 1.55256798138 |
|||
55 1.73847475415 |
|||
60 1.99177718552 |
|||
65 2.35157928673 |
|||
70 2.89478919419 |
|||
75 3.79512945489 |
|||
80 5.53784169364 |
|||
85 10.0771111633 |
|||
90 34.3235064081 } |
|||
</pre> |
|||
=={{header|Rust}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="rust">const RE: f64 = 6371000.0; // Earth radius in meters |
|||
const DD: f64 = 0.001; // integrate in this fraction of the distance already covered |
|||
const FIN: f64 = 10000000.0; // integrate only to a height of 10000km, effectively infinity |
|||
fn rho(a: f64) -> f64 { |
|||
// the density of air as a function of height above sea level |
|||
(-a / 8500.0).exp() |
|||
} |
|||
fn height(a: f64, z: f64, d: f64) -> f64 { |
|||
// a = altitude of observer |
|||
// z = zenith angle (in degrees) |
|||
// d = distance along line of sight |
|||
let aa = RE + a; |
|||
let hh = (aa * aa + d * d - 2.0 * d * aa * (180.0 - z).to_radians().cos()).sqrt(); |
|||
hh - RE |
|||
} |
|||
fn column_density(a: f64, z: f64) -> f64 { |
|||
// integrates density along the line of sight |
|||
let mut sum = 0.0; |
|||
let mut d = 0.0; |
|||
while d < FIN { |
|||
// adaptive step size to avoid it taking forever |
|||
let mut delta = DD * d; |
|||
if delta < DD { |
|||
delta = DD; |
|||
} |
|||
sum += rho(height(a, z, d + 0.5 * delta)) * delta; |
|||
d += delta; |
|||
} |
|||
sum |
|||
} |
|||
fn airmass(a: f64, z: f64) -> f64 { |
|||
column_density(a, z) / column_density(a, 0.0) |
|||
} |
|||
fn main() { |
|||
println!("Angle 0 m 13700 m"); |
|||
println!("------------------------------------"); |
|||
for a in (0..=90).step_by(5) { |
|||
let z = a as f64; |
|||
println!( |
|||
"{:2} {:11.8} {:11.8}", |
|||
z, |
|||
airmass(0.0, z), |
|||
airmass(13700.0, z) |
|||
); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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 |
|||
</pre> |
|||
=={{header|Seed7}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight lang="seed7">$ include "seed7_05.s7i"; |
|||
include "float.s7i"; |
|||
include "math.s7i"; |
|||
const float: DEG is 0.017453292519943295769236907684886127134; #degrees to radians |
|||
const float: RE is 6371000.0; #Earth radius in meters |
|||
const float: dd is 0.001; #integrate in this fraction of the distance already covered |
|||
const float: FIN is 10000000.0; #integrate only to a height of 10000km, effectively infinity |
|||
const func float: rho (in float: a) is |
|||
#the density of air as a function of height above sea level |
|||
return exp(-a / 8500.0); |
|||
const func float: height (in float: a, in float: z, in float: d) is func |
|||
#a is altitude of observer |
|||
#z is zenith angle (in degrees) |
|||
#d is distance along line of sight |
|||
result |
|||
var float: r is 0.0; |
|||
local |
|||
var float: AA is 0.0; |
|||
var float: HH is 0.0; |
|||
begin |
|||
AA := RE + a; |
|||
HH := sqrt( AA ** 2.0 + d ** 2.0 - 2.0 * d * AA * cos((180.0 - z) * DEG) ); |
|||
r := HH - RE; |
|||
end func; |
|||
const func float: columnDensity (in float: a, in float: z) is func |
|||
#integrates density along line of sight |
|||
result |
|||
var float: sum is 0.0; |
|||
local |
|||
var float: d is 0.0; |
|||
var float: delta is 0.0; |
|||
begin |
|||
while d < FIN do |
|||
delta := max(dd, dd * d); #adaptive step size to avoid taking it forever |
|||
sum +:= rho(height(a, z, d + 0.5 * delta)) * delta; |
|||
d +:= delta; |
|||
end while; |
|||
end func; |
|||
const func float: airmass (in float: a, in float: z) is |
|||
return columnDensity(a, z) / columnDensity(a, 0.0); |
|||
const proc: main is func |
|||
local |
|||
var integer: zz is 0; |
|||
var float: z is 0.0; |
|||
begin |
|||
writeln("Angle 0 m 13700 m"); |
|||
writeln("------------------------------------"); |
|||
for zz range 0 to 90 step 5 do |
|||
z := flt(zz); |
|||
write(z lpad 4); |
|||
write(airmass(0.0, z) digits 8 lpad 15); |
|||
writeln(airmass(13700.0, z) digits 8 lpad 17); |
|||
end for; |
|||
end func;</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Angle 0 m 13700 m |
|||
------------------------------------ |
|||
0.0 1.00000000 1.00000000 |
|||
5.0 1.00380963 1.00380965 |
|||
10.0 1.01538466 1.01538475 |
|||
15.0 1.03517744 1.03517765 |
|||
20.0 1.06399053 1.06399093 |
|||
25.0 1.10305937 1.10306005 |
|||
30.0 1.15418974 1.15419083 |
|||
35.0 1.21998076 1.21998246 |
|||
40.0 1.30418931 1.30419190 |
|||
45.0 1.41234169 1.41234567 |
|||
50.0 1.55280404 1.55281025 |
|||
55.0 1.73875921 1.73876915 |
|||
60.0 1.99212000 1.99213665 |
|||
65.0 2.35199740 2.35202722 |
|||
70.0 2.89531368 2.89537287 |
|||
75.0 3.79582352 3.79596149 |
|||
80.0 5.53885809 5.53928113 |
|||
85.0 10.07896219 10.08115981 |
|||
90.0 34.32981136 34.36666557 |
|||
</pre> |
|||
=={{header|Swift}}== |
|||
{{trans|Go}} |
|||
<syntaxhighlight lang="swift">import Foundation |
|||
extension Double { |
|||
var radians: Double { self * .pi / 180 } |
|||
} |
|||
func columnDensity(_ a: Double, _ z: Double) -> Double { |
|||
func rho(_ a: Double) -> Double { |
|||
exp(-a / 8500) |
|||
} |
|||
func height(_ d: Double) -> Double { |
|||
let aa = 6_371_000 + a |
|||
let hh = aa * aa + d * d - 2 * d * aa * cos((180 - z).radians) |
|||
return hh.squareRoot() - 6_371_000 |
|||
} |
|||
var sum = 0.0 |
|||
var d = 0.0 |
|||
while d < 1e7 { |
|||
let delta = max(0.001, 0.001 * d) |
|||
sum += rho(height(d + 0.5 * delta)) * delta |
|||
d += delta |
|||
} |
|||
return sum |
|||
} |
|||
func airMass(altitude: Double, zenith: Double) -> Double { |
|||
return columnDensity(altitude, zenith) / columnDensity(altitude, 0) |
|||
} |
|||
print("Angle 0 m 13700 m") |
|||
print("------------------------------------") |
|||
for z in stride(from: 0.0, through: 90.0, by: 5.0) { |
|||
let air = String( |
|||
format: "%2.0f %11.8f %11.8f", |
|||
z, |
|||
airMass(altitude: 0, zenith: z), |
|||
airMass(altitude: 13700, zenith: z) |
|||
) |
|||
print(air) |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
Line 351: | Line 1,661: | ||
{{libheader|Wren-math}} |
{{libheader|Wren-math}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
< |
<syntaxhighlight lang="wren">import "./math" for Math |
||
import "/fmt" for Fmt |
import "./fmt" for Fmt |
||
// constants |
// constants |
||
Line 360: | Line 1,670: | ||
// The density of air as a function of height above sea level. |
// The density of air as a function of height above sea level. |
||
var rho = Fn.new { |a| |
var rho = Fn.new { |a| (-a/8500).exp } |
||
// a = altitude of observer |
// a = altitude of observer |
||
Line 376: | Line 1,686: | ||
var d = 0 |
var d = 0 |
||
while (d < FIN) { |
while (d < FIN) { |
||
var delta = |
var delta = DD.max(DD * d) // adaptive step size to avoid it taking forever |
||
sum = sum + rho.call(height.call(a, z, d + 0.5 * delta)) * delta |
sum = sum + rho.call(height.call(a, z, d + 0.5 * delta)) * delta |
||
d = d + delta |
d = d + delta |
||
Line 391: | Line 1,701: | ||
Fmt.print("$2d $11.8f $11.8f", z, airmass.call(0, z), airmass.call(13700, z)) |
Fmt.print("$2d $11.8f $11.8f", z, airmass.call(0, z), airmass.call(13700, z)) |
||
z = z + 5 |
z = z + 5 |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|XPL0}}== |
|||
{{trans|FreeBASIC}} |
|||
<syntaxhighlight 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.; |
|||
] |
|||
]</syntaxhighlight> |
|||
{{out}} |
{{out}} |
Latest revision as of 09:42, 6 November 2023
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 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.
11l
V DEG = 0.017453292519943295769236907684886127134
V RE = 6371000
V dd = 0.001
V FIN = 10000000
F rho(a)
‘ the density of air as a function of height above sea level ’
R exp(-a / 8500.0)
F height(Float a; z, d)
‘
a = altitude of observer
z = zenith angle (in degrees)
d = distance along line of sight
’
R sqrt((:RE + a) ^ 2 + d ^ 2 - 2 * d * (:RE + a) * cos((180 - z) * :DEG)) - :RE
F column_density(a, z)
‘ integrates density along the line of sight ’
V (dsum, d) = (0.0, 0.0)
L d < :FIN
V delta = max(:dd, (:dd) * d)
dsum += rho(height(a, z, d + 0.5 * delta)) * delta
d += delta
R dsum
F airmass(a, z)
R column_density(a, z) / column_density(a, 0)
print("Angle 0 m 13700 m\n "(‘-’ * 36))
L(z) (0.<91).step(5)
print(f:‘{z:3} {airmass(0, z):12.7} {airmass(13700, z):12.7}’)
- Output:
Angle 0 m 13700 m ------------------------------------ 0 1.0000000 1.0000000 5 1.0038096 1.0038096 10 1.0153847 1.0153848 15 1.0351774 1.0351776 20 1.0639905 1.0639909 25 1.1030594 1.1030601 30 1.1541897 1.1541908 35 1.2199808 1.2199825 40 1.3041893 1.3041919 45 1.4123417 1.4123457 50 1.5528040 1.5528102 55 1.7387592 1.7387692 60 1.9921200 1.9921366 65 2.3519974 2.3520272 70 2.8953137 2.8953729 75 3.7958235 3.7959615 80 5.5388581 5.5392811 85 10.0789622 10.0811598 90 34.3298114 34.3666656
Ada
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
procedure Main is
subtype double is Long_Float;
package double_io is new Ada.Text_IO.Float_IO (double);
use double_io;
package Elementary_Double is new Ada.Numerics.Generic_Elementary_Functions
(Float_Type => double);
use Elementary_Double;
-- degrees to radians
Deg : constant := 0.017_453_292_519_943_295_769_236_907_684_886_127_134;
-- Earth radius in meters
Re : constant := 6_371_000.0;
-- integrate in this fraction of the distance already covered
Dd : constant := 0.001;
-- integrate only to a height of 10000km. efectively infinity
Fin : constant := 10_000_000.0;
function rho (a : double) return double is (Exp (-a / 8_500.0));
function height (a : double; z : double; d : double) return double is
aa : double := Re + a;
hh : double :=
Sqrt (aa * aa + d * d - 2.0 * d * aa * Cos ((180.0 - z) * Deg));
begin
return hh - Re;
end height;
function column_density (a : double; z : double) return double is
sum : double := 0.0;
d : double := 0.0;
d_delta : double;
begin
while d < Fin loop
-- adaptive step size to avoid it taking forever
d_delta := Dd * d;
if d_delta < Dd then
d_delta := Dd;
end if;
sum := sum + rho (height (a, z, d + 0.5 * d_delta)) * d_delta;
d := d + d_delta;
end loop;
return sum;
end column_density;
function air_mass (a : double; z : double) return double is
(column_density (a, z) / column_density (a, 0.0));
z : double := 0.0;
begin
Put_Line ("Angle 0 m 13700 m");
Put_Line ("------------------------------------");
while z <= 90.0 loop
Put(Item => Integer(z), Width => 2);
Put (Item => air_mass (0.0, z), Fore => 8, Aft => 8, Exp => 0);
Put (Item => air_mass (13_700.0, z), Fore => 8, Aft => 8, Exp => 0);
New_Line;
z := z + 5.0;
end loop;
end Main;
- 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
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) }
- 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
BASIC
BASIC256
global RE, dd, LIM
RE = 6371000 #Earth radius in meters
dd = 0.001 #integrate in this fraction of the distance already covered
LIM = 10000000 #integrate only to a height of 10000km, effectively inLIMity
print "Angle 0 m 13700 m"
print "------------------------------------"
for z = 0 to 90 step 5
print rjust(z,2); " "; ljust(airmass(0, z),13,"0"); " "; ljust(airmass(13700, z),13,"0")
next z
end
function max(a, b)
if a > b then return a else return b
end function
function rho(a)
#the density of air as a function of height above sea level
return exp(-a/8500.0)
end function
function height(a, z, d)
#a = altitude of observer
#z = zenith angle (in degrees)
#d = distance along line of sight
AA = RE + a
HH = sqr(AA^2 + d^2 - 2*d*AA*cos(radians(180-z)))
return HH - RE
end function
function column_density(a, z)
#integrates density along the line of sight
sum = 0.0
d = 0.0
while d < LIM
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
end while
return sum
end function
function airmass(a, z)
return column_density(a, z) / column_density(a, 0)
end function
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
- 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
True BASIC
FUNCTION max(a, b)
IF a > b then LET max = a else LET max = b
END FUNCTION
FUNCTION rho(a)
!the density of air as a function of height above sea level
LET rho = exp(-a/8500)
END FUNCTION
FUNCTION height(a, z, d)
!a = altitude of observer
!z = zenith angle (in degrees)
!d = distance along line of sight
LET aa = re+a
LET hh = sqr(aa^2+d^2-2*d*aa*cos((180-z)*deg))
LET height = hh-re
END FUNCTION
FUNCTION columndensity(a, z)
!integrates density along the line of sight
LET sum = 0
LET d = 0
DO while d < lim
LET delta = max(dd, (dd)*d) !adaptive step size to avoid it taking forever:
LET sum = sum+rho(height(a, z, d+.5*delta))*delta
LET d = d+delta
LOOP
LET columndensity = sum
END FUNCTION
FUNCTION airmass(a, z)
LET airmass = columndensity(a, z)/columndensity(a, 0)
END FUNCTION
LET deg = .0174532925199433 !degrees to radians
LET re = 6371000 !Earth radius in meters
LET dd = .001 !integrate in this fraction of the distance already covered
LET lim = 10000000 !integrate only to a height of 10000km, effectively infinity
PRINT "Angle 0 m 13700 m"
PRINT "------------------------------------"
FOR z = 0 to 90 step 5
PRINT using "## ##.######## ##.########": z, airmass(0, z), airmass(13700, z)
NEXT z
END
C
#include <math.h>
#include <stdio.h>
#define DEG 0.017453292519943295769236907684886127134 // degrees to radians
#define RE 6371000.0 // Earth radius in meters
#define DD 0.001 // integrate in this fraction of the distance already covered
#define FIN 10000000.0 // integrate only to a height of 10000km, effectively infinity
static double rho(double a) {
// the density of air as a function of height above sea level
return exp(-a / 8500.0);
}
static double height(double a, double z, double d) {
// a = altitude of observer
// z = zenith angle (in degrees)
// d = distance along line of sight
double aa = RE + a;
double hh = sqrt(aa * aa + d * d - 2.0 * d * aa * cos((180 - z) * DEG));
return hh - RE;
}
static double column_density(double a, double z) {
// integrates density along the line of sight
double sum = 0.0, d = 0.0;
while (d < FIN) {
// adaptive step size to avoid it taking forever
double delta = DD * d;
if (delta < DD)
delta = DD;
sum += rho(height(a, z, d + 0.5 * delta)) * delta;
d += delta;
}
return sum;
}
static double airmass(double a, double z) {
return column_density(a, z) / column_density(a, 0.0);
}
int main() {
puts("Angle 0 m 13700 m");
puts("------------------------------------");
for (double z = 0; z <= 90; z+= 5) {
printf("%2.0f %11.8f %11.8f\n",
z, airmass(0.0, z), airmass(13700.0, z));
}
}
- 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
Delphi
const RE = 6371000; { radius of earth in meters}
const DD = 0.001; { integrate in this fraction of the distance already covered}
const FIN = 1e7; { integrate only to a height of 10000km, effectively infinity}
function rho(a: double): double;
{ The density of air as a function of height above sea level.}
begin
Result:=Exp(-a / 8500);
end;
function Radians(degrees: double): double;
{ Converts degrees to radians}
begin
Result:= degrees * Pi / 180
end;
function Height(A, Z, D: double): double;
{ a = altitude of observer}
{ z = zenith angle (in degrees)}
{ d = distance along line of sight}
var AA,HH: double;
begin
AA := RE + A;
HH := Sqrt(AA*AA + D*D - 2*D*AA*Cos(Radians(180-z)));
Result:= HH - RE;
end;
function ColumnDensity(A, Z: double): double;
{ Integrates density along the line of sight.}
var Sum,D,Delta: double;
begin
Sum := 0.0;
D := 0.0;
while D < FIN do
begin
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;
end;
Result:= Sum;
end;
function AirMass(A, Z: double): double;
begin
Result:= ColumnDensity(A, Z) / ColumnDensity(a, 0);
end;
procedure ShowAirMass(Memo: TMemo);
var Z: integer;
begin
Memo.Lines.Add('Angle 0 m 13700 m');
Memo.Lines.Add('------------------------------------');
Z:=0;
while Z<=90 do
begin
Memo.Lines.Add(Format('%2d %11.8f %11.8f', [z, airmass(0, Z), airmass(13700, Z)]));
Z:=Z+5;
end;
end;
- 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 Elapsed Time: 189.304 ms.
EasyLang
func rho a .
return pow 2.718281828459 (-a / 8500)
.
func height a z d .
AA = 6371000 + a
HH = sqrt (AA * AA + d * d - 2 * d * AA * cos (180 - z))
return HH - 6371000
.
func density a z .
while d < 10000000
delta = higher 0.001 (0.001 * d)
sum += rho height a z (d + 0.5 * delta) * delta
d += delta
.
return sum
.
func airmass a z .
return density a z / density a 0
.
numfmt 8 2
print "Angle 0 m 13700 m"
print "------------------------"
for z = 0 step 5 to 90
print z & " " & airmass 0 z & " " & airmass 13700 z
.
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
- 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
Go
package main
import (
"fmt"
"math"
)
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
)
// The density of air as a function of height above sea level.
func rho(a float64) float64 { return math.Exp(-a / 8500) }
// Converts degrees to radians
func radians(degrees float64) float64 { return degrees * math.Pi / 180 }
// a = altitude of observer
// z = zenith angle (in degrees)
// d = distance along line of sight
func height(a, z, d float64) float64 {
aa := RE + a
hh := math.Sqrt(aa*aa + d*d - 2*d*aa*math.Cos(radians(180-z)))
return hh - RE
}
// Integrates density along the line of sight.
func columnDensity(a, z float64) float64 {
sum := 0.0
d := 0.0
for d < FIN {
delta := math.Max(DD, DD*d) // adaptive step size to avoid it taking forever
sum += rho(height(a, z, d+0.5*delta)) * delta
d += delta
}
return sum
}
func airmass(a, z float64) float64 {
return columnDensity(a, z) / columnDensity(a, 0)
}
func main() {
fmt.Println("Angle 0 m 13700 m")
fmt.Println("------------------------------------")
for z := 0; z <= 90; z += 5 {
fz := float64(z)
fmt.Printf("%2d %11.8f %11.8f\n", z, airmass(0, fz), airmass(13700, fz))
}
}
- 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
Java
public class AirMass {
public static void main(String[] args) {
System.out.println("Angle 0 m 13700 m");
System.out.println("------------------------------------");
for (double z = 0; z <= 90; z+= 5) {
System.out.printf("%2.0f %11.8f %11.8f\n",
z, airmass(0.0, z), airmass(13700.0, z));
}
}
private static double rho(double a) {
// the density of air as a function of height above sea level
return Math.exp(-a / 8500.0);
}
private static double height(double a, double z, double d) {
// a = altitude of observer
// z = zenith angle (in degrees)
// d = distance along line of sight
double aa = RE + a;
double hh = Math.sqrt(aa * aa + d * d - 2.0 * d * aa * Math.cos(Math.toRadians(180 - z)));
return hh - RE;
}
private static double columnDensity(double a, double z) {
// integrates density along the line of sight
double sum = 0.0, d = 0.0;
while (d < FIN) {
// adaptive step size to avoid it taking forever
double delta = Math.max(DD * d, DD);
sum += rho(height(a, z, d + 0.5 * delta)) * delta;
d += delta;
}
return sum;
}
private static double airmass(double a, double z) {
return columnDensity(a, z) / columnDensity(a, 0.0);
}
private static final double RE = 6371000.0; // Earth radius in meters
private static final double DD = 0.001; // integrate in this fraction of the distance already covered
private static final double FIN = 10000000.0; // integrate only to a height of 10000km, effectively infinity
}
- 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
jq
Adapted from Wren
Works with gojq, the Go implementation of jq
Preliminaries
def pi: 4 * (1|atan);
def radians: . * pi / 180;
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
# Input: a number
# Output: a string with $digits fractional decimal digits, with proper rounding
def fmt($width; $digits):
. as $in
| tostring
| index(".") as $ix
| if test("[eE]") then .
elif $ix
then pow(10; $digits) as $p
| ($in * $p | round | tostring) as $s
| if test("[eE]") then $s
else ($s | index(".")) as $ix
| if $ix then $s[:$ix + 1] + $s[$ix+1: $ix+1+$digits]
else $s[:-$digits] + "." + $s[-$digits:]
end
end
else . + "." + "0" * digits
end
| lpad($width);
Physics
# constants
def RE: 6371000; # radius of earth in meters
def DD: 0.001; # integrate in this fraction of the distance already covered
def FIN: 1e7; # integrate only to a height of 10000km, effectively infinity
# The density of air as a function of height above sea level.
def rho: (-./8500) | exp;
# a = altitude of observer (in m)
# z = zenith angle (in degrees)
# d = distance along line of sight (in m)
def height($a; $z; $d):
(RE + $a) as $aa
| (($aa * $aa + $d * $d - 2 * $d * $aa * ((180-$z)|radians|cos) )|sqrt ) - RE;
# Integrates density along the line of sight.
def columnDensity($a; $z):
{ sum: 0, d: 0 }
| until (.d >= FIN;
([DD, DD * .d] | max) as $delta # adaptive step size to avoid it taking forever
| .sum = .sum + ((height($a; $z; .d + 0.5 * $delta))|rho) * $delta
| .d += $delta )
| .sum ;
def airmass(a; z): columnDensity(a; z) / columnDensity(a; 0);
"Angle 0 m 13700 m",
"------------------------------------",
( range(0; 91; 5)
| "\(lpad(2)) \(airmass(0; .)|fmt(11;8)) \(airmass(13700; .)|fmt(11;8))" )
- 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
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
- 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
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
- 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
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);
}
- 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
Python
""" Rosetta Code task: Air_mass """
from math import sqrt, cos, exp
DEG = 0.017453292519943295769236907684886127134 # degrees to radians
RE = 6371000 # Earth radius in meters
dd = 0.001 # integrate in this fraction of the distance already covered
FIN = 10000000 # integrate only to a height of 10000km, effectively infinity
def rho(a):
""" the density of air as a function of height above sea level """
return exp(-a / 8500.0)
def height(a, z, d):
"""
a = altitude of observer
z = zenith angle (in degrees)
d = distance along line of sight
"""
return sqrt((RE + a)**2 + d**2 - 2 * d * (RE + a) * cos((180 - z) * DEG)) - RE
def column_density(a, z):
""" integrates density along the line of sight """
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
return dsum
def airmass(a, z):
return column_density(a, z) / column_density(a, 0)
print('Angle 0 m 13700 m\n', '-' * 36)
for z in range(0, 91, 5):
print(f"{z: 3d} {airmass(0, z): 12.7f} {airmass(13700, z): 12.7f}")
- Output:
Angle 0 m 13700 m ------------------------------------ 0 1.0000000 1.0000000 5 1.0038096 1.0038096 10 1.0153847 1.0153848 15 1.0351774 1.0351776 20 1.0639905 1.0639909 25 1.1030594 1.1030601 30 1.1541897 1.1541908 35 1.2199808 1.2199825 40 1.3041893 1.3041919 45 1.4123417 1.4123457 50 1.5528040 1.5528102 55 1.7387592 1.7387692 60 1.9921200 1.9921366 65 2.3519974 2.3520272 70 2.8953137 2.8953729 75 3.7958235 3.7959615 80 5.5388581 5.5392811 85 10.0789622 10.0811598 90 34.3298114 34.3666656
Raku
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);
}
- 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
/*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
- 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 ─────┴─────────────────────────────────────────────────────────────
RPL
≪ → a ≪ a NEG 8500 / EXP ≫ ≫ ‘RHO’ STO ≪ 'RHO(√(aa^2+D^2-2*aa*D*COS(180-Z))-re)' EVAL ≫ ‘COLD’ STO ≪ 6371000 3 PICK OVER + → a z re aa ≪ DEG z 'Z' STO 'COLD' { D 0 1E7 } 1E-7 ∫ DROP 0 'Z' STO 'COLD' { D 0 1E7 } 1E-7 ∫ DROP / ≫ ≫ ‘AM’ STO
≪ { } 0 90 FOR z z + z AM + 5 STEP ≫
- Output:
1: { 0 1 5 1.00380686363 10 1.01537368745 15 1.03515302646 20 1.06394782383 25 1.10299392042 30 1.15409753978 35 1.21985818096 40 1.30403285254 45 1.41214767977 50 1.55256798138 55 1.73847475415 60 1.99177718552 65 2.35157928673 70 2.89478919419 75 3.79512945489 80 5.53784169364 85 10.0771111633 90 34.3235064081 }
Rust
const RE: f64 = 6371000.0; // Earth radius in meters
const DD: f64 = 0.001; // integrate in this fraction of the distance already covered
const FIN: f64 = 10000000.0; // integrate only to a height of 10000km, effectively infinity
fn rho(a: f64) -> f64 {
// the density of air as a function of height above sea level
(-a / 8500.0).exp()
}
fn height(a: f64, z: f64, d: f64) -> f64 {
// a = altitude of observer
// z = zenith angle (in degrees)
// d = distance along line of sight
let aa = RE + a;
let hh = (aa * aa + d * d - 2.0 * d * aa * (180.0 - z).to_radians().cos()).sqrt();
hh - RE
}
fn column_density(a: f64, z: f64) -> f64 {
// integrates density along the line of sight
let mut sum = 0.0;
let mut d = 0.0;
while d < FIN {
// adaptive step size to avoid it taking forever
let mut delta = DD * d;
if delta < DD {
delta = DD;
}
sum += rho(height(a, z, d + 0.5 * delta)) * delta;
d += delta;
}
sum
}
fn airmass(a: f64, z: f64) -> f64 {
column_density(a, z) / column_density(a, 0.0)
}
fn main() {
println!("Angle 0 m 13700 m");
println!("------------------------------------");
for a in (0..=90).step_by(5) {
let z = a as f64;
println!(
"{:2} {:11.8} {:11.8}",
z,
airmass(0.0, z),
airmass(13700.0, z)
);
}
}
- 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
Seed7
$ include "seed7_05.s7i";
include "float.s7i";
include "math.s7i";
const float: DEG is 0.017453292519943295769236907684886127134; #degrees to radians
const float: RE is 6371000.0; #Earth radius in meters
const float: dd is 0.001; #integrate in this fraction of the distance already covered
const float: FIN is 10000000.0; #integrate only to a height of 10000km, effectively infinity
const func float: rho (in float: a) is
#the density of air as a function of height above sea level
return exp(-a / 8500.0);
const func float: height (in float: a, in float: z, in float: d) is func
#a is altitude of observer
#z is zenith angle (in degrees)
#d is distance along line of sight
result
var float: r is 0.0;
local
var float: AA is 0.0;
var float: HH is 0.0;
begin
AA := RE + a;
HH := sqrt( AA ** 2.0 + d ** 2.0 - 2.0 * d * AA * cos((180.0 - z) * DEG) );
r := HH - RE;
end func;
const func float: columnDensity (in float: a, in float: z) is func
#integrates density along line of sight
result
var float: sum is 0.0;
local
var float: d is 0.0;
var float: delta is 0.0;
begin
while d < FIN do
delta := max(dd, dd * d); #adaptive step size to avoid taking it forever
sum +:= rho(height(a, z, d + 0.5 * delta)) * delta;
d +:= delta;
end while;
end func;
const func float: airmass (in float: a, in float: z) is
return columnDensity(a, z) / columnDensity(a, 0.0);
const proc: main is func
local
var integer: zz is 0;
var float: z is 0.0;
begin
writeln("Angle 0 m 13700 m");
writeln("------------------------------------");
for zz range 0 to 90 step 5 do
z := flt(zz);
write(z lpad 4);
write(airmass(0.0, z) digits 8 lpad 15);
writeln(airmass(13700.0, z) digits 8 lpad 17);
end for;
end func;
- Output:
Angle 0 m 13700 m ------------------------------------ 0.0 1.00000000 1.00000000 5.0 1.00380963 1.00380965 10.0 1.01538466 1.01538475 15.0 1.03517744 1.03517765 20.0 1.06399053 1.06399093 25.0 1.10305937 1.10306005 30.0 1.15418974 1.15419083 35.0 1.21998076 1.21998246 40.0 1.30418931 1.30419190 45.0 1.41234169 1.41234567 50.0 1.55280404 1.55281025 55.0 1.73875921 1.73876915 60.0 1.99212000 1.99213665 65.0 2.35199740 2.35202722 70.0 2.89531368 2.89537287 75.0 3.79582352 3.79596149 80.0 5.53885809 5.53928113 85.0 10.07896219 10.08115981 90.0 34.32981136 34.36666557
Swift
import Foundation
extension Double {
var radians: Double { self * .pi / 180 }
}
func columnDensity(_ a: Double, _ z: Double) -> Double {
func rho(_ a: Double) -> Double {
exp(-a / 8500)
}
func height(_ d: Double) -> Double {
let aa = 6_371_000 + a
let hh = aa * aa + d * d - 2 * d * aa * cos((180 - z).radians)
return hh.squareRoot() - 6_371_000
}
var sum = 0.0
var d = 0.0
while d < 1e7 {
let delta = max(0.001, 0.001 * d)
sum += rho(height(d + 0.5 * delta)) * delta
d += delta
}
return sum
}
func airMass(altitude: Double, zenith: Double) -> Double {
return columnDensity(altitude, zenith) / columnDensity(altitude, 0)
}
print("Angle 0 m 13700 m")
print("------------------------------------")
for z in stride(from: 0.0, through: 90.0, by: 5.0) {
let air = String(
format: "%2.0f %11.8f %11.8f",
z,
airMass(altitude: 0, zenith: z),
airMass(altitude: 13700, zenith: z)
)
print(air)
}
- 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
Wren
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| (-a/8500).exp }
// 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 = DD.max(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
}
- 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
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.;
]
]
- 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