#DEFINITIONS:  -*-sh-*-

# --------------------------------------------------------------------
# !!! WARNING !!! Untested !!! USE WITH CARE !!!
# Test batch calculations for SEVIRI.
# Based on HIRS test by Viju Oommen John and Stefan Buehler, August to
# October 2008.
# --------------------------------------------------------------------


Arts2 {

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

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

# (standard) emission calculation
Copy( iy_main_agenda, iy_main_agenda__Emission )

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# standard surface agenda (i.e., make use of surface_rtprop_agenda)
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )
#Copy( ppath_step_agenda,ppath_step_agenda__RefractedPath )
#Copy( refr_index_air_agenda, refr_index_air_agenda__NoRefrac)

# blackbody surface with skin temperature interpolated from t_surface field
Copy( surface_rtprop_agenda, surface_rtprop_agenda__Blackbody_SurfTFromt_field )


StringCreate(satellite)
ArrayOfIndexCreate(channels)
ArrayOfIndexCreate(views)
StringCreate(hitran_file)
NumericCreate(f_grid_spacing)


# Select here which satellite you want
# ---
StringSet(satellite, "MET9")


# Select here which channels you want
# ---
#
# MVIRI Channel-1 is WV2 
#              -2 is IR1
#              -3 is VIS
#
# WATCH OUT! We use the usual zero-based ARTS indexing here, so index
# 1 is actually channel 2, etc. 
ArrayOfIndexSet(channels, [3,4,5,6,7,8,9,10,11])

# Select here which views you want
# ---
ArrayOfIndexSet(views, [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25])

# You have to give here the location of the HITRAN catalogue
# ---
# FIXME: This has to be replaced by a little piece of the catalog
# explicitly included here, or?
# We have to discuss with Oliver (and perhaps Patrick) how to handle
# this.
# We can also include a pre-calculated absorption table for HIRS with
# ARTS. 
#StringSet(hitran_file,"/storage3/data/catalogue/hitran/hitran2004/HITRAN04.par")
StringSet(hitran_file,"/scratch/uni/u237/data/catalogue/hitran/hitran2012_140407/HITRAN2012.par")
#StringSet(hitran_file,"/home/patrick/Data/HITRAN_2004/HITRAN04.par")


# Set frequency grid spacing
# (See comments in hirs_reference.arts concerning useful values here)
# ---
StringCreate(f_grid_spacing_str)
NumericSet(f_grid_spacing   , 6e9 ) # default 5e8
#StringSet(f_grid_spacing_str,"6e8")
StringSet(f_grid_spacing_str,"6e9_fast")


# Basic settings (already needed in sensor part)
# ---
# This example assumes 1D
AtmosphereSet1D
# scalar RT
IndexSet( stokes_dim, 1 )


INCLUDE "seviri_fast.arts"
#INCLUDE "seviri_reference.arts"

# Set up absorption
# =================

# Atmospheric profiles
ReadXML( batch_atm_fields_compact, "testdata/garand_profiles.xml.gz" )

# add constant profiles for O2 and N2
batch_atm_fields_compactAddConstant( name="abs_species-O2", value=0.2095 )
batch_atm_fields_compactAddConstant( name="abs_species-N2", value=0.7808 )

# Set parameters for lookup table
# ---
# Arguments omitted for better maintainability of this test file.
abs_lookupSetupBatch


# Optional manual override of T and VMR perturbations
# ---
# If your input data contains extreme outliers, the lookup table will
# get unreasonably large. It is suggested that you instead set them
# manually here. The Chevallier 91L data (humidity set) contains
# temperature variations from -70 to +55 (K) and humidity variations from
# 0 to 6 (fractional units). This should be the reasonable range of
# atmospheric variability. You will get failures from individual jobs in
# the batch, that are outside this range.
#VectorLinSpace( abs_t_pert, -70, 55, 5 )
#VectorLinSpace( abs_nls_pert, 0, 6, 0.5 )


# Create the lookup table
# ---
abs_xsec_agenda_checkedCalc
abs_lookupCalc

#WriteXML("binary", abs_lookup)

# Test (and print) lookup table accuracy
# ---
# This method is not necessary, it just tests and prints the lookup
# table accuracy. Comment it out if you do not want this
# information. The test should take approximately as much time as the
# lookup table generation, perhaps even more. So, it is not cheap!
#abs_lookupTestAccuracy

# Set propmat_clearsky_agenda to use lookup table
# ---
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable )


# Set up RT calculation
# =====================


# Definition of sensor position and LOS
# ---

# Optionally set sensor_pos
# ---
# The sensor altitude is predefined in amsu.arts to 850 km above the geoid. 
# Comment out the next line if you want to set it to something else.
MatrixSetConstant( sensor_pos, 26, 1, 36000e3 )


# Set the agenda for batch calculations:
# ---
#
AgendaSet( ybatch_calc_agenda ){

  
  
  # Extract the atmospheric profiles for this case:
  Extract(
    atm_fields_compact, 
    batch_atm_fields_compact, 
    ybatch_index
  )

  # Split up *atm_fields_compact* to
  # generate p_grid, t_field, z_field, vmr_field:
  AtmFieldsFromCompact

  # Optionally set Jacobian parameters.
  # uncomment this for NO jacobian calculations
  jacobianOff
  # Uncomment this block if you want Jacobians. Attention, it slows down the
  # computation a lot.
  # Also, you can add other Jacobians here, for example for temperature.
#
#   jacobianInit
#   jacobianAddAbsSpecies( jacobian_quantities, jacobian_agenda,
#                          atmosphere_dim,
#                          p_grid, lat_grid, lon_grid, 
#                          p_grid, lat_grid, lon_grid,
#                          "H2O, H2O-SelfContCKDMT100, H2O-ForeignContCKDMT100",
#                          "analytical",
#                          "rel",
#                          0.01 )
#   jacobianClose

  # No scattering
  cloudboxOff

  # get some surface properties from corresponding atmospheric fields
  Extract( z_surface, z_field, 0 )
  Extract( t_surface, t_field, 0 )

  #sensorOff

  # Perform RT calculations
  # ---
  atmfields_checkedCalc(bad_partition_functions_ok=1)
  atmgeom_checkedCalc
  cloudbox_checkedCalc
  sensor_checkedCalc
  yCalc

  # Optionally write out jacobian:
  #  WriteXMLIndexed( output_file_format, ybatch_index, jacobian, "" )

  StringSet( iy_unit, "PlanckBT" )
  yApplyUnit
}


# Set number of batch cases:
nelemGet( ybatch_n, batch_atm_fields_compact )
#IndexSet(ybatch_n, 1)

# Execute the batch calculations:
# ---
propmat_clearsky_agenda_checkedCalc

StringCreate(out_file_name_sensor_response)
StringCreate(out_file_name_ybatch)
StringCreate(out_file_name_f_grid)

StringSet(out_file_name_sensor_response, "TestSEVIRI.sensor_response_")
StringSet(out_file_name_ybatch, "TestSEVIRI.ybatch_")
StringSet(out_file_name_f_grid, "TestSEVIRI.f_grid_")

Append(out_file_name_sensor_response, satellite)
Append(out_file_name_ybatch, satellite)
Append(out_file_name_f_grid, satellite)

StringCreate(extent)
StringSet(extent,"_")
StringSet(dummy,".xml")

Append(extent,f_grid_spacing_str)
Append(extent,dummy)

Append(out_file_name_sensor_response, extent)
Append(out_file_name_ybatch, extent)
Append(out_file_name_f_grid, extent)

WriteXML( "ascii", sensor_response, out_file_name_sensor_response)
ybatchCalc


# Store result matrix:
# ---
WriteXML( "ascii", ybatch , out_file_name_ybatch)
WriteXML( "ascii", f_grid, out_file_name_f_grid)


# Verify results
ArrayOfVectorCreate( ybatch_ref )
ReadXML( ybatch_ref, "TestSEVIRI.ybatch_MET9_6e9_fastREFERENCE.xml" )

Compare( ybatch, ybatch_ref, 0.01,
         "Total BT should be close to the reference values")
}