Air mass: Difference between revisions

=={{header|AWK}}==

<lang AWK>

# syntax: GAWK -f AIR_MASS.AWK

# converted from FreeBASIC

BEGIN {

dd = 0.001 # integrate in this fraction of the distance already covered

DEG = 0.017453292519943295769236907684886127134 # degrees to radians

RE = 6371000 # Earth radius in meters

print("Angle 0 m 13700 m")

for (z=0; z<=90; z+=5) {

printf("%5d %12.8f %12.8f\n",z,am_airmass(0,z),am_airmass(13700,z))

}

exit(0)

}

function am_airmass(a,z) {

return am_column_density(a,z) / am_column_density(a,0)

}

function am_column_density(a,z, d,delta,sum) { # integrates density along the line of sight

while (d < 10000000) { # integrate only to a height of 10000km, effectively infinity

delta = max(dd,(dd)*d) # adaptive step size to avoid it taking forever

sum += am_rho(am_height(a,z,d+0.5*delta))*delta

d += delta

}

return(sum)

}

function am_height(a,z,d, aa,hh) {

# a - altitude of observer

# z - zenith angle in degrees

# d - distance along line of sight

aa = RE + a

hh = sqrt(aa^2 + d^2 - 2*d*aa*cos((180-z)*DEG))

return(hh-RE)

}

function am_rho(a) { # density of air as a function of height above sea level

return exp(-a/8500.0)

}

function max(x,y) { return((x > y) ? x : y) }

</lang>

{{out}}

<pre>

Angle 0 m 13700 m

0 1.00000000 1.00000000

5 1.00380963 1.00380965

10 1.01538466 1.01538475

15 1.03517744 1.03517765

20 1.06399053 1.06399093

25 1.10305937 1.10306005

30 1.15418974 1.15419083

35 1.21998076 1.21998246

40 1.30418931 1.30419190

45 1.41234169 1.41234567

50 1.55280404 1.55281025

55 1.73875921 1.73876915

60 1.99212000 1.99213665

65 2.35199740 2.35202722

70 2.89531368 2.89537287

75 3.79582352 3.79596149

80 5.53885809 5.53928113

85 10.07896219 10.08115981

90 34.32981136 34.36666557

</pre>