#
# 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 Mars and specifically tests
#
#  - surface altitude:
#      - reading in raw 3D data, preprocess into z_surface variable within a 3D
#         global case (!but no atmfields_checkedCalc - as, on a global
#         scale this fails since altitude reference of altitude data and
#         atmospheric data is inconsistent. with plenty of negative altitudes in
#         surface altitude data and atm data starting only at 5m ARTS rejects
#         the data due to gap between atmo and surface!) (CASE A-1)
#      - reading in raw 3D data, preprocess into z_surface variable within a 3D
#         regional case containing only surface altitudes, and perform
#         atmfields_checkedCalc (CASE A-2, B)
#  - surface refractive index data (CASE B, C):
#      - 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 for regional (CASE B) and global (CASEs C) cases
#  - surface temperature (CASE C): 
#      - reading in raw 3D data, preprocess into t_surface variable, and within
#        surface_rt_prop_agenda derive surface_skin_t from t_surface
#      - loop over all scenarios (4seasons x 2daytimes x 3dustloads)
#      - all in global 3D
#      
# Jana Mendrok 2013-02-26

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

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

# the 3D geo grid for global case
VectorCreate( lat_grid3D_glob )
VectorCreate( lon_grid3D_glob )
VectorLinSpace( lat_grid3D_glob, -90, 90, 5 )
VectorLinSpace( lon_grid3D_glob, 0, 360, 10 )

# the 3D geo grid for regional case (z_surface>5m)
VectorCreate( lat_grid3D_reg )
VectorCreate( lon_grid3D_reg )
VectorLinSpace( lat_grid3D_reg, -49, -18, 1 )
VectorLinSpace( lon_grid3D_reg, 208, 298, 2 )


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

StringCreate( basename )
StringSet( basename, "planets/Mars/MPS/" )

StringCreate( zsurfname )
StringSet( zsurfname, "Mars.z_surface" )
StringCreate( risurfname )
StringSet( risurfname, "Mars.surface_complex_refr_index_field" )
StringCreate( atmcase )
StringCreate( caseext )
StringCreate( casefull )

StringCreate( tsurfname )
StringSet( tsurfname, ".t_surface" )

StringCreate( atmext )
StringSet( atmext, ".sol-avg" )

# Arrays with (sub)case names
ArrayOfStringCreate( seasoncasearray )
ArrayOfStringSet( seasoncasearray, ["Mars.Ls0", "Mars.Ls90", "Mars.Ls180", "Mars.Ls270"] )
ArrayOfStringCreate( timecasearray )
ArrayOfStringSet( timecasearray, [".day", ".night"] )
ArrayOfStringCreate( dustcasearray )
ArrayOfStringSet( dustcasearray, [".dust-high", ".dust-low", ".dust-medium"] )


# a vector for holding reference RT results
VectorCreate(yREFERENCE)

# 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" )

# we manually select a minumim set of basic atm data (main atm constituents)
abs_speciesSet( species=["CO2"] )
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, -19. )
Append( sensor_pos, mtmp, "trailing" )
# sensor longitutde
MatrixSetConstant( mtmp, itmp, 1, 210. )
Append( sensor_pos, mtmp, "trailing" )



#####
# CASE A
# applying surface altitude. fixed surface reflectivity, t_surface from t_field.
# first part global, then reducing to region with z_surface>5m
#####

# some stuff to get a basic atmosphere
#####
Extract( casefull, seasoncasearray, 0 )
Copy( atmcase, casefull )
Extract( casefull, timecasearray, 0 )
Append( atmcase, casefull )
Extract( casefull, dustcasearray, 0 )
Append( atmcase, casefull )
Copy( casefull, basename )
Append( casefull, atmcase )
StringSet( caseext, "/" )
Append( casefull, caseext )
Append( atmcase, atmext )
Append( casefull, atmcase )
Append( casefull, caseext )
Append( casefull, atmcase )
#Print( casefull, 0 )

AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw,
            abs_species, casefull )
p_gridFromZRaw(	p_grid, z_field_raw, 0 )


AtmosphereSet3D

#####
# CASE A-1
#####
Copy( lat_grid, lat_grid3D_glob )
Copy( lon_grid, lon_grid3D_glob )
AtmFieldsCalcExpand1D

# reading the surface altitude field
Copy( casefull, basename )
Append( casefull, zsurfname )
ReadXML( gf2tmp, casefull )
GriddedFieldLatLonRegrid( gf2tmp, lat_grid, lon_grid, gf2tmp )
FieldFromGriddedField( z_surface, p_grid, lat_grid, lon_grid, gf2tmp )

#####
# CASE A-2
#####
Copy( lat_grid, lat_grid3D_reg )
Copy( lon_grid, lon_grid3D_reg )
AtmFieldsCalcExpand1D

# reading the surface altitude field
Copy( casefull, basename )
Append( casefull, zsurfname )
ReadXML( gf2tmp, casefull )
GriddedFieldLatLonRegrid( gf2tmp, lat_grid, lon_grid, gf2tmp )
FieldFromGriddedField( z_surface, p_grid, lat_grid, lon_grid, gf2tmp )

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

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 )   
   #Print( surface_skin_t, 0 )
   surfaceFlatScalarReflectivity
}
yCalc
#Print( y, 0 )



#####
# CASE B
# applying surface altitude and refractive index. t_surface from t_field.
# still regional only.
#####

# reading surface refractive index data (GriddedField5). 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 )

AgendaSet( surface_rtprop_agenda ){
   specular_losCalc
   InterpAtmFieldToPosition( out=surface_skin_t, field=t_field )
   #Print( surface_skin_t, 0 )
   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, 0 )
#WriteXML( in=y )
ReadXML( out=yREFERENCE, filename="TestSurf_Mars.y.xml" )
Compare( y, yREFERENCE, 1e-6 )


#####
# CASE C
# applying refractive index and t_surface. z_surface from lowest atmo-z.
# returning to global case
#####

AtmosphereSet3D
Copy( lat_grid, lat_grid3D_glob )
Copy( lon_grid, lon_grid3D_glob )

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


# we go with several nested foorloop through the different cases.

AgendaCreate( forloop_agenda_dust )
AgendaSet( forloop_agenda_dust ){
  # construct atmcase name III (Mars.LsXX.YY.dust-ZZ)
  Extract( casefull, dustcasearray, forloop_index )
  Append( atmcase, casefull )
  # keep the casestring till dust and make upper-level folder name
  Append( basename, atmcase )
  StringSet( caseext, "/" )
  Append( basename, caseext )
  Copy( casefull, basename )
  Append( casefull, atmcase )

  # reading the surface temperature field
  Append( casefull, tsurfname )
  ReadXML( gf2tmp, casefull )
  GriddedFieldLatLonRegrid( gf2tmp, lat_grid, lon_grid, gf2tmp )

  # reading atmospheric field data
  Append( atmcase, atmext )
  Append( basename, atmcase )
  Append( basename, caseext )
  Append( basename, atmcase )
  AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw,
              abs_species, basename )

  p_gridFromZRaw( p_grid, z_field_raw, 0 )
  AtmFieldsCalcExpand1D
  Extract( z_surface, z_field, 0 )
  FieldFromGriddedField( t_surface, p_grid, lat_grid, lon_grid, gf2tmp )

  atmfields_checkedCalc
  atmgeom_checkedCalc
  cloudbox_checkedCalc
  sensor_checkedCalc

  #Print( casefull, 0 )
  yCalc
  #Print( y, 0 )
}

AgendaCreate( forloop_agenda_time )
AgendaSet( forloop_agenda_time ){
  # construct atmcase name II (Mars.LsXX.d/n)
  Extract( casefull, timecasearray, forloop_index )
  Append( atmcase, casefull )
  Copy( forloop_agenda, forloop_agenda_dust )
  nelemGet( ncases, dustcasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1 )
}

AgendaCreate( forloop_agenda_season )
AgendaSet( forloop_agenda_season ){
  # construct atmcase name I (Mars.LsXX)
  Extract( casefull, seasoncasearray, forloop_index )
  Copy( atmcase, casefull )
  Copy( forloop_agenda, forloop_agenda_time )
  nelemGet( ncases, timecasearray )
  IndexStepDown( ncases, ncases )
  ForLoop( forloop_agenda, 0, ncases, 1 )
}

nelemGet( ncases, seasoncasearray )
IndexStepDown( ncases, ncases )
Copy( forloop_agenda, forloop_agenda_season )
ForLoop( forloop_agenda, 0, ncases, 1  )

}