#DEFINITIONS:  -*-sh-*-
#
# Demonstration and test of simple ARTS 1D, 2D and 3D clear sky calculations.
#
# Observations from a satellite is treated, with three viewing directions:
#   1: Cold space (ie. above the model atmosphere)
#   2: Limb sounding
#   3: Downward observation.
#
# For the test sequence (ie. "make check") this is a first test on that the
# full chain around yCalc is working. Calculation of optical depth as an
# auxilary variable is included and tested.
#
# Author: Patrick Eriksson


Arts2 {

INCLUDE "general/general.arts"
INCLUDE "general/continua.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/planet_earth.arts"

# Agenda for scalar gas absorption calculation
Copy(abs_xsec_agenda, abs_xsec_agenda__noCIA)

# (standard) emission calculation
Copy( iy_main_agenda, iy_main_agenda__Emission )

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# standard surface agenda (i.e., make use of surface_rtprop_agenda)
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )

# on-the-fly absorption
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )


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

# No jacobian calculation
#
jacobianOff

# Clearsky = No scattering
#
cloudboxOff


# Read a line file and a matching small frequency grid
# ---
ReadXML( abs_lines, "abs_lines.xml" )
VectorNLinSpace( f_grid, 5, 320e9, 322e9 )


# A pressure grid rougly matching 0 to 80 km, in steps of 2 km.
# ---
VectorNLogSpace( p_grid, 41, 1000e2, 1 )


# Definition of species
# ---
abs_speciesSet( species=
            ["H2O-SelfContStandardType, H2O-ForeignContStandardType, H2O", 
             "N2-SelfContStandardType",
             "O3"] )


# Sort the line file according to species
# ---
abs_lines_per_speciesCreateFromLines

# Atmospheric scenario
# ---
AtmRawRead( basename = "testdata/tropical" )


# Weakly reflecting surface
# ---
VectorSetConstant( surface_scalar_reflectivity, 1, 0.8 )
Copy( surface_rtprop_agenda,
      surface_rtprop_agenda__Specular_NoPol_ReflFix_SurfTFromt_surface )

# No sensor properties
# ---
sensorOff


# We select here to use Rayleigh-Jean brightness temperatures
# ---
StringSet( iy_unit, "RJBT" )


# Extract radiative background and optical depth as auxiliary variables
# ---
ArrayOfStringSet( iy_aux_vars, [ "Optical depth", "Radiative background" ] )
# Create vector container for the optical depth
VectorCreate( odepth )



#########################################################################
# 1D
#########################################################################

# Atmosphere and surface
# ---
AtmosphereSet1D
AtmFieldsCalc
Extract( z_surface, z_field, 0 )
Extract( t_surface, t_field, 0 )

# Definition of sensor position and LOS
# ---
MatrixSetConstant( sensor_pos, 3, 1, 600e3 )
MatrixSet( sensor_los, [ 95; 113; 135] )

# Perform RT calculations
# ---
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc
yCalc

# OK?
# ---
Extract( odepth, y_aux, 0 )
#WriteXML( "ascii", y, "yREFERENCE_1D.xml" )
#WriteXML( "ascii", odepth, "y_auxREFERENCE_1D.xml" )
VectorCreate( yREFERENCE )
ReadXML( yREFERENCE, "yREFERENCE_1D.xml" )
Compare( y, yREFERENCE, 1e-2 )
ReadXML( yREFERENCE, "y_auxREFERENCE_1D.xml" )
Compare( odepth, yREFERENCE, 1e-3 )



#########################################################################
# 2D
#########################################################################

# Atmosphere and surface
# ---
AtmosphereSet2D
VectorLinSpace( lat_grid, -45, 45, 1 )
AtmFieldsCalcExpand1D
Extract( z_surface, z_field, 0 )
Extract( t_surface, t_field, 0 )
refellipsoidEarth( refellipsoid, "WGS84" )

# Definition of sensor position and LOS
# ---
MatrixCreate( zeros )
MatrixSetConstant( zeros, 3, 1, 0 )
Append( sensor_pos, zeros, "trailing" )

# Perform RT calculations
# ---
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc
yCalc

# OK?
# ---
Extract( odepth, y_aux, 0 )
#WriteXML( "ascii", y, "yREFERENCE_2D.xml" )
#WriteXML( "ascii", odepth, "y_auxREFERENCE_2D.xml" )
ReadXML( yREFERENCE, "yREFERENCE_2D.xml" )
Compare( y, yREFERENCE, 1e-2 )
ReadXML( yREFERENCE, "y_auxREFERENCE_2D.xml" )
Compare( odepth, yREFERENCE, 1e-3 )





#########################################################################
# 3D
#########################################################################

# Atmosphere and surface
# ---
AtmosphereSet3D
VectorLinSpace( lon_grid, -45, 45, 1 )
AtmFieldsCalcExpand1D
Extract( z_surface, z_field, 0 )
Extract( t_surface, t_field, 0 )

# Definition of sensor position and LOS
# ---
Append( sensor_pos, zeros, "trailing" )
Append( sensor_los, zeros, "trailing" )

# Perform RT calculations
# ---
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc
yCalc

# OK?
# ---
Extract( odepth, y_aux, 0 )
#WriteXML( "ascii", y, "yREFERENCE_3D.xml" )
#WriteXML( "ascii", odepth, "y_auxREFERENCE_3D.xml" )
ReadXML( yREFERENCE, "yREFERENCE_3D.xml" )
Compare( y, yREFERENCE, 1e-2 )
ReadXML( yREFERENCE, "y_auxREFERENCE_3D.xml" )
Compare( odepth, yREFERENCE, 1e-3 )

} # End of Main