atm: profile_func seems to ignore observer altitude

The following was reported via email (by Robert Miesen): [translated]

pycraf slant attenuation values seem to differ from ITU's MatLab version (see here), and appear mostly independent on the observers altitude. It could be associated with line 1336 in atm.py, where the call to profile_func is done with the heights variable, which is not accounting for the observer height.

Here's an example:

import numpy as np
from astropy import units as u
from pycraf import atm

elevation = [5, 10, 50, 90] * u.deg
altitude = [10, 100, 1000, 10000, 100000] * u.m
frequency_q = 26 * u.GHz

# profile = atm.profile_standard
# profile = atm.profile_lowlat
profile = atm.profile_midlat_summer
# profile = atm.profile_highlat_summer

attenuation = np.zeros((len(elevation) * len(altitude), 1))

print('Frequency | station height | elevation | attenuation')

for el_idx in range(len(elevation)):
    for alt_idx in range(len(altitude)):
        total_atten, refraction, tebb = atm.atten_slant_annex1(
            frequency_q, elevation[el_idx], altitude[alt_idx], profile, t_bg=2.73 * u.K)
        attenuation[alt_idx * len(elevation) + el_idx] = total_atten.value
        print('{:9.1f}'.format(frequency_q.value) +
              '{:17.1f}'.format(altitude[alt_idx].value) +
              '{:12.1f}'.format(elevation[el_idx].value) +
              '{:14.2f}'.format(total_atten.value[0]))

# Output:

# Frequency | station height | elevation | attenuation
#      26.0             10.0         5.0          5.54
#      26.0            100.0         5.0          5.54
#      26.0           1000.0         5.0          5.54
#      26.0          10000.0         5.0          5.54
#      26.0         100000.0         5.0          5.54
#      26.0             10.0        10.0          2.84
#      26.0            100.0        10.0          2.84
#      26.0           1000.0        10.0          2.84
#      26.0          10000.0        10.0          2.84
#      26.0         100000.0        10.0          2.84
#      26.0             10.0        50.0          0.65
#      26.0            100.0        50.0          0.65
#      26.0           1000.0        50.0          0.65
#      26.0          10000.0        50.0          0.65
#      26.0         100000.0        50.0          0.65
#      26.0             10.0        90.0          0.50
#      26.0            100.0        90.0          0.50
#      26.0           1000.0        90.0          0.50
#      26.0          10000.0        90.0          0.50
#      26.0         100000.0        90.0          0.50

But it should be:

Atmosphere Frequency St_Height Elevation Attenuation
1.00       26.00       0.01     5.00     5.44
1.00       26.00       0.01    10.00     2.70
1.00       26.00       0.01    50.00     0.61
1.00       26.00       0.01    90.00     0.47
1.00       26.00       0.10     5.00     5.20
1.00       26.00       0.10    10.00     2.59
1.00       26.00       0.10    50.00     0.58
1.00       26.00       0.10    90.00     0.45
1.00       26.00       1.00     5.00     3.28
1.00       26.00       1.00    10.00     1.65
1.00       26.00       1.00    50.00     0.37
1.00       26.00       1.00    90.00     0.29
1.00       26.00      10.00     5.00     0.10
1.00       26.00      10.00    10.00     0.05
1.00       26.00      10.00    50.00     0.01
1.00       26.00      10.00    90.00     0.01
1.00       26.00     100.00     5.00      NaN
1.00       26.00     100.00    10.00      NaN
1.00       26.00     100.00    50.00      NaN
1.00       26.00     100.00    90.00      NaN
(St_height in km)

I've introduced a fix (the for-loop in atm.py:1373 is now starting at the layer index that fits the observer's altitude). New output:

Frequency | station height | elevation | attenuation
     26.0             10.0         5.0          5.51
     26.0             10.0        10.0          2.83
     26.0             10.0        50.0          0.65
     26.0             10.0        90.0          0.49
     26.0            100.0         5.0          5.27
     26.0            100.0        10.0          2.70
     26.0            100.0        50.0          0.62
     26.0            100.0        90.0          0.47
     26.0           1000.0         5.0          3.34
     26.0           1000.0        10.0          1.72
     26.0           1000.0        50.0          0.39
     26.0           1000.0        90.0          0.30
     26.0          10000.0         5.0          0.10
     26.0          10000.0        10.0          0.05
     26.0          10000.0        50.0          0.01
     26.0          10000.0        90.0          0.01
     26.0         100000.0         5.0          0.00
     26.0         100000.0        10.0          0.00
     26.0         100000.0        50.0          0.00
     26.0         100000.0        90.0          0.00

There is still a slight discrepancy between the MatLab values and pycraf results. Not sure, where this originates from.

bwinkel / pycraf

atm: profile_func seems to ignore observer altitude #1