#DEFINITIONS:  -*-sh-*-
#
# This file tests temperature Jacobian calculations for fully polarised
# simulations. 3D calculations applied, to allow Zeeman.
#
# Two frequencies are done, one with high suface sensitivity and one at 118 GHz
# with a significant Zeeman signature for all Stokes elements. 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_ZeemanPreCalc  )

# 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,118.751e9] )


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


# Line data
#
ReadXML(abs_lines, "line_118ghz.xml")
abs_lines_per_speciesCreateFromLines
Wigner6Init(wigner_initialized, 40000, 100)
zeeman_linerecord_precalcCreateFromLines


# Atmosphere
#
AtmosphereSet3D
VectorNLogSpace( p_grid, 201, 1013e2, 0.05 )
VectorSet( lat_grid, [-10,10] )
Copy( lon_grid, lat_grid)
AtmRawRead( basename = "testdata/tropical" )
#
AtmFieldsCalcExpand1D


# Magnetic field
#
# Craete a synthetic field, where all componenets have the same value
#
Copy( mag_u_field, t_field )
Tensor3Scale( mag_u_field, mag_u_field, 0 )
Tensor3AddScalar( mag_u_field, mag_u_field, 25e-6 )
Copy( mag_v_field, mag_u_field )
Copy( mag_w_field, mag_u_field )


# 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,0,0] )
MatrixSet( sensor_los, [140,45] )


# 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


# Check y against reference
#
VectorCreate( yref )
ReadXML( yref, "yREF4.xml" )
Compare( y, yref, 1e-5,
         "Calculated *y* does not agree with saved reference values." )


# Copy Jacobian
#
MatrixCreate( jcopy )
Copy( jcopy, jacobian )


# Save
#
#output_file_formatSetAscii
#WriteXML( output_file_format, f_grid, "f.xml" )
#WriteXML( output_file_format, z_field, "z.xml" )
#WriteXML( output_file_format, y, "y.xml" )
#WriteXML( output_file_format, jacobian, "Ja.xml" )


# Re-do by perturbations
#
jacobianInit
jacobianAddTemperature( g1=p_grid, g2=lat_grid, g3=lon_grid, hse="on",
                        method="perturbation", dt=0.1 )
jacobianClose
#
yCalc
#
#WriteXML( output_file_format, jacobian, "Jp.xml" )


# To make test more stringent, apply different tolerances for different Stokes
#
VectorCreate( v1 )
VectorCreate( v2 )


# Compare frequency 1, Stokes element 1
#
VectorExtractFromMatrix( v1, jcopy, 0, "row" )
VectorExtractFromMatrix( v2, jacobian, 0, "row" )
Compare( v1, v2, 1e-4,
         "Disagreement for first frequncy and first Stokes element" )

# Compare frequency 1, Stokes element 2
#
VectorExtractFromMatrix( v1, jcopy, 1, "row" )
VectorExtractFromMatrix( v2, jacobian, 1, "row" )
Compare( v1, v2, 1e-5,
         "Disagreement for first frequncy and second Stokes element" )

# Compare frequency 1, Stokes element 3
#
VectorExtractFromMatrix( v1, jcopy, 2, "row" )
VectorExtractFromMatrix( v2, jacobian, 2, "row" )
Compare( v1, v2, 1e-5,
         "Disagreement for first frequncy and third Stokes element" )

# Compare frequency 1, Stokes element 4
#
VectorExtractFromMatrix( v1, jcopy, 3, "row" )
VectorExtractFromMatrix( v2, jacobian, 3, "row" )
Compare( v1, v2, 1e-5,
         "Disagreement for first frequncy and fourth Stokes element" )
}