# --------------------------------------------------------------------
# Test batch calculations for AVHRR. See also TestHIRS.arts.
# Gerrit Holl, July 2011
# Based on TestHIRS.arts
# --------------------------------------------------------------------

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, "NOAA19")


# Select here which channels you want
# ---
#
# ARTS 0 --> AVHRR 3B
# ARTS 1 --> AVHRR 4
# ARTS 2 --> AVHRR 5
ArrayOfIndexSet(channels, [1,2])

# Select here which views you want.
# AVHRR LAC/FRAC has 1024 views
# ---
ArrayOfIndexSet(views, [0, 256, 512, 768, 1023])

StringSet(hitran_file,"/storage3/data/catalogue/hitran/hitran2004/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 "avhrr_reference.arts"


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

ArrayOfMatrixCreate( arrayofmatrix_1 )

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

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

# Setup lookup table
# ---
abs_lookupSetupBatch
abs_lookupCalc

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

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

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

  AtmFieldsFromCompact

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

  # No scattering
  cloudboxOff

  # No jacobian calculations
  jacobianOff

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

  # Convert the measurement from radiance units to Planck Tb:
  StringSet( iy_unit, "PlanckBT" )
  yApplyUnit
}

# Set number of batch cases (uncomment according to your needs):
# this one for all cases in atmosphere batch
#nelemGet( ybatch_n, batch_atm_fields_compact )
# for testing, ALL cases take too long. so we do only the first 2 (still takes
# ~15min on an 8-node system.
IndexSet(ybatch_n, 2)

# Execute the batch calculations:
# ---
ybatchCalc

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


}