#DEFINITIONS:  -*-sh-*-
# This is a test of weighting function calculations.
#
# The test is baded on the Odin-SMR 501 GHz case found in another folder.
#
# Author: Patrick Eriksson


Arts2 {


##############################################################################
#
# Initial part
#
##############################################################################


# Select frequency band here:
#
INCLUDE "odin/odinsmr_501.arts"



# ---- Atmospheric scenario --------------------------------------------------

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

# 1D atmosphere
#
AtmosphereSet1D

# Atmospheric profiles here taken from Fascod
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, "../atmosphere_data/tropical" )
#
AtmFieldsCalc

# Get ground altitude (z_surface) from z_field
Extract( z_surface, z_field, 0 )

basics_checkedCalc
cloudbox_checkedCalc



# ---- Create absorption table -----------------------------------------------

abs_lines_per_speciesCreateFromLines

AbsInputFromAtmFields
SpeciesSet( abs_nls, [] )
VectorSet( abs_nls_pert, [] )
VectorSet( abs_t_pert, [] )

abs_lookupCreate

AgendaSet( abs_scalar_gas_agenda ){
  abs_scalar_gasExtractFromLookup
}



# ---- Sensor position and LOS -----------------------------------------------

# Number of tangent altitudes
IndexCreate( n_tan )
IndexSet( n_tan, 3 )

# Sensor position, with platform altitude set to 600 km
MatrixSetConstant( sensor_pos, n_tan, 1, 600e3 )
sensor_posAddRgeoid

# LOS, specified by the corresponding geometrical tangent altitudes
# Tangent altitudes will be equally spaced between 50 and 20 km
VectorCreate( z_tan )
VectorNLinSpace( z_tan, n_tan, 50e3, 20e3 )
VectorCreate( za )
VectorZtanToZa1D( za, sensor_pos, r_geoid, atmosphere_dim, z_tan )
Matrix1ColFromVector( sensor_los, za )





##############################################################################
#
# Absorption
#
##############################################################################

MatrixCreate( Ja )
MatrixCreate( Jp )

# Species (just O3)
#
# Retrieve for a grid rougly matching 16 to 64 km in 2 km steps.
#
VectorCreate( retrieval_grid )
VectorNLogSpace( retrieval_grid, 25, 100e2, 10 )

jacobianInit
jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "O3", "analytical", "rel", 0.01 )
jacobianClose
yCalc
Copy( Ja, jacobian )


# Same with perturbations
jacobianInit
jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "O3", "perturbation", "rel", 0.01 )
jacobianClose
yCalc
Copy( Jp, jacobian )

# Compare
MatrixCompare( Ja, Jp, 0.005 )





##############################################################################
#
# Temperature, without HSE
#
##############################################################################

# Dummy values
NumericSet( p_hse, 1000e2 )
NumericSet( z_hse_accuracy, 1 )

jacobianInit
jacobianAddTemperature( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "off", "analytical", 0.1 )
jacobianClose
yCalc
Copy( Ja, jacobian )


# Same with perturbations
jacobianInit
jacobianAddTemperature( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "off", "perturbation", 0.1 )
jacobianClose
yCalc
Copy( Jp, jacobian )

# Compare
MatrixCompare( Ja, Jp, 0.001 )





##############################################################################
#
# Pointing
#
##############################################################################

# Sensor time must be specified here
nrowsGet( nrows, sensor_pos )
VectorNLinSpace( sensor_time, nrows, 0, 1 )

jacobianInit
jacobianAddPointingZa( jacobian_quantities, jacobian_agenda, 
                       sensor_pos, sensor_time, 0, "recalc", 0.001 )
jacobianClose
yCalc
Copy( Ja, jacobian )

jacobianInit
jacobianAddPointingZa( jacobian_quantities, jacobian_agenda, 
                       sensor_pos, sensor_time, 0, "interp", 0.001 )
jacobianClose
yCalc
Copy( Jp, jacobian )

#WriteXML( "ascii", Ja, "Ja.xml" )
#WriteXML( "ascii", Jp, "Jp.xml" )

# Compare (Note that the WF is for 1 deg change, corresponding to about
# 60 km change in tangent altitude, and 10 K/deg accuracy is OK)
MatrixCompare( Ja, Jp, 10 )





##############################################################################
#
# Just check that remaining weighting functions don't cause any error
#
##############################################################################

jacobianInit

# Another test of jacobianAddFreqShiftAndStretch is found in Groundbased
jacobianAddFreqShiftAndStretch( jacobian_quantities, jacobian_agenda, 
                                f_grid, "interp", 50e3, 1 )

jacobianAddPolyfit( jacobian_quantities, jacobian_agenda, 
    sensor_response_pol_grid, sensor_response_f_grid, 
    sensor_response_za_grid, sensor_pos, 1, 0, 0, 0 )

jacobianClose


}