#DEFINITIONS:  -*-sh-*-
#
# Demonstration of a radio occultation calculation.
#
# The link between two satellites is simulated. Just a few satellite positions
# are included for speed considerations.
#
# 2013-02-07, 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 ] )


# 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, nlte_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 )

#
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, 200 )


# Radiative transfer agendas and variables
#
Copy( iy_transmitter_agenda, iy_transmitter_agenda__UnitUnpolIntensity )
AgendaSet( iy_main_agenda ){
  Ignore( iy_unit )
  Ignore( rte_los )
  Ignore( nlte_field )
  Ignore( diy_dx )
  Ignore( iy_id )
  iyRadioLink( defocus_method = 2 ) 
}



# Sensor/receiver and transmitter
#
IndexCreate( npos )
VectorCreate( transmitter_lat )
#
IndexSet( npos, 6 )
VectorNLinSpace( transmitter_lat, npos, 47, 48.5 )
#
MatrixSetConstant( sensor_pos,      npos, 1, 600e3 )
MatrixSetConstant( sensor_los,      npos, 1, 0     )    # Dummy value
MatrixSetConstant( transmitter_pos, npos, 1, 600e3 )
Append( transmitter_pos, transmitter_lat, "trailing" )
#
sensorOff
sensor_checkedCalc


# Auxilary variables
#
# These are the ones that can be obtained by yCalc.
#
ArrayOfStringSet( iy_aux_vars, [ 
  "Bending angle", 
  "Extra path delay",
  "Free space loss", 
  "Atmospheric loss", 
  "Defocusing loss" ] )

abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc

# Calculate
#
yCalc


# Extract bending angle
#
VectorCreate( ba )
Extract( ba, y_aux, 0 )

#WriteXML( "ascii", ba, "baREFERENCE.xml" )



# OK?
# ---
VectorCreate( baREFERENCE )
ReadXML( baREFERENCE, "baREFERENCE.xml" )
Compare( ba, baREFERENCE, 1e-4)


}