#DEFINITIONS:  -*-sh-*-
#
# Demonstration and test of calculation including Faraday rotation 
#
# The rotation for a signal transmitted along the zenith direction is
# calculated.
#
# 2013-03-20, Patrick Eriksson

Arts2 {

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


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


# Frequency grid 
#
VectorNLogSpace( f_grid, 101, 0.1e9, 5e9 )
#VectorSet( f_grid, [0.1e9] )


# Dimensionality of the atmosphere
#
AtmosphereSet1D


# A pressure grid rougly matching 0 to 1000 km, in steps of 2km.
#
VectorNLogSpace( p_grid, 501, 1013e2, 1e-80 )


# Tempature and z covering the ionosphere
#
ReadXML( t_field_raw, "testdata/tropical.expanded.t.xml" )
ReadXML( z_field_raw, "testdata/tropical.expanded.z.xml" )


#
# Include som gas species
#
# (For demonstration only, no absorption set)

abs_speciesInit

GriddedField3Create( raw3 )

# N2
#
abs_speciesAdd( species = [ "N2" ] )
ReadXML( raw3, "testdata/tropical.N2.xml"  )
Append( vmr_field_raw, raw3 )

# O2
#
abs_speciesAdd( species = [ "O2" ] )
ReadXML( raw3, "testdata/tropical.O2.xml"  )
Append( vmr_field_raw, raw3 )

# H2O
#
abs_speciesAdd( species = [ "H2O" ] )
ReadXML( raw3, "testdata/tropical.H2O.xml"  )
Append( vmr_field_raw, raw3 )


#
# Add free electrons
#

abs_speciesAdd( species = [ "free_electrons" ] )
ReadXML( raw3, "testdata/ne_iri_solmax_spring_12UTC_0latlon.xml" )
Append( vmr_field_raw, raw3 )


# Interpolate to p_grid (VMR is "zero padded")
#
AtmFieldsCalc( vmr_zeropadding=1 )


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


# Apply HSE
#
# (this roughly recreates the original altitudes for the input electron density
#  and magnetic field)
#
VectorSet( lat_true, [0] )
VectorSet( lon_true, [0] )
NumericSet( p_hse, 1013e2 )
NumericSet( z_hse_accuracy, 10 )
#
atmfields_checkedCalc
atmgeom_checkedCalc
z_fieldFromHSE


# Magnetic field components
ReadXML( raw3, "testdata/bu_igrf11_2000_0latlon.xml" )
GriddedFieldPRegrid( raw3, p_grid, raw3 )
FieldFromGriddedField( mag_u_field, p_grid, lat_grid, lon_grid, raw3 )
ReadXML( raw3, "testdata/bv_igrf11_2000_0latlon.xml" )
GriddedFieldPRegrid( raw3, p_grid, raw3 )
FieldFromGriddedField( mag_v_field, p_grid, lat_grid, lon_grid, raw3 )
ReadXML( raw3, "testdata/bw_igrf11_2000_0latlon.xml" )
GriddedFieldPRegrid( raw3, p_grid, raw3 )
FieldFromGriddedField( mag_w_field, p_grid, lat_grid, lon_grid, raw3 )



#
# Absorption including Faraday rotation
#

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

# No line data needed here
# 
abs_lines_per_speciesSetEmpty

Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly_Faraday )




# No scattering
#
cloudboxOff

# No jacobian calculations
#
jacobianOff

# Check model atmosphere
# 
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc


# Propagation path agendas and variables
#
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )
NumericSet( ppath_lmax, 10e3 )



# Radiative transfer agendas and variables
#
AgendaSet( iy_transmitter_agenda ){
  Ignore( rtp_pos )
  Ignore( rtp_los )
  iy_transmitterSinglePol
}
Copy( iy_main_agenda, iy_main_agenda__Transmission )


# Sensor/receiver and transmitter
#
MatrixSet( sensor_pos, [ 0 ] )
MatrixSet( sensor_los, [ 0 ] )    
MatrixSet( transmitter_pos, [] ) # Dummy value
ArrayOfIndexSet( sensor_pol, [ 5 ] )
#
sensorOff
sensor_checkedCalc


# Auxilary variables
#
ArrayOfStringSet( iy_aux_vars, [ 
  "Faraday rotation" ] )


# Calculate 
#
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
yCalc


# Extraxt total Faraday rotation
VectorCreate( farrot_total )
Extract( farrot_total, y_aux, 0 )


WriteXML( "ascii", f_grid,      "f.xml"        )
WriteXML( "ascii", vmr_field,   "vmr.xml"      )
WriteXML( "ascii", z_field,     "z_field.xml"  )
WriteXML( "ascii", mag_w_field, "bw_field.xml" )
WriteXML( "ascii", farrot_total, "farrot.xml"  )

#WriteXML( "ascii", y, "yREFERENCE.xml" )
#WriteXML( "ascii", farrot_total, "farrot_totalREFERENCE.xml" )


# Expected results
#
VectorCreate( yREFERENCE )
VectorCreate( frtotREFERENCE )
#
ReadXML( yREFERENCE, "yREFERENCE.xml" )
ReadXML( frtotREFERENCE, "farrot_totalREFERENCE.xml" )


# Check
#
Compare( y, yREFERENCE, 1e-4 )
Compare( farrot_total, frtotREFERENCE, 1e-2 )

}