#DEFINITIONS:  -*-sh-*-
#
# filename: simpleDOIT.arts
#
# Demonstration of a DOIT scattering calculation
#
# Author: Claudia Emde
# 
# A big part of the control file consists of definitions for the 
# radiative transfer calculations. Parts, where a calculation is 
# performed are marked with ===start=== and ===stop===.
#

Arts2 {

INCLUDE "general.arts"
INCLUDE "continua.arts"


# Frequency grid 
# --------------
# This example is only monochromatic. Note: The frequency must be contained
# in the gas absorption lookup table.  
VectorSet( f_grid, [230e9] )


# Number of Stokes components to be computed
#-------------------------------------------
IndexSet( stokes_dim, 4 )

# Definition of the atmosphere
# ----------------------------

# Pressure grid
ReadXML( p_grid, "p_grid.xml" )

# Definition of species
SpeciesSet( abs_species, [ "H2O", "N2", "O2"] )

# Atmospheric profiles
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, "../atmosphere_data/tropical" )

# Gas absorption from lookup table 
# ---------------------------------

# Lookup tables can be computed using the MATLAB (please contact Stefan)
ReadXML( abs_lookup, "gas_abs_lookup.xml" )
   	
abs_lookupAdapt

AgendaSet( abs_scalar_gas_agenda ){ 
  abs_scalar_gasExtractFromLookup 
}

AtmFieldsCalc


# Definition of Earth surface
# ----------------------------
# spherical geoid 
#r_geoidWGS84

r_geoidSpherical( r_geoid, atmosphere_dim, lat_grid, lon_grid, -1 )

# Ground altitude (measured from geoid)
nrowsGet( nrows, r_geoid )
ncolsGet( ncols, r_geoid )
MatrixSetConstant( z_surface, nrows, ncols, 500 )

# Properties of surface

AgendaSet( surface_prop_agenda ){
  Ignore( rte_pos )
  InterpAtmFieldToRteGps( surface_skin_t, atmosphere_dim,
                          rte_gp_p, rte_gp_lat, rte_gp_lon, t_field )
  complex_nWaterLiebe93( complex_n, f_grid, surface_skin_t )
  surfaceFlatRefractiveIndex
}


# Definition of sensor position and LOS
#--------------------------------------
# This file holds the viewing angles of the sensor:
IndexSet( nelem, 1 )
VectorCreate( vector_1 )
VectorSetConstant( vector_1, nelem, 99.7841941981 )

#IndexSet( nelem, 19 )
#VectorNLinSpace( vector_1, 0, 180 )

# Sensor altitude from earth surface
nelemGet( nelem, vector_1 )
VectorCreate( vector_2 )
VectorSetConstant( vector_2, nelem, 95000.1 )

# Add Earth Radius
VectorAddScalar( vector_2, vector_2, 6.378e6 )

Matrix1ColFromVector( sensor_pos, vector_2 )
Matrix1ColFromVector( sensor_los, vector_1 )

# SensorOff means that the result of the calculation are the radiances,
# which are not modified by sensor properties
sensorOff

# This should be sufficiently small to assume single scattering in 
# one propagation path step
AgendaSet( ppath_step_agenda ){
  ppath_stepGeometric( ppath_step, atmosphere_dim, p_grid, lat_grid, lon_grid,
                       z_field, r_geoid, z_surface, -1 )
}

# Set the cloudbox limits (pressure units)
# Alternative: cloudboxSetManuallyAltitude (specification in [km])
# ----------------------------------------------------------------
cloudboxSetManually( cloudbox_on, cloudbox_limits, atmosphere_dim, p_grid,
                     lat_grid, lon_grid, 71617.7922264, 17111.6808705,
                     0, 0, 0, 0 )

# Specification of cloud
# -----------------------
ParticleTypeInit

# Only one particle type is added in this example 
ParticleTypeAdd( scat_data_raw, pnd_field_raw, atmosphere_dim, f_grid, p_grid,
                 lat_grid, lon_grid, cloudbox_limits,
                 "p30f229-231T214-225r100NP-1ar1_5ice.xml",
                 "pnd_field_1D.xml" )

pnd_fieldCalc


# Select interpolation method ('linear' or 'polynomial'):
# ----------------------------------------------------

# For limb calculations is is very important to have a fine resolution     
# about 90°.
# If "polynomial" is selected one has to use an optimized grid. Please     
# use *doit_za_grid_optCalc* to optimize the grid.
doit_za_interpSet( doit_za_interp, atmosphere_dim, "linear" )

# As one needs for the RT calculations in limb direction a much finer
# zenith angle grid resolution as for the computation of the scattering 
# integral, it is necessary to define two zenith angle grids. For the scattering 
# integral equidistant grids are created from N_za_grid, N_aa_grid, which give 
# the number of grid points. An optimized grid for the RT calculation needs to 
# be given by a file. If no filename is specified (za_grid_opt_file = ""), the 
# equidistant grids are used for both, scattering integral and RT calculation. 
# This option can only be used for down-looking geometries. 
DoitAngularGridsSet( doit_za_grid_size, scat_aa_grid, scat_za_grid,
                     19, 10, "doit_za_grid_opt.xml" )

# Main agenda for DOIT calculation
# --------------------------------
#
# Input: incoming field on the cloudbox boundary
# Ouput: the scattered field on the cloudbox boundary
AgendaSet( doit_mono_agenda ){
  # Prepare scattering data for DOIT calculation (Optimized method):
  DoitScatteringDataPrepare
  # Alternative method (needs less memory):
  #scat_data_monoCalc
  # Set first guess field:		
  doit_i_fieldSetClearsky
  # Perform iterations: 1. scattering integral. 2. RT calculations with 
  # fixed scattering integral field, 3. convergence test 
  doit_i_fieldIterate
  # Put solution into interface for clearsky calculation   
  DoitCloudboxFieldPut
  # Write the radiation field inside the cloudbox:
  WriteXML( output_file_format, doit_i_field )
}

# Definitions for methods used in *i_fieldIterate*:
#----------------------------------------------------

# 1. Scattering integral
# --------------------------

# Calculation of the phase matrix
AgendaSet( pha_mat_spt_agenda ){
  # Optimized option:
  pha_mat_sptFromDataDOITOpt
  # Alternative option:
  #pha_mat_sptFromMonoData
}

AgendaSet( doit_scat_field_agenda ){
  doit_scat_fieldCalcLimb
  # Alternative: use the same za grids in RT part and scattering integral part
  #doit_scat_fieldCalc
}

# 2. Radiative transfer with fixed scattering integral term
# ---------------------------------------------------------

AgendaSet( doit_rte_agenda ){
  # Sequential update for 1D
  doit_i_fieldUpdateSeq1D
  # Alternatives:
  # Plane parallel approximation for determination of propagation path steps
  #doit_i_fieldUpdateSeq1DPP
  # Without sequential update (not efficient):
  #doit_i_fieldUpdate1D
  # 3D atmosphere:
  #doit_i_fieldUpdateSeq3D
}

# Calculate opticle properties of particles and add particle absorption
# and extiction to the gaseous properties to get total extinction and
# absorption:

AgendaSet( spt_calc_agenda ){
  opt_prop_sptFromMonoData
}

AgendaSet( opt_prop_part_agenda ){
  ext_matInit
  abs_vecInit
  ext_matAddPart
  abs_vecAddPart
}
	

# 3. Convergence test
# ----------------------
AgendaSet( doit_conv_test_agenda ){
  # Give limits for all Stokes components in Rayleigh Jeans BT:
  doit_conv_flagAbsBT( doit_conv_flag, doit_iteration_counter, doit_i_field,
                       doit_i_field_old, f_grid, f_index,
                       [0.1, 0.01, 0.01, 0.01] )

  # Alternative: Give limits in radiances
  #doit_conv_flagAbs( doit_conv_flag, doit_iteration_counter, doit_i_field,
  #                   doit_i_field_old ){
  #  epsilon = [0.1e-15, 0.1e-18, 0.1e-18, 0.1e-18]
  #}
  # If you want to investigat several iteration fields, for example 
  # to investigate the convergence behavior, you can use
  # the following method:
  #DoitWriteIterationFields( doit_iteration_counter, doit_i_field ){
  #  iterations = [2, 4]
  #}
}

AgendaSet( iy_cloudbox_agenda ){
  Ignore (rte_pos)
  Ignore (rte_los)
  iyInterpCloudboxField
}

#==================start==========================

# Perform scattering calculation 
AgendaSet( iy_clearsky_agenda ){
  iyEmissionStandardClearsky
}
basics_checkedCalc
cloudbox_checkedCalc
CloudboxGetIncoming

# Initialize Doit variables
DoitInit

# Executes doit_mono_agenda for all frequencies
ScatteringDoit

# Calculate RT from cloudbox boundary to the sensor
#
StringSet( y_unit, "RJBT" )
#
yCalc

WriteXML( "ascii", y )
WriteXML( "ascii", y_aux )
WriteXML( "ascii", y_error )



#==================stop==========================

} # End of Main