Air mass

From Rosetta Code
Air mass is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

Translation of: Python
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

Translation of: C
Works with: Ada 2012
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

Translation of: 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

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

Translation of: FreeBASIC
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

Translation of: FreeBASIC
#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

Works with: Delphi version 6.0
Translation of: GO
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

Translation of: FreeBASIC
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

Translation of: FreeBASIC
Works with: Factor version 0.99 2021-02-05
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

Translation of: FreeBASIC
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

Translation of: FreeBASIC
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: jq

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

Translation of: FreeBASIC
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

Translation of: Wren
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

Translation of: Raku
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

Translation of: FreeBASIC
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² +  - 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

Translation of: FreeBASIC
/*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

Translation of: 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
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

Translation of: FreeBASIC
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

Translation of: FreeBASIC
$ 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

Translation of: Go
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

Translation of: FreeBASIC
Library: Wren-math
Library: Wren-fmt
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

Translation of: 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
 
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