#
# Testing functionality (meeting format requirements, etc.) of surface related
#  data.
#
# General test setup: reading in raw data (including a basic atmosphere),
#  extracting/regridding, executing standard pre-RT calc internal test method
#  atmfields_checkedCalc, and performing some RT simulations in order to apply data
#  that has no dedicated check method (e.g., when only used through agendas like
#  surface reflectivity).
#
#
# This case is for Venus and specifically tests
#
#  - surface temperature (CASE A): reading in raw 3D data, preprocess into
#     t_surface variable, and within surface_rt_prop_agenda derive surface_skin_t
#     from t_surface (CASE A-0 first does without t_surface, taking
#     surface_skin_t from atmospheric t_field). all in 3D.
#  - surface refractive index data (CASE B):
#      - reading 3D raw data
#      - in surface_rt_prop_agenda deriving surface emission/reflection field
#         from complex refractive index data (using surfaceFlatRefractiveIndex)
#      - all in 3D only
#      
# Jana Mendrok 2013-02-26

Arts2 {
INCLUDE "general/general.arts"
INCLUDE "general/continua.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/planet_venus.arts"

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

# the 3D geo grid to test
VectorCreate( lat_grid3D )
VectorCreate( lon_grid3D )
VectorLinSpace( lat_grid3D, -90, 90, 18 )
VectorLinSpace( lon_grid3D, -20, 350, 18 )


GriddedField2Create( gf2tmp )
MatrixCreate( mtmp )
IndexCreate( itmp )


StringCreate( basename )
StringSet( basename, "planets/Venus/MPS/" )
StringCreate( atmcase )
StringSet( atmcase, "Venus.spicav.night/Venus.spicav.night" )
StringCreate( tsurfname )
StringSet( tsurfname, "Venus.t_surface" )
StringCreate( risurfname )
StringSet( risurfname, "Venus.surface_complex_refr_index_field" )
StringCreate( casefull )


# some stuff to get a basic atmosphere
#####
# we manually select a minumim set of basic atm data (main atm constituents)
abs_speciesSet( species=["CO2"] )
Copy( casefull, basename )
Append( casefull, atmcase )
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw,
            abs_species, casefull )
p_gridFromZRaw(	p_grid, z_field_raw, 0 )


# and some further settings in order to be able to do an RT calc
#####
cloudboxOff
jacobianOff
IndexSet( stokes_dim, 1 )
VectorSet( f_grid, [300e9] )
sensorOff
StringSet( iy_unit, "PlanckBT" )
abs_linesReadFromSplitArtscat(abs_lines, abs_species, "spectroscopy/Perrin/", 0, 1e12)
abs_lines_per_speciesCreateFromLines

# and agenda settings needed for RT calc
#####
Copy( iy_main_agenda, iy_main_agenda__Emission )
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )
Copy( blackbody_radiation_agenda, blackbody_radiation_agenda__Planck )
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )




# sensor placed over Maxwell Montes region scanning from the high to low surface
#  RI region
# LOS zenith angle
MatrixSet( sensor_los, [180;130;115;113.8] )
# LOS azimuth angle
#MatrixSet( mtmp,       [])
nrowsGet( itmp, sensor_los )
MatrixSetConstant( mtmp, itmp, 1, 120. )
Append( sensor_los, mtmp, "trailing" )

# sensor altitude
MatrixSetConstant( sensor_pos, itmp, 1, 600e3 )
# sensor latitude
MatrixSetConstant( mtmp, itmp, 1, 60. )
Append( sensor_pos, mtmp, "trailing" )
# sensor longitutde
MatrixSetConstant( mtmp, itmp, 1, 60. )
Append( sensor_pos, mtmp, "trailing" )


#####
# CASEs A: first only using t_surface data. reflectivity fixed.
#####

AtmosphereSet3D
Copy( lat_grid, lat_grid3D )
Copy( lon_grid, lon_grid3D )
AtmFieldsCalcExpand1D( vmr_zeropadding=1 )
Extract( z_surface, z_field, 1 )

# reading and regridding the surface temperature field
Copy( casefull, basename )
Append( casefull, tsurfname )
ReadXML( gf2tmp, casefull )
GriddedFieldLatLonRegrid( gf2tmp, lat_grid, lon_grid, gf2tmp )
FieldFromGriddedField( t_surface, p_grid, lat_grid, lon_grid, gf2tmp )

# reading surface refractive index data (GriddedField4). no regridding here. we
#  need to do that in form of 1D only inside the surface_rtprop_agenda, as we
#  only know then which exact point(s) we need. for better traceability (and
#  since here this isn't just a temporary field), we create a dedicated
#  workspace variable for this data.
Copy( casefull, basename )
Append( casefull, risurfname )
GriddedField5Create( ri_surface )
ReadXML( ri_surface, casefull )

abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc
atmgeom_checkedCalc

# now we need to do some RT calc in order to APPLY the reflectivity data
cloudbox_checkedCalc
sensor_checkedCalc


# CASE A-0 (we actually don't use t_surface, but surface_skin_t from t_field
#####

VectorSet( surface_scalar_reflectivity, [0.4] )

# surface temp from atmospheric t_field
AgendaSet( surface_rtprop_agenda ){
   specular_losCalc
   InterpAtmFieldToPosition( out=surface_skin_t, field=t_field )
   surfaceFlatScalarReflectivity
}
#yCalc
#Print( y )


# CASE A
#####

# surface temp from t_surface field
AgendaSet( surface_rtprop_agenda ){
  specular_losCalc
  InterpSurfaceFieldToPosition( out=surface_skin_t, field=t_surface )  
  surfaceFlatScalarReflectivity
} 
yCalc
#Print( y )


#####
# CASE B: now also applying surface reflection properties from surface complex
#  refractive index data
#####

AgendaSet( surface_rtprop_agenda ){
   specular_losCalc
   InterpSurfaceFieldToPosition( out=surface_skin_t, field=t_surface )  
   Select( lat_true, rtp_pos, [1] )
   Select( lon_true, rtp_pos, [2] )
   surface_complex_refr_indexFromGriddedField5( complex_refr_index_field=ri_surface )
   surfaceFlatRefractiveIndex
} 
yCalc
#Print( y )


}