#
# Demonstration and test of adjusting z_field according to hydrostatic equilibrium 
#
# Jana Mendrok 2013-01-29

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 )


# Dimensionality of the atmosphere
#
AtmosphereSet1D


#
# Some settings resulting in no gas absorption
#

# on-the-fly absorption
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )


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

# No line data needed here
# 
abs_lines_per_speciesSetEmpty




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


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


# Neutral atmosphere, up to about 95 km
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "testdata/tropical" )

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

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


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


# No scattering
#
cloudboxOff

# No jacobian calculations
#
jacobianOff

# Check model atmosphere
# 
atmfields_checkedCalc
cloudbox_checkedCalc



#WriteXML( "ascii", z_field, "z_fieldFromHSE_REFERENCE.xml" )


# Expected results
#
Tensor3Create( zREFERENCE )
#
ReadXML( zREFERENCE, "z_fieldFromHSE_REFERENCE.xml" )

# Check
#
Compare( z_field, zREFERENCE, 1e-4 )


# SIDE ISSUE: testing effect of H2O for z_fieldFromHSE
# (result: dev<50m for z<1400km - basically indicates that H2O is not that important...)
#
# we need to re-prepare the atmo setup without H2O
#abs_speciesSet( species=
#                ["N2",
#                 "O2"] )
#AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
#            "testdata/tropical" )
#ReadXML( t_field_raw, "testdata/tropical.expanded.t.xml" )
#ReadXML( z_field_raw, "testdata/tropical.expanded.z.xml" )
#AtmFieldsCalc( t_field, z_field, vmr_field, p_grid, 
#               lat_grid, lon_grid, t_field_raw, z_field_raw, 
#               vmr_field_raw, atmosphere_dim )
#atmfields_checkedCalc
#atmgeom_checkedCalc
#z_fieldFromHSE
#WriteXML( "ascii", z_field, "z_fieldFromHSE_noH2O.xml" )
}