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


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


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

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

# Dimensionality of the atmosphere
AtmosphereSet1D

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

# Definition of species
#abs_speciesSet( species=[ "H2O", "N2", "O2"] )
abs_speciesSet( species=[ "H2O-PWR98",
                          "O2-PWR93",
                          "N2-SelfContStandardType" ] )
#as we use full abs models, we do not need line data. but we need to initialize
# the line data variable.
abs_lines_per_speciesSetEmpty

# Atmospheric profiles
AtmRawRead( basename="testdata/tropical" )
AtmFieldsCalc
atmfields_checkedCalc


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

# in case an abs_lookup table need to be calculated, uncomment the below lines.
#abs_xsec_agenda_checkedCalc
#abs_lookupSetup
#abs_lookupCalc
#WriteXML( output_file_format="binary", in=abs_lookup,
#           filename="gas_abs_lookup.xml" )

ReadXML( abs_lookup, "gas_abs_lookup.xml" )
abs_lookupAdapt

# absorption from LUT
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable )


# Definition of Earth surface
# ----------------------------
MatrixSetConstant( z_surface, 1, 1, 500 )

# Properties of surface:
# - surface reflectivity from Liebe model
# - surface skin temperature interpolated from atmospheric t_field
#
AgendaSet( surface_rtprop_agenda ){
  specular_losCalc
  InterpAtmFieldToPosition( 
    out   = surface_skin_t, 
    field = t_field )
  VectorCreate( n_t_grid )
  VectorSetConstant( n_t_grid, 1, surface_skin_t )
  complex_refr_indexWaterLiebe93( 
    complex_refr_index = surface_complex_refr_index,
    data_f_grid             = f_grid,
    data_T_grid             = n_t_grid )
  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 )

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


# Set the cloudbox limits (pressure units)
# Alternative: cloudboxSetManuallyAltitude (specification in [km])
# ----------------------------------------------------------------
cloudboxSetManually( p1=71617.7922264, p2=17111.6808705,
                     lat1=0, lat2=0, lon1=0, lon2=0 )

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

# Only one particle type is added in this example 
# Here actually both added elements are indentical. however, for testing and for
# demonstration purposed, having 2 elements is better.
ParticleTypeAdd( filename_scat_data="p30f229-231T214-225r100NP-1ar1_5ice.xml",
                 filename_pnd_field="pnd_field_1D.xml" )
ParticleTypeAdd( filename_scat_data="p30f229-231T214-225r100NP-1ar1_5ice.xml",
                 filename_pnd_field="pnd_field_1D.xml" )
scat_data_arrayCheck
pnd_fieldCalc
# rescaling back to pnd equivalent to only one scattering element. in a true
# case, this is NOT needed!
Tensor4Scale( out=pnd_field, in=pnd_field, value=0.5 )

# 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( interp_method="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( N_za_grid=19, N_aa_grid=37,
                     za_grid_opt_file="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_array_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:
  #WriteXMLIndexed( in=doit_i_field, file_index=f_index )
}

# 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( normalize=1,
                           norm_error_threshold=0.05 )
  # 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( epsilon=[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 ){
  iyInterpCloudboxField
}

# No jacobian calculations
#
jacobianOff

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


propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc

# Initialize Doit variables
DoitInit

# Calculate incoming radiation field at cloudbox boundary 
CloudboxGetIncoming

# Executes doit_mono_agenda for all frequencies
ScatteringDoit
# Writing out the radition field solution from inside cloudbox for
# initialization of perturbed pnd_field calculation from non-clearsky first
# guess
WriteXML( in=doit_i_field1D_spectrum )

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

WriteXML( in=y )

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

#==================check==========================

VectorCreate(yREFERENCE)
ReadXML( yREFERENCE, "yREFERENCE_DOIT.xml" )
Compare( y, yREFERENCE, 2e-5 )

} # End of Main