#DEFINITIONS:  -*-sh-*-
#
# This file tests consistency between Stokes 1, 2, 3 and 4, involving a
# polaristaion signature, across the Stokes elements, created by the surface
#
# One frequency is done. The surface is modelled by FASTEM.
#
# 2018-11-18, 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, with Zeeman
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# Geometrical path calculation (i.e., refraction neglected)
#
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

# Standard RT agendas
#
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )
Copy( iy_main_agenda, iy_main_agenda__Emission )


# Frequencies and Stokes dim.
#
IndexSet( stokes_dim, 4 )
VectorSet( f_grid, [35e9] )


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


# No line data needed here
# 
abs_lines_per_speciesSetEmpty


# Atmosphere
#
AtmosphereSet1D
VectorNLogSpace( p_grid, 101, 1013e2, 100e2 )
AtmRawRead( basename = "testdata/tropical" )
#
AtmFieldsCalc


# Surface
#
# Don't use interpolation of t_field to set surface temperature. That will
# cause a difference between analytical and perturbation Jacobian
#
Extract( z_surface, z_field, 0 )
Extract( t_surface, t_field, 0 )
VectorCreate( trv )
nelemGet( v=f_grid )
VectorSetConstant( trv, nelem, 0.9 )
AgendaSet( surface_rtprop_agenda ){
  specular_losCalc
  InterpSurfaceFieldToPosition( out=surface_skin_t, field=t_surface )
  surfaceFastem( wind_speed=5, wind_direction=45, fastem_version=5,
                 transmittance=trv )
}


# Sensor pos and los
#
MatrixSet( sensor_pos, [820e3] )
MatrixSet( sensor_los, [140] )


# Define analytical Jacobian
#
jacobianInit
jacobianAddTemperature( g1=p_grid, g2=lat_grid, g3=lon_grid, hse="on" )
jacobianClose


# Deactive parts not used
#
cloudboxOff
sensorOff


# Checks
#
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc( bad_partition_functions_ok = 1 )
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc


# HSE
#
VectorSet( lat_true, [0] )
VectorSet( lon_true, [0] )
#
Extract( p_hse, p_grid, 0 )
NumericSet( z_hse_accuracy, 0.5 )
z_fieldFromHSE


# Run RT calcs
#
StringSet( iy_unit, "RJBT" )
#
yCalc


# Copy results
#
NumericCreate( i4 )
NumericCreate( q4 )
NumericCreate( u4 )
Extract( i4, y, 0 )
Extract( q4, y, 1 )
Extract( u4, y, 2 )
#
VectorCreate( ji4 )
VectorCreate( jq4 )
VectorCreate( ju4 )
VectorExtractFromMatrix( ji4, jacobian, 0, "row" )
VectorExtractFromMatrix( jq4, jacobian, 1, "row" )
VectorExtractFromMatrix( ju4, jacobian, 2, "row" )


# Required agreement
#
NumericCreate( delta )
NumericSet( delta, 1e-9 )


# Stokes = 3
#
IndexSet( stokes_dim, 3 )
sensorOff
yCalc
#
NumericCreate( i )
NumericCreate( q )
NumericCreate( u )
Extract( i, y, 0 )
Extract( q, y, 1 )
Extract( u, y, 2 )
Compare( i, i4, delta, "Disagreement for I beteen Stokes 4 and 3." )
Compare( q, q4, delta, "Disagreement for Q beteen Stokes 4 and 3." )
Compare( u, u4, delta, "Disagreement for U beteen Stokes 4 and 3." )
#
VectorCreate( ji )
VectorCreate( jq )
VectorCreate( ju )
VectorExtractFromMatrix( ji, jacobian, 0, "row" )
VectorExtractFromMatrix( jq, jacobian, 1, "row" )
VectorExtractFromMatrix( ju, jacobian, 2, "row" )
Compare( ji, ji4, delta, "Disagreement for I-part of Jacobian beteen Stokes 4 and 3." )
Compare( jq, jq4, delta, "Disagreement for Q-part of Jacobian beteen Stokes 4 and 3." )
Compare( ju, ju4, delta, "Disagreement for U-part of Jacobian beteen Stokes 4 and 3." )


# Stokes = 2
#
IndexSet( stokes_dim, 2 )
sensorOff
yCalc
#
Extract( i, y, 0 )
Extract( q, y, 1 )
Compare( i, i4, delta, "Disagreement for I beteen Stokes 4 and 2." )
Compare( q, q4, delta, "Disagreement for Q beteen Stokes 4 and 2." )
#
VectorExtractFromMatrix( ji, jacobian, 0, "row" )
VectorExtractFromMatrix( jq, jacobian, 1, "row" )
Compare( ji, ji4, delta, "Disagreement for I-part of Jacobian beteen Stokes 4 and 2." )
Compare( jq, jq4, delta, "Disagreement for Q-part of Jacobian beteen Stokes 4 and 2." )


# Stokes = 1
#
IndexSet( stokes_dim, 1 )
sensorOff
yCalc
#
Extract( i, y, 0 )
Compare( i, i4, delta, "Disagreement for I beteen Stokes 4 and 2." )
#
VectorExtractFromMatrix( ji, jacobian, 0, "row" )
Compare( ji, ji4, delta, "Disagreement for I-part of Jacobian beteen Stokes 4 and 1." )
}