I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

# Air mass

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 ${\displaystyle \rho }$(h(a,z,x)) where ${\displaystyle \rho }$(h) is the atmospheric density for a given height above sea level, and h(a,z,x) is the height above sea level for a point at distance x along the line of sight. Determining this last function requires some trigonometry.

For this task you can assume:

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

•   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.017453292519943295769236907684886127134V RE = 6371000V dd = 0.001V 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


## AWK

 # syntax: GAWK -f AIR_MASS.AWK# converted from FreeBASICBEGIN {    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, LIMRE  = 6371000  #Earth radius in metersdd  = 0.001    #integrate in this fraction of the distance already coveredLIM = 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 zend function max(a, b)    if a > b then return a else return bend 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 - REend 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 sumend 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 - REend 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 sumend 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 = bEND 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-reEND 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 = sumEND FUNCTION FUNCTION airmass(a, z)    LET airmass = columndensity(a, z)/columndensity(a, 0)END FUNCTION LET deg = .0174532925199433       !degrees TO radiansLET re = 6371000                  !Earth radius in metersLET dd = .001                     !integrate in this fraction of the distance already coveredLET lim = 10000000                !integrate only TO a height of 10000km, effectively infinityPRINT "Angle     0 m              13700 m"PRINT "------------------------------------"FOR z = 0 TO 90 STEP 5    PRINT  USING "##      ##.########      ##.########": z, airmass(0, z), airmass(13700, z)NEXT zEND

## 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


## Factor

Translation of: FreeBASIC
Works with: Factor version 0.99 2021-02-05
USING: formatting io kernel math math.functions math.ordermath.ranges math.trig sequences ; CONSTANT: RE 6,371,000     ! Earth's radius in metersCONSTANT: dd 0.001         ! integrate in this fraction of the distance already coveredCONSTANT: 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"------------------------------------" print0 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 radiansfunc radians(degrees float64) float64 { return degrees * math.Pi / 180 } // a = altitude of observer// z = zenith angle (in degrees)// d = distance along line of sightfunc 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

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 roundingdef 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 # constantsdef RE: 6371000; # radius of earth in metersdef DD: 0.001; # integrate in this fraction of the distance already covereddef 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 radiansconst RE = 6371000 # Earth radius in metersconst dd = 0.001 # integrate in this fraction of the distance already coveredconst 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 dsumend 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.0while 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 radiansuse constant RE => 6371000; # Earth radius in metersuse constant dd => 0.001; # integrate in this fraction of the distance already covereduse constant FIN => 10000000; # integrate only to a height of 10000km, effectively infinity # Density of air as a function of height above sea levelsub 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 sightsub 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 radiansRE = 6371000 # Earth radius in metersdd = 0.001 # integrate in this fraction of the distance already coveredFIN = 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 radiansconstant RE = 6371000; # Earth radius in metersconstant dd = 0.001; # integrate in this fraction of the distance already coveredconstant FIN = 10000000; # integrate only to a height of 10000km, effectively infinity #| Density of air as a function of height above sea levelsub 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 sightsub 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 pirho: procedure; parse arg a; return exp(-a / 8500)r2r: return arg(1) // (pi() * 2) /*normalize radians ──► a unit circle. */e: e= 2.718281828459045235360287471352662497757247093699959574966967627724; return ecommas: 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 ─────┴─────────────────────────────────────────────────────────────  ## Rust Translation of: FreeBASIC const RE: f64 = 6371000.0; // Earth radius in metersconst DD: f64 = 0.001; // integrate in this fraction of the distance already coveredconst 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 radiansconst float: RE is 6371000.0;   #Earth radius in metersconst float: dd is 0.001;       #integrate in this fraction of the distance already coveredconst 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 Mathimport "/fmt" for Fmt // constantsvar RE  = 6371000  // radius of earth in metersvar DD  = 0.001    // integrate in this fraction of the distance already coveredvar 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 sightvar 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 = 0while (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 radiansdefine RE = 6371000.;   \Earth radius in metersdefine DD = 0.001;      \integrate in this fraction of the distance already covereddefine 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 sightreal 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 sightreal 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