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


Arts2 {


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


# Select frequency band here:
#

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

# 1D atmosphere
#
AtmosphereSet1D

jacobianOff

INCLUDE "instruments/odinsmr/odinsmr_501.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 )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

# absorption from LUT
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable )
#Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )


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

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

# Atmospheric profiles here taken from Fascod
AtmRawRead( basename = "testdata/tropical" )
#
AtmFieldsCalc

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



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

abs_lines_per_speciesCreateFromLines

abs_lineshapeDefine( shape="Faddeeva_Algorithm_916",
                     forefactor="VVH",
                     cutoff=750e9 )

AbsInputFromAtmFields
abs_speciesSet( abs_species=abs_nls, species=[] )
VectorSet( abs_nls_pert, [] )
VectorSet( abs_t_pert, [] )

abs_xsec_agenda_checkedCalc
abs_lookupCalc

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

# 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, refellipsoid, atmosphere_dim, z_tan )
Matrix1ColFromVector( sensor_los, za )

sensor_checkedCalc

propmat_clearsky_agenda_checkedCalc



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

MatrixCreate( Ja )
MatrixCreate( Jp )
StringCreate( info )

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

StringSet( info, "O3 rel analytical" )
Print( info, 0 )
jacobianInit
jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "O3", "analytical", "rel", 1, 0.01 )
jacobianClose
#WriteXML( in=jacobian_quantities, filename="Jq_O3_analytical.xml" )


# No scattering
#
cloudboxOff
#
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc


yCalc
WriteXML( in=y, filename="y_O3_analytical.xml" )
WriteXML( in=jacobian, filename="J_O3_analytical.xml" )
Copy( Ja, jacobian )


# Same with perturbations
StringSet( info, "O3 rel perturbation" )
Print( info, 0 )
jacobianInit
jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "O3", "perturbation", "rel", 1, 0.01 )
jacobianClose
#WriteXML( in=jacobian_quantities, filename="Jq_O3_perturbation.xml" )
yCalc
#WriteXML( in=y, filename="y_O3_perturbation.xml" )
#WriteXML( in=jacobian, filename="J_O3_perturbation.xml" )
Copy( Jp, jacobian )

# Compare
Compare( Ja, Jp, 0.005 )


# Analytical but now using VMR, though for all "O2" species
StringSet( info, "O2 vmr analytical" )
Print( info, 0 )
jacobianInit
jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
    atmosphere_dim, p_grid, lat_grid, lon_grid, 
    retrieval_grid, lat_grid, lon_grid, 
    "O2", "analytical", "vmr", 0, 0.01 )
jacobianClose
#WriteXML( in=jacobian_quantities, filename="Jq_O2_analytical.xml" )
yCalc
#WriteXML( in=y, filename="y_O2_analytical.xml" )
WriteXML( in=jacobian, filename="J_O2_analytical.xml" )
#Copy( Ja, jacobian )



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

# For limb sounding, the analytical expressions do not cover all effects
# related to HSE and refraction and HSE must be "off" here to get consistent
# results. 

# Stuff needed for HSE
NumericSet( p_hse, 1000e2 )
NumericSet( z_hse_accuracy, 1 )
VectorSet( lat_true, [15] )
VectorSet( lon_true, [123] )

StringSet( info, "T analytical" )
Print( info, 0 )
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
#WriteXML( in=jacobian_quantities, filename="Jq_T_analytical.xml" )
yCalc
#WriteXML( in=y, filename="y_T_analytical.xml" )
WriteXML( in=jacobian, filename="J_T_analytical.xml" )
Copy( Ja, jacobian )


# Same with perturbations
StringSet( info, "T perturbation" )
Print( info, 0 )
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
#WriteXML( in=jacobian_quantities, filename="Jq_T_perturbation.xml" )
yCalc
#WriteXML( in=y, filename="y_T_perturbation.xml" )
#WriteXML( in=jacobian, filename="J_T_perturbation.xml" )
Copy( Jp, jacobian )

# Compare
Compare( Ja, Jp, 0.001 )



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

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

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

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

#WriteXML( "ascii", Ja, "J_pointing_recalc.xml" )
#WriteXML( "ascii", Jp, "J_pointing_interp.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)
Compare( Ja, Jp, 10 )





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

StringSet( info, "Others: Winds" )
Print( info, 0 )
jacobianInit

jacobianAddWind( jacobian_quantities, jacobian_agenda, atmosphere_dim, 
    p_grid, lat_grid, lon_grid, retrieval_grid, lat_grid, lon_grid, "v" )

jacobianAddWind( jacobian_quantities, jacobian_agenda, atmosphere_dim, 
    p_grid, lat_grid, lon_grid, retrieval_grid, lat_grid, lon_grid, "w", 0.1 )
    
jacobianClose
#WriteXML( in=jacobian_quantities, filename="Jq_Winds.xml" )

IndexSet( abs_f_interp_order, 1 )
yCalc
#WriteXML( in=y, filename="y_Winds.xml" )
WriteXML( in=jacobian, filename="J_Winds.xml" )



StringSet( info, "Others: Freqs, Baseline" )
Print( info, 0 )
jacobianInit

jacobianAddFreqShift( df = 50e3 )
jacobianAddFreqStretch( df = 50e3 )

jacobianAddPolyfit( jacobian_quantities, jacobian_agenda, 
    sensor_response_pol_grid, sensor_response_dlos_grid, sensor_pos, 
    1, 0, 0, 0 )

jacobianAddSinefit( jacobian_quantities, jacobian_agenda, 
    sensor_response_pol_grid, sensor_response_dlos_grid, sensor_pos, 
    [200e6,40e6], 0, 0, 0 )

jacobianClose
#WriteXML( in=jacobian_quantities, filename="Jq_Other.xml" )

IndexSet( abs_f_interp_order, 1 )
yCalc
#WriteXML( in=y, filename="y_Other.xml" )
WriteXML( in=jacobian, filename="J_Other.xml" )


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

StringSet( info, "Transmission analytical" )
Print( info, 0 )
# (standard) transmission calculation
Copy( iy_main_agenda, iy_main_agenda__Transmission )
Copy( iy_transmitter_agenda , iy_transmitter_agenda__UnitUnpolIntensity )

# Calculate with analytical 
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
StringSet( info, "Transmission perturbation" )
Print( info, 0 )
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
Compare( Ja, Jp, 0.001 )
}