#DEFINITIONS:  -*-sh-*-
#
# ARTS control file for testing 1D propagation path calculations
# with different refractive index calculation methods (None, Thayer, MWgeneral)
# for the different planets (particularly dependence on atmospheric composition)
#
# Jana Mendrok 2013-02-25

Arts2{

INCLUDE "general/general.arts"
INCLUDE "general/agendas.arts"

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

# Use sensor_pos/los to store each case
VectorCreate( za )
VectorCreate( ztan )
VectorSet( ztan, [50e3, 46e3, 34e3, 10e3, 5e3, 2.5e3] )
nelemGet( nelem, ztan )
MatrixSetConstant( sensor_pos, nelem, 1, 600e3 )
VectorSet( rte_pos2, [] ) 

IndexCreate( ilast )
nrowsGet( ilast, sensor_pos )
IndexStepDown( ilast, ilast )


Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )
IndexSet( stokes_dim, 1 )
AtmosphereSet1D
cloudboxOff
NumericSet( ppath_lmax, 5e3 )

# A dummy frequency grid
VectorSet( f_grid, [10e9] )


IndexCreate( tpoutlev )
IndexSet( tpoutlev, 1 )
IndexCreate( outlev )
IndexSet( outlev, 0 )
NumericCreate( tanh1 )
NumericCreate( tanlat1 )
VectorCreate( tp1 )
NumericCreate( tanh2 )
NumericCreate( tanlat2 )
VectorCreate( tp2 )
NumericCreate( diff )
StringCreate( infostring )
AgendaSet( forloop_agenda ){
  VectorExtractFromMatrix( rte_pos, sensor_pos, forloop_index, "row" )
  VectorExtractFromMatrix( rte_los, sensor_los, forloop_index, "row" )

  ###
  # CASE 0-0: geometric path (no refraction)
  Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )
  ppathCalc
  TangentPointExtract( tan_pos=tp1 )
  TangentPointPrint( level=tpoutlev )

  ###
  # CASE 0-1: refracted path with n==1 (no refraction)
  Copy( ppath_step_agenda, ppath_step_agenda__RefractedPath )
  Copy( refr_index_air_agenda, refr_index_air_agenda__NoRefrac )
  ppathCalc
  TangentPointExtract( tan_pos=tp2 )
  TangentPointPrint( level=tpoutlev )

  # Compare geom and n==1 results
  Extract( tanh1, tp1, 0 )
  Extract( tanh2, tp2, 0 )
  Compare( tanh1, tanh2, 1.,
           "Geometric and n==1 tangent point altitude differ too much." )
  Extract( tanlat1, tp1, 1 )
  Extract( tanlat2, tp2, 1 )
  Compare( tanlat1, tanlat2, 0.04,
           "Geometric and n==1 tangent point latitude differ too much." )

  ###
  # CASE 1: refracted path with n from Thayer
  Copy( refr_index_air_agenda, refr_index_air_agenda__GasThayer )
  ppathCalc
  #WriteXMLIndexed( "ascii", forloop_index, ppath, "thayer" )
  TangentPointExtract( tan_pos=tp1 )
  TangentPointPrint( level=tpoutlev )

  ###
  # CASE 2: refracted path with n from MWgeneral
  Copy( refr_index_air_agenda, refr_index_air_agenda__GasMWgeneral )
  ppathCalc
  #WriteXMLIndexed( "ascii", forloop_index, ppath, "mwgeneral" )
  TangentPointExtract( tan_pos=tp2 )
  TangentPointPrint( level=tpoutlev )
  
  # Print diff of Thayer and MWgeneral
  Extract( tanh1, tp1, 0 )
  Extract( tanh2, tp2, 0 )
  NumericScale( tanh2, tanh2, -1 )
  NumericAdd( diff, tanh1, tanh2 )
  StringSet( infostring, "Tangent altitude diff of Thayer and MWgeneral [m]: " )
  Print( infostring, outlev )
  Print( diff, outlev )

  #Extract( tanlat1, tp1, 1 )
  #Extract( tanlat2, tp2, 1 )
  #NumericScale( tanlat2, tanlat2, -1 )
  #NumericAdd( diff, tanlat1, tanlat2 )
  #StringSet( infostring, "Tangent latitude diff of Thayer and MWgeneral [deg]: " )
  #Print( infostring, outlev )
  #Print( diff, outlev )
}


#####
# CASE A: Earth
#####

# Main atmospheric contributors (water vapour considered by Thayer refraction.
#  However, when comparing effects of planetary composition, we don't want
#  further confusion by H2O. so we leave it out here.)
#abs_speciesSet( species=["H2O","N2","O2"] )
abs_speciesSet( species=["N2","O2","CO2"] )

refellipsoidEarth( refellipsoid, "Sphere" )
# convert (geometric) tangent altitudes to zenith angles at observer
VectorZtanToZa1D( za, sensor_pos, refellipsoid, atmosphere_dim, ztan )
Matrix1ColFromVector( sensor_los, za )

# Pressure grid rougly matching 0 to 80 km.
VectorNLogSpace( p_grid, 41, 1000e2, 1 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "planets/Earth/Fascod/tropical/tropical" )
AtmFieldsCalc
Extract( z_surface, z_field, 0 )

atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc

ForLoop( forloop_agenda, 0, ilast, 1  )


#####
# CASE B: Mars
#####

# Main atmospheric contributors
abs_speciesSet( species=["N2","O2","CO2","H2"] )

refellipsoidMars( refellipsoid, "Sphere" )
# convert (geometric) tangent altitudes to zenith angles at observer
VectorZtanToZa1D( za, sensor_pos, refellipsoid, atmosphere_dim, ztan )
Matrix1ColFromVector( sensor_los, za )

# Pressure grid rougly matching 0 to 80 km.
VectorNLogSpace( p_grid, 41, 765, 0.08 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "planets/Mars/MPS/Mars.Ls0.day.dust-medium/Mars.Ls0.day.dust-medium.sol-avg/Mars.Ls0.day.dust-medium.sol-avg" )
AtmFieldsCalc
Extract( z_surface, z_field, 0 )

atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc

ForLoop( forloop_agenda, 0, ilast, 1  )


#####
# CASE C: Venus
#####

# Main atmospheric contributors
abs_speciesSet( species=["N2","O2","CO2"] )

refellipsoidVenus( refellipsoid, "Sphere" )
# convert (geometric) tangent altitudes to zenith angles at observer
VectorZtanToZa1D( za, sensor_pos, refellipsoid, atmosphere_dim, ztan )
Matrix1ColFromVector( sensor_los, za )

# Pressure grid rougly matching 0 to 80 km.
VectorNLogSpace( p_grid, 41, 9.2e6, 4.9e2 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "planets/Venus/MPS/Venus.vira.day/Venus.vira.day" )
AtmFieldsCalc( vmr_zeropadding=1 )
Extract( z_surface, z_field, 0 )

atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc

ForLoop( forloop_agenda, 0, ilast, 1  )


#####
# CASE D: Jupiter
#####

# Main atmospheric contributors
abs_speciesSet( species=["CO2","H2","He"] )

refellipsoidJupiter( refellipsoid, "Sphere" )
# convert (geometric) tangent altitudes to zenith angles at observer
VectorZtanToZa1D( za, sensor_pos, refellipsoid, atmosphere_dim, ztan )
Matrix1ColFromVector( sensor_los, za )

# Pressure grid rougly matching 0 to 80 km.
VectorNLogSpace( p_grid, 41, 1.0e5, 1.8e3 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "planets/Jupiter/MPS/Jupiter.mean/Jupiter.mean" )
AtmFieldsCalc( vmr_zeropadding=1 )
Extract( z_surface, z_field, 0 )

atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc

ForLoop( forloop_agenda, 0, ilast, 1  )

}