Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Xia Mengqi]

Hi Freddy,


Thank you for the write-up. If I understand correctly you have the delta 
functions in the flat surface case to specify the single reflection direction. 
If we consider the light coming from the star only has one single direction, 
which I think it's a reasonable assumption for the sun; then if the sun is at 0 
latitude and 0 longtitude and the sensor looking down then the ratio should be 
ratio of two rhos in your equation, which is pi I think? It seems that in ARTS 
the sun is modeled as an area light source, but the solid angle it subtends is 
pretty small. So even if we integrate the directions I would think the ratio is 
still close to pi. Maybe I'm missing something here.


Thanks,

Mengqi


From: Manfred Brath 
Sent: Thursday, October 20, 2022 9:57:46 AM
To: Xia Mengqi; arts_users.mi@lists.uni-hamburg.de
Subject: Re: Fwd: [arts-users] Direct radiation is subtracted when it should not


Hello Mengqi,


I think, I understand your problem now. Please have a look in the attached pdf 
(I hope my handwriting is readable).


Cheers,

Freddy


Am 19.10.2022 um 20:49 schrieb Xia Mengqi:

Hi Freddy,


Thanks for getting back to me. My question is about the result of using a 
specular reflective surface versus using a Lambertian surface. I was thinking 
the ratio should be pi (the specular one is larger) because when there is no 
scattering the only difference is the BRDF evaluation. However, I observed the 
Lambertian result is much smaller, more than 1e4 times smaller than the 
specular result. I tried "Test_iySurfaceLambertian.py" and it is also the case 
there when I change the Lambertian surface to specular.


Thanks,

Mengqi


From: Manfred Brath 

Sent: Wednesday, October 19, 2022 5:19:20 PM
To: Xia Mengqi; 
arts_users.mi@lists.uni-hamburg.de
Subject: Re: Fwd: [arts-users] Direct radiation is subtracted when it should not


Hi Mengqi,


I do not understand your problem from the information of your email, but maybe 
you can have a look in "Test_iySurfaceLambertian.py", which you can find in 
your arts-folder:


controlfiles-python/artscomponents/surface/


Another thing, how do you calculate your pi-ratio? When I use the above 
mentioned test script and set the sun at 0 lat and 0 lon and divide the top of 
the atmosphere (TOA) solar radiation F by pi I got the same value as from yCalc.


F=ws.stars.value[0].spectrum.value*scale_factor

scale_factor=sun_radius**2/(sun_radius**2+sun2TOA_distance**2)


The scale factor is needed, because the star irradiance spectrum in ARTS is 
defined at the surface of the sun and not at TOA.


Cheers,

Freddy


Am 18.10.22 um 12:56 schrieb Xia Mengqi:

Thank you Freddy! This makes sense.


I tried with flat reflective surface and lambertian surface with absorption 
only. I extended the 1D atmosphere to 3D so I can provide surface temperature 
(set it to a really small value) but in theory since there is no scattering the 
path that contributes is unchanged. I found that the flat reflective surface 
result is correct but the lambertian one is much smaller than the expected pi 
ratio. I tried to print ppath and I found los changing quite a lot. I am 
wondering if this is the expected behavior and maybe there is something about 
the 3D setup I did not understand properly. I copied the main part of the code 
below.


Thanks!

Mengqi


@arts_agenda
def propmat_clearsky_agenda(ws):
ws.Ignore(ws.rtp_mag)
ws.Ignore(ws.rtp_los)
ws.propmat_clearskyInit()
ws.propmat_clearskyAddConts()
ws.propmat_clearskyAddLines()

@arts_agenda
def gas_scattering_agenda(ws):
ws.Ignore(ws.rtp_vmr)
ws.gas_scattering_coefXsecConst(ConstXsec=4.65e-31)
ws.gas_scattering_matIsotropic()

# surface scattering agenda
# lambertian
@arts_agenda
def iy_surface_agenda(ws):
ws.iySurfaceInit()
ws.Ignore(ws.dsurface_rmatrix_dx)
ws.Ignore(ws.dsurface_emission_dx)

ws.iySurfaceLambertian()
ws.iySurfaceLambertianDirect()

# flat reflective surface
# @arts_agenda
# def iy_surface_agenda(ws):
# ws.iySurfaceInit()
# ws.iySurfaceFlatReflectivity()
# ws.iySurfaceFlatReflectivityDirect()

# generate atmosphere data
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)


# =
# open workspace
# =

ws = Workspace()
ws.verbositySetScreen(level=2)

# =
# generate atmosphere data
# =
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_pat

Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Xia Mengqi]

Hi Freddy,


Thanks for getting back to me. My question is about the result of using a 
specular reflective surface versus using a Lambertian surface. I was thinking 
the ratio should be pi (the specular one is larger) because when there is no 
scattering the only difference is the BRDF evaluation. However, I observed the 
Lambertian result is much smaller, more than 1e4 times smaller than the 
specular result. I tried "Test_iySurfaceLambertian.py" and it is also the case 
there when I change the Lambertian surface to specular.


Thanks,

Mengqi


From: Manfred Brath 
Sent: Wednesday, October 19, 2022 5:19:20 PM
To: Xia Mengqi; arts_users.mi@lists.uni-hamburg.de
Subject: Re: Fwd: [arts-users] Direct radiation is subtracted when it should not


Hi Mengqi,


I do not understand your problem from the information of your email, but maybe 
you can have a look in "Test_iySurfaceLambertian.py", which you can find in 
your arts-folder:


controlfiles-python/artscomponents/surface/


Another thing, how do you calculate your pi-ratio? When I use the above 
mentioned test script and set the sun at 0 lat and 0 lon and divide the top of 
the atmosphere (TOA) solar radiation F by pi I got the same value as from yCalc.


F=ws.stars.value[0].spectrum.value*scale_factor

scale_factor=sun_radius**2/(sun_radius**2+sun2TOA_distance**2)


The scale factor is needed, because the star irradiance spectrum in ARTS is 
defined at the surface of the sun and not at TOA.


Cheers,

Freddy


Am 18.10.22 um 12:56 schrieb Xia Mengqi:

Thank you Freddy! This makes sense.


I tried with flat reflective surface and lambertian surface with absorption 
only. I extended the 1D atmosphere to 3D so I can provide surface temperature 
(set it to a really small value) but in theory since there is no scattering the 
path that contributes is unchanged. I found that the flat reflective surface 
result is correct but the lambertian one is much smaller than the expected pi 
ratio. I tried to print ppath and I found los changing quite a lot. I am 
wondering if this is the expected behavior and maybe there is something about 
the 3D setup I did not understand properly. I copied the main part of the code 
below.


Thanks!

Mengqi


@arts_agenda
def propmat_clearsky_agenda(ws):
ws.Ignore(ws.rtp_mag)
ws.Ignore(ws.rtp_los)
ws.propmat_clearskyInit()
ws.propmat_clearskyAddConts()
ws.propmat_clearskyAddLines()

@arts_agenda
def gas_scattering_agenda(ws):
ws.Ignore(ws.rtp_vmr)
ws.gas_scattering_coefXsecConst(ConstXsec=4.65e-31)
ws.gas_scattering_matIsotropic()

# surface scattering agenda
# lambertian
@arts_agenda
def iy_surface_agenda(ws):
ws.iySurfaceInit()
ws.Ignore(ws.dsurface_rmatrix_dx)
ws.Ignore(ws.dsurface_emission_dx)

ws.iySurfaceLambertian()
ws.iySurfaceLambertianDirect()

# flat reflective surface
# @arts_agenda
# def iy_surface_agenda(ws):
# ws.iySurfaceInit()
# ws.iySurfaceFlatReflectivity()
# ws.iySurfaceFlatReflectivityDirect()

# generate atmosphere data
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)


# =
# open workspace
# =

ws = Workspace()
ws.verbositySetScreen(level=2)

# =
# generate atmosphere data
# =
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)

ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)

# =
# select/define agendas
# =

ws.LegacyContinuaInit()
ws.PlanetSet(option="Earth")

# cosmic background radiation
ws.iy_space_agendaSet( option="CosmicBackground" )

# sensor-only path
ws.ppath_agendaSet( option="FollowSensorLosPath" )

# no refraction
ws.ppath_step_agendaSet( option="GeometricPath" )

# main agenda
ws.iy_main_agendaSet( option="Clearsky")

# water agenda
ws.water_p_eq_agendaSet()

# surface agenda
ws.iy_surface_agenda = iy_surface_agenda

ws.ArrayOfStringSet( ws.iy_aux_vars,
[ "Optical depth",
"Radiative background"
] )

ws.propmat_clearsky_agenda=propmat_clearsky_agenda

# gas scattering agenda
ws.gas_scattering_agenda = gas_scattering_agenda

ws.NumericSet( ws.ppath_lmax, 1e10)

# =
# basic conditions
# ===

Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Manfred Brath]

Hi Mengqi,


I do not understand your problem from the information of your email, but 
maybe you can have a look in "Test_iySurfaceLambertian.py", which you 
can find in your arts-folder:



controlfiles-python/artscomponents/surface/


Another thing, how do you calculate your pi-ratio? When I use the above 
mentioned test script and set the sun at 0 lat and 0 lon and divide the 
top of the atmosphere (TOA) solar radiation F by pi I got the same value 
as from yCalc.


    F=ws.stars.value[0].spectrum.value*scale_factor

scale_factor=sun_radius**2/(sun_radius**2+sun2TOA_distance**2)


    The scale factor is needed, because the star irradiance spectrum in 
ARTS is defined at the surface of the sun and not at TOA.



Cheers,

Freddy


Am 18.10.22 um 12:56 schrieb Xia Mengqi:


Thank you Freddy! This makes sense.


I tried with flat reflective surface and lambertian surface with 
absorption only. I extended the 1D atmosphere to 3D so I can provide 
surface temperature (set it to a really small value) but in theory 
since there is no scattering the path that contributes is unchanged. I 
found that the flat reflective surface result is correct but the 
lambertian one is much smaller than the expected pi ratio. I tried to 
print ppath and I found los changing quite a lot. I am wondering if 
this is the expected behavior and maybe there is something about the 
3D setup I did not understand properly. I copied the main part of the 
code below.



Thanks!

Mengqi



@arts_agenda
def propmat_clearsky_agenda(ws):
    ws.Ignore(ws.rtp_mag)
    ws.Ignore(ws.rtp_los)
    ws.propmat_clearskyInit()
    ws.propmat_clearskyAddConts()
    ws.propmat_clearskyAddLines()

@arts_agenda
def gas_scattering_agenda(ws):
    ws.Ignore(ws.rtp_vmr)
    ws.gas_scattering_coefXsecConst(ConstXsec=4.65e-31)
    ws.gas_scattering_matIsotropic()

# surface scattering agenda
# lambertian
@arts_agenda
def iy_surface_agenda(ws):
    ws.iySurfaceInit()
    ws.Ignore(ws.dsurface_rmatrix_dx)
    ws.Ignore(ws.dsurface_emission_dx)

    ws.iySurfaceLambertian()
    ws.iySurfaceLambertianDirect()

# flat reflective surface
# @arts_agenda
# def iy_surface_agenda(ws):
#     ws.iySurfaceInit()
#     ws.iySurfaceFlatReflectivity()
#     ws.iySurfaceFlatReflectivityDirect()

# generate atmosphere data
dataset_path = 
'/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'

save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)


# 
=

# open workspace
# 
=


ws = Workspace()
ws.verbositySetScreen(level=2)

# 
=

# generate atmosphere data
# 
=
dataset_path = 
'/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'

save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)
ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)


# 
=

# select/define agendas
# 
=


ws.LegacyContinuaInit()
ws.PlanetSet(option="Earth")

# cosmic background radiation
ws.iy_space_agendaSet( option="CosmicBackground" )

# sensor-only path
ws.ppath_agendaSet( option="FollowSensorLosPath" )

# no refraction
ws.ppath_step_agendaSet( option="GeometricPath" )

# main agenda
ws.iy_main_agendaSet( option="Clearsky")

# water agenda
ws.water_p_eq_agendaSet()

# surface agenda
ws.iy_surface_agenda = iy_surface_agenda

ws.ArrayOfStringSet( ws.iy_aux_vars,
[ "Optical depth",
"Radiative background"
] )

ws.propmat_clearsky_agenda=propmat_clearsky_agenda

# gas scattering agenda
ws.gas_scattering_agenda = gas_scattering_agenda

ws.NumericSet( ws.ppath_lmax, 1e10)

# 
=

# basic conditions
# 
=

# Postion and line-of-sight of sensor
sensor_pos = np.array([[600e+3, 0., 0.]])
# Sensor looking direction in zenith angle (0 = upwards, 180 = 
downward) and

# azimuth angle ( 0 = North, 90 = east)
sensor_los = np.array([[180, 0]])
ws.sensor_pos = sensor_pos
ws.sensor_los = sensor_los
ws.VectorSet(ws.rte_pos2, [])

# define environment
# 
=

# Number of Stokes components to be computed
ws.IndexSet(ws.stokes_dim, 1)

# Read the spectroscopic line data from the ARTS catalogue and
# create the workspace variable `lines'.
ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)

ws.abs

Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Xia Mengqi]

Thank you Freddy! This makes sense.


I tried with flat reflective surface and lambertian surface with absorption 
only. I extended the 1D atmosphere to 3D so I can provide surface temperature 
(set it to a really small value) but in theory since there is no scattering the 
path that contributes is unchanged. I found that the flat reflective surface 
result is correct but the lambertian one is much smaller than the expected pi 
ratio. I tried to print ppath and I found los changing quite a lot. I am 
wondering if this is the expected behavior and maybe there is something about 
the 3D setup I did not understand properly. I copied the main part of the code 
below.


Thanks!

Mengqi


@arts_agenda
def propmat_clearsky_agenda(ws):
ws.Ignore(ws.rtp_mag)
ws.Ignore(ws.rtp_los)
ws.propmat_clearskyInit()
ws.propmat_clearskyAddConts()
ws.propmat_clearskyAddLines()

@arts_agenda
def gas_scattering_agenda(ws):
ws.Ignore(ws.rtp_vmr)
ws.gas_scattering_coefXsecConst(ConstXsec=4.65e-31)
ws.gas_scattering_matIsotropic()

# surface scattering agenda
# lambertian
@arts_agenda
def iy_surface_agenda(ws):
ws.iySurfaceInit()
ws.Ignore(ws.dsurface_rmatrix_dx)
ws.Ignore(ws.dsurface_emission_dx)

ws.iySurfaceLambertian()
ws.iySurfaceLambertianDirect()

# flat reflective surface
# @arts_agenda
# def iy_surface_agenda(ws):
# ws.iySurfaceInit()
# ws.iySurfaceFlatReflectivity()
# ws.iySurfaceFlatReflectivityDirect()

# generate atmosphere data
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)


# =
# open workspace
# =

ws = Workspace()
ws.verbositySetScreen(level=2)

# =
# generate atmosphere data
# =
dataset_path = '/home/mandy/Github/arts/build_new/afgl_1986-us_standard.nc'
save_path = '/home/mandy/Github/arts/controlfiles/testdata/'
data = generate_atmos_arts(dataset_path, save_path)

ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)

# =
# select/define agendas
# =

ws.LegacyContinuaInit()
ws.PlanetSet(option="Earth")

# cosmic background radiation
ws.iy_space_agendaSet( option="CosmicBackground" )

# sensor-only path
ws.ppath_agendaSet( option="FollowSensorLosPath" )

# no refraction
ws.ppath_step_agendaSet( option="GeometricPath" )

# main agenda
ws.iy_main_agendaSet( option="Clearsky")

# water agenda
ws.water_p_eq_agendaSet()

# surface agenda
ws.iy_surface_agenda = iy_surface_agenda

ws.ArrayOfStringSet( ws.iy_aux_vars,
[ "Optical depth",
"Radiative background"
] )

ws.propmat_clearsky_agenda=propmat_clearsky_agenda

# gas scattering agenda
ws.gas_scattering_agenda = gas_scattering_agenda

ws.NumericSet( ws.ppath_lmax, 1e10)

# =
# basic conditions
# =
# Postion and line-of-sight of sensor
sensor_pos = np.array([[600e+3, 0., 0.]])
# Sensor looking direction in zenith angle (0 = upwards, 180 = downward) and
# azimuth angle ( 0 = North, 90 = east)
sensor_los = np.array([[180, 0]])
ws.sensor_pos = sensor_pos
ws.sensor_los = sensor_los
ws.VectorSet(ws.rte_pos2, [])

# define environment
# =
# Number of Stokes components to be computed
ws.IndexSet(ws.stokes_dim, 1)

# Read the spectroscopic line data from the ARTS catalogue and
# create the workspace variable `lines'.
ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)
ws.abs_linesNormalization(ws.abs_lines, "VVH")

# Frequency grid
c = 299792458
# Set frequency
#wavelengths = np.linspace(1250e-9, 1177e-9, 74)
#wavelengths = np.linspace(700e-9, 400e-9, 74)
wavelengths = np.linspace(1178e-9, 1177e-9, 2)
f_grid = c / wavelengths
ws.f_grid = f_grid


# set a simple blackbody sun
ws.Touch(ws.stars)
ws.starsAddSingleBlackbody(distance=1.495978707e11, latitude=0., longitude=0.)

ws.Print(ws.stars, 2)


# Reference ellipsoid
ws.refellipsoidEarth(ws.refellipsoid, "Sphere")

# A pressure grid rougly matching 0 to 80 km, in steps of 2 km.
#ws.p_grid = data.p.to_numpy()
ws.p_grid = np.array([1.01040e+05, 2.54000e-03]) # 0 and 120km

# Atmospheric dimensionality and lat/lon grids
nlat = 3
nlon = 5
ws.VectorNLinSpace(ws.lat_grid, nlat, -90., 90.)
ws.VectorNLinSpace(ws.lon_grid, nlon,

Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Manfred Brath]

Hello Mengqi,


Am 17.10.22 um 16:33 schrieb Xia Mengqi:


Hi Freddy,


Thank you so much and this is very helpful! Just to make sure I fully 
understand how to use it in the code -- According to ARTS output, it 
seems I need to provide the skin temperature as well. I noticed that 
in one provided example it has "


ws.ArrayOfStringSet(ws.surface_props_names, ["Skin temperature"])
    ws.Tensor3SetConstant(ws.surface_props_data, 1, nlat, nlon, 
ws.t_field.value[0, 0, 0])"


I'm wondering if there is a simple way to do this for 1D.



so far you cannot use a 1d atmosphere for simulation with a direct 
source like the sun because the atmosphere dimension and (radiative 
transfer) geometry are coupled in ARTS. Since you have a zenith and an 
azimuth dependency for simulation with a direct source, you have to use 
a 3d atmosphere in ARTS as simulation of 1d atmospheres have only zenith 
dependency in ARTS.




Also, is it possible to just not include surface emission?
No, but you can set the surface temperature to a small value greter than 
0 K. This would have a similar effect.



Cheers,

Freddy


--
--
Dr. Manfred Brath
Radiation and Remote Sensing
Meteorological Institute
Universität Hamburg
Bundesstraße 55
D-20146 Hamburg
Room 1535

Tel: +49 40 42838-8786

Re: [arts-users] Fwd: Direct radiation is subtracted when it should not

[Xia Mengqi]

Hi Freddy,


Thank you so much and this is very helpful! Just to make sure I fully 
understand how to use it in the code -- According to ARTS output, it seems I 
need to provide the skin temperature as well. I noticed that in one provided 
example it has "

ws.ArrayOfStringSet(ws.surface_props_names, ["Skin temperature"])
ws.Tensor3SetConstant(ws.surface_props_data, 1, nlat, nlon, 
ws.t_field.value[0, 0, 0])"

I'm wondering if there is a simple way to do this for 1D.
Also, is it possible to just not include surface emission? I'm getting a 
different result using the new setup (a correct surface and simulation in 3D) 
and I'm trying to understand where the difference comes from.

Thanks,
Mengqi


From: Manfred Brath 
Sent: Monday, October 17, 2022 2:38:24 PM
To: Xia Mengqi
Cc: Oliver Lemke
Subject: Re: Fwd: [arts-users] Direct radiation is subtracted when it should not


Hello Xia,

since the star related methods are still under development, some documentation 
is still missing. Therefore, you chose unfortunately the wrong surface agenda. 
The iySurfaceRtpropAgenda calculates only the surface reflection/scattering 
from the diffusive radiation field. No direct radiation is intended to be 
considered by this method.

If stars are present you should use following surface methods:

for a Lambertian surface with predefined reflectivity

  *   iySurfaceLambertian
  *   iySurfaceLambertianDirect

for a flat surface with predefined reflectivity

  *   iySurfaceFlatReflectivity
  *   iySurfaceFlatReflectivityDirect

for a flat surface (Fresnel reflection) with defined refractive index

  *   iySurfaceFlatRefractiveIndex
iySurfaceFlatRefractiveIndexDirect


Here is a example for a surface agenda with a Lambertian surface

@arts_agenda
def iy_surface_agenda(ws):
ws.iySurfaceInit()
ws.Ignore(ws.dsurface_rmatrix_dx)
ws.Ignore(ws.dsurface_emission_dx)
ws.iySurfaceLambertian()
ws.iySurfaceLambertianDirect()


and for a flat surface with Fresnel reflection

@arts_agenda
def iy_surface_agenda(ws):
ws.iySurfaceInit()
ws.iySurfaceFlatRefractiveIndex()
ws.iySurfaceFlatRefractiveIndexDirect()


I hope this could help you.

Cheers,

Freddy

Am 17.10.22 um 13:06 schrieb Lemke, Oliver:


Begin forwarded message:

From: Xia Mengqi mailto:mengqi@epfl.ch>>
Subject: [arts-users] Direct radiation is subtracted when it should not
Date: 16 October 2022 at 14:05:16 CEST
To: 
"arts_users.mi@lists.uni-hamburg.de" 
mailto:arts_users.mi@lists.uni-hamburg.de>>

Hi,

I'm new to ARTS and I set up a simple 1D test where there is sun, absorption 
only gas and a specular surface with reflectivity 1. The sun and sensor is both 
at 0 latitude and 0 longitude, with the sensor looking down (los=180). For 
simplicity I force it to be just 1 layer by providing the pressures on the top 
of the atmosphere and at 0km. I noticed that the radiance output is almost 
zero, which is unexpected. However, if I comment out the operation that 
subtracts the direct radiation (iy-=iy_aux[0]) in iySurfaceRtpropAgenda the 
result looks correct (and validated against libRadtran). I'm puzzled why in 
this case the subtraction is needed; or maybe I made simple mistakes in the set 
up. It would be great if you could help me with this problem. Thank you!


The main part of the code is copied here:

ws = Workspace(verbosity=2)
ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)

ws.LegacyContinuaInit()
ws.water_p_eq_agendaSet()
ws.gas_scattering_agendaSet()
ws.PlanetSet(option="Earth")
ws.iy_main_agendaSet( option="Clearsky" )
ws.iy_space_agendaSet( option="CosmicBackground" )
ws.iy_surface_agendaSet()

ws.ArrayOfStringSet( ws.iy_aux_vars, [ "Optical depth", "Radiative background"] 
)

ws.propmat_clearsky_agenda=propmat_clearsky_agenda

ws.NumericSet( ws.ppath_lmax, 1e10)
ws.ppath_step_agendaSet( option="GeometricPath" )
ws.ppath_agendaSet( option="FollowSensorLosPath" )

# define environment
# =
# Number of Stokes components to be computed
ws.IndexSet(ws.stokes_dim, 1)

# Read the spectroscopic line data from the ARTS catalogue and
# create the workspace variable `lines'.
ws.ReadHITRAN(filename='/home/mandy/Github/MiAtmosphere/HITRAN/ALL.par', 
hitran_type="Online", abs_lines=ws.abs_lines)
ws.abs_linesNormalization(ws.abs_lines, "VVH")

# Frequency grid
c = 299792458
wavelengths = np.linspace(1178e-9, 1177e-9, 2)
f_grid = c / wavelengths
ws.f_grid = f_grid


# set a simple blackbody sun
ws.starsAddSingleBlackbody(latitude=sun_pos[0], longitude=sun_pos[1])
ws.Print(ws.stars, 2)

# A pressure grid rougly matching 0 to 80 km, in steps of 2 km.
#ws.p_grid = data.p.to_numpy()
ws.p_grid = np.array([1.01040e+05, 2.54000e-03]) # 0 and 120km

# Definition of species
ws.abs_speciesSet(species=
  ["H2O", "O