#DEFINITIONS:  -*-sh-*-
#
# Demonstration and test of a single radio link calculations, using iyCalc.
#
# The link between a satellite and a surface point is simulated. The refracted
# path is at about 5 degress above the horizon (at the surface point). 
#
# 2012-08-21, 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)

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


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

# Frequency grid 
#
VectorSet( f_grid, [ 1e9, 5e9, 10e9 ] )


# A pressure grid rougly matching 0 to 80 km, in steps of 250.
#
VectorNLogSpace( p_grid, 321, 1013e2, 1 )


# Definition of species
# 
abs_speciesSet( species=
                ["H2O-PWR98", 
                 "N2-SelfContStandardType",
                 "O2-PWR93"] )


# No line data needed here
# 
abs_lines_per_speciesSetEmpty


# Dimensionality of the atmosphere
#
AtmosphereSet1D

# Atmospheric profiles
# 
AtmRawRead( basename = "testdata/tropical" )
#
AtmFieldsCalc( t_field, z_field, vmr_field, nlte_field, p_grid, 
               lat_grid, lon_grid, t_field_raw, z_field_raw, 
               vmr_field_raw, t_field_raw, atmosphere_dim, 3 )


# Surface altitude
MatrixSetConstant( z_surface, 1, 1, 0 )


# No jacobian calculations
#
jacobianOff

# No scattering
#
cloudboxOff

# Check model atmosphere
# 
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc


# Propagation path agendas and variables
#

# transmitter-receiver path
Copy( ppath_agenda, ppath_agenda__TransmitterReceiverPath )

# refracted path
Copy( ppath_step_agenda, ppath_step_agenda__RefractedPath )

# The reference values are calculated with Thayer's values for k1, k2 and k3:
AgendaSet( refr_index_air_agenda ){
  Ignore( f_grid )
  NumericSet( refr_index_air, 1.0 )
  NumericSet( refr_index_air_group, 1.0 )
  refr_index_airMicrowavesEarth( k1=77.6e-8, k2=64.8e-8, k3=3.776e-3 )
}


#
NumericSet( ppath_lmax, 10e3 )
NumericSet( ppath_lraytrace, 100 )


# Radiative transfer agendas and variables
#
VectorSet( rte_los, [] )   # Dummy value
#
Copy( iy_transmitter_agenda, iy_transmitter_agenda__UnitUnpolIntensity )
Copy( iy_main_agenda, iy_main_agenda__Radiolink )


# Postion of sensor/receiver and transmitter
#
VectorSet( rte_pos, [ 80e3 ] )
VectorSet( rte_pos2, [ 0, 5.1 ] )


# Auxilary variables
#
ArrayOfStringSet( iy_aux_vars, [ 
    "Pressure", "Temperature", "VMR, species 0", 
    "Atmospheric loss", "Absorption, summed", 
    "Free space loss", "Free space attenuation", 
    "Defocusing loss", "Extra path delay", "Bending angle" ] )


abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc

# Run and save
#
iyCalc


# Set "Atmospheric loss" is main variable:
# (If yCalc is used, iyReplaceFromAux shall be included in iy_main_agenda)
#
iyReplaceFromAux( iy, iy_aux, iy_aux_vars, jacobian_do, "Atmospheric loss" )


#WriteXML( "ascii", iy, "iyREFERENCE.xml" )
#WriteXML( "ascii", ppath, "ppath.xml" )


# OK?
# ---
MatrixCreate( iyREFERENCE )
ReadXML( iyREFERENCE, "iyREFERENCE.xml" )
Compare( iy, iyREFERENCE, 1e-5)

}