#DEFINITIONS:  -*-sh-*-

# --------------------------------------------------------------------
# Test batch calculations for HIRS.
# 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 )

# Planck as blackbody radiation
Copy( blackbody_radiation_agenda, blackbody_radiation_agenda__Planck )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

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


# Select here which channels you want
# ---
#
# HIRS Channels 13-19 are shortwave channels. Simulating them with
# ARTS for thermal radiation only is pointless.
#
# WATCH OUT! We use the usual zero-based ARTS indexing here, so index
# 11 is actually channel 12, etc. 
#ArrayOfIndexSet(channels, [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18])
ArrayOfIndexSet(channels, [0,1,2,3,4,5,6,7,8,9,10,11])
#ArrayOfIndexSet(channels, [11])

# Select here which views you want
# ---
ArrayOfIndexSet(views, [0, 7, 14, 21, 27])


# 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,"/home/patrick/Data/HITRAN_2004/HITRAN04.par")


# Set frequency grid spacing
# (See comments in hirs_reference.arts concerning useful values here)
# ---
NumericSet(f_grid_spacing, 5e8)


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


INCLUDE "hirs_fast.arts"


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

# Atmospheric profiles
# ---
# Atmospheric profiles are stored in an ArrayOfMatrix.
# It contains one matrix for each atmospheric state.
# Each matrix row corresponds to one pressure level. The
# meaning of the columns is:
# p[Pa]        T[K]         z[m]         H2O[VMR]	O3[VMR]
#
ArrayOfMatrixCreate( arrayofmatrix_1 )
#ReadXML( arrayofmatrix_1,
#         "/storage2/home/sbuehler/checkouts/arts-xml-data/atmosphere/chevallier_91L/chevallierl91_clear_q.xml" )
#ReadXML( arrayofmatrix_1, "../AMSU/chevallierl91_clear_q_extract.xml" )

ReadXML( arrayofmatrix_1, "testdata/garand_profiles.xml.gz" )

# Storage in an array of matrix is handy for profile data, because
# it is very compact. However, the more general internal
# representation of the data is in batch_atm_fields_compact.

# Convert to batch_atm_fields_compact
# ---
# The values taken for O2 and N2 are from Wallace&Hobbs, 2nd edition.
batch_atm_fields_compactFromArrayOfMatrix(
   batch_atm_fields_compact,
   atmosphere_dim,
   arrayofmatrix_1,
   ["T", "z", "H2O", "O3", "CO2", "N2O", "CO", "CH4"],
   ["O2", "N2"],
   [0.2095, 0.7808]
)

# Delete original data array to conserve memory:
# ---
Delete( arrayofmatrix_1 )

# Set parameters for lookup table
# ---
# Arguments omitted for better maintainability of this test file.
abs_lookupSetupBatch
# abs_lookupSetupBatch( abs_p, 
#                       abs_t, 
#                       abs_t_pert, 
#                       abs_vmrs, 
#                       abs_nls, 
#                       abs_nls_pert, 
#                       abs_species, 
#                       batch_atm_fields_compact, 
#                       abs_t_interp_order, 
#                       abs_nls_interp_order, 
#                       , 
#                       , 
#                       , 
#                       )


# 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, 1, 1, 850e3 )


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

  # Perform RT calculations
  # ---
  atmfields_checkedCalc
  atmgeom_checkedCalc
  cloudbox_checkedCalc
  sensor_checkedCalc
  yCalc

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

  # Convert the measurement from radiance units to Planck Tb:
  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
ybatchCalc


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


}