#DEFINITIONS:  -*-sh-*-
# A simple test and demonstration of *yCalcAppend*.
#
# A retrieval based on a combination of three sources is assumed:
# 1. Emission spectrum of 110.8 GHz ozone transition
# 2. A solar occultation spectrum of the same transition
# 3. Data from a 2-channel tropospheric microwave radiometer
#
# The set-up and selection of retrieval quantities are not realistic, rather
# selected to test the code in various ways.
#
# Author: Patrick Eriksson


Arts2 {

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

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

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

# Planck as blackbody radiation
Copy( blackbody_radiation_agenda, blackbody_radiation_agenda__Planck )
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

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



# ---- Basic and common stuff -----------------------------------------------

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

# Frequency grid
VectorNLinSpace( f_grid, 201, 110e9, 111e9 )

# Dimensionality of the atmosphere
AtmosphereSet1D

# Species and absorption
abs_speciesSet( species=[ 
        "O3",
        "H2O-PWR98",
        "N2-SelfContStandardType",
        "O2-PWR98" ] )
ReadXML( abs_lines, "ozone_line.xml" )
abs_lines_per_speciesCreateFromLines

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

# A pressure grid rougly matching 0 to 80 km in 500 m steps.
VectorNLogSpace( p_grid, 160, 1.013e5, 1 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
                                  "testdata/tropical" )
AtmFieldsCalc


# No cloudbox
cloudboxOff

abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc



# ---- Part 1: An emission measurement at 110.8 GHz --------------------------

Copy( iy_main_agenda, iy_main_agenda__Emission )


MatrixSetConstant( sensor_pos, 1, 1, 0 )
MatrixSetConstant( sensor_los, 1, 1, 0 )
sensorOff
sensor_checkedCalc

StringSet( iy_unit, "RJBT" )

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

jacobianInit
jacobianAddAbsSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid,
                       species="O3", method="analytical", unit="rel" )
jacobianAddPolyfit( poly_order=1 )
jacobianClose

yCalc



# ---- Part 2: Corresponding transmission (at another angle) -----------------

Copy( iy_transmitter_agenda, iy_transmitter_agenda__UnitUnpolIntensity )
Copy( iy_main_agenda, iy_main_agenda__Transmission )

StringSet( iy_unit, "1" )

MatrixSetConstant( sensor_los, 1, 1, 45 )
sensorOff
sensor_checkedCalc

ArrayOfStringSet( iy_aux_vars, ["Optical depth"] )

ArrayOfRetrievalQuantityCreate( jacobian_quantities_copy )
Copy( jacobian_quantities_copy, jacobian_quantities )
ArrayOfArrayOfIndexCreate( jacobian_indices_copy )
Copy( jacobian_indices_copy, jacobian_indices )

jacobianInit
jacobianAddAbsSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid,
                       species="H2O-PWR98", method="analytical", unit="rel" )

VectorCreate( rgrid ) # Here just to test check of consistency between grids
#VectorNLogSpace( rgrid, 159, 1.013e5, 1 )
Copy( rgrid, p_grid )
jacobianAddAbsSpecies( g1=rgrid, g2=lat_grid, g3=lon_grid,
                       species="O3", method="analytical", unit="rel" )
jacobianAddPolyfit( poly_order=0 )
jacobianClose

yCalcAppend( jacobian_quantities_copy=jacobian_quantities_copy, 
             jacobian_indices_copy=jacobian_indices_copy,
             append_instrument_wfs=0 )



# ---- Part 3: Data from tropospheric water vapour radiometer -----------------

VectorSet( f_grid, [ 23.4e9, 31e9 ] )

Copy( iy_main_agenda, iy_main_agenda__Emission )

MatrixSetConstant( sensor_los, 1, 1, 0 )
sensorOff
sensor_checkedCalc

StringSet( iy_unit, "RJBT" )

ArrayOfStringSet( iy_aux_vars, [] )

Copy( jacobian_quantities_copy, jacobian_quantities )
Copy( jacobian_indices_copy, jacobian_indices )

jacobianInit
jacobianAddAbsSpecies( g1=p_grid, g2=lat_grid, g3=lon_grid,
                       species="H2O-PWR98", method="analytical", unit="rel" )
jacobianClose

yCalcAppend( jacobian_quantities_copy=jacobian_quantities_copy, 
             jacobian_indices_copy=jacobian_indices_copy )



# Save reults
WriteXML( output_file_format, y,     "TestYCalcAppend.y.xml" )
WriteXML( output_file_format, y_f,   "TestYCalcAppend.y_f.xml" )
WriteXML( output_file_format, y_pol, "TestYCalcAppend.y_pol.xml" )
WriteXML( output_file_format, y_pos, "TestYCalcAppend.y_pos.xml" )
WriteXML( output_file_format, y_los, "TestYCalcAppend.y_los.xml" )
WriteXML( output_file_format, y_aux, "TestYCalcAppend.y_aux.xml" )
#
WriteXML( output_file_format, jacobian, 
                              "TestYCalcAppend.jacobian.xml" )
WriteXML( output_file_format, jacobian_quantities, 
                              "TestYCalcAppend.jacobian_quantities.xml" )
WriteXML( output_file_format, jacobian_indices, 
                              "TestYCalcAppend.jacobian_indices.xml" )
}