# example: limb case using an optimized polar angle grid.
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 )

# No jacobian calculations
#
jacobianOff


# Frequency grid 
# --------------
# Note: The frequencies must be contained in the gas absorption lookup table.
VectorSet( f_grid, [229.5e9,230.5e9] )

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

# Radiance output units
#----------------------
StringSet( iy_unit, "RJBT" )


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

# Dimensionality of the atmosphere
AtmosphereSet1D

# Pressure grid
ReadXML( p_grid, "testdata/testdoit_p_grid.xml" )

# Definition of species
abs_speciesSet( species=[ "H2O-PWR98",
                          "O2-PWR93",
                          "N2-SelfContStandardType" ] )

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


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

# in case an abs_lookup table need to be calculated, uncomment the block below
# (and replace filename for abs_lookup table to read in ReadXML call).
#atmfields_checkedCalc
#abs_lines_per_speciesSetEmpty
#abs_xsec_agenda_checkedCalc
#abs_lookupSetup
#abs_lookupCalc
#WriteXML( output_file_format="binary", in=abs_lookup, filename="my_gas_abs_lookup.xml" )

ReadXML( abs_lookup, "testdata/testdoit_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 )

# Emission and reflection properties:
# - specular reflecting water surface at ambient atmospheric temperature
# - water reflective indices from Liebe model
# - surface skin temperature interpolated from atmospheric t_field
#
VectorCreate( n_t_grid )
AgendaSet( surface_rtprop_agenda ){
  specular_losCalc
  InterpAtmFieldToPosition( 
    out   = surface_skin_t, 
    field = t_field )
  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
#--------------------------------------
VectorCreate( vector_1 )
VectorCreate( vector_2 )

# Sensor viewing angles
VectorSet( vector_1, [99.7841941981] )
#VectorNLinSpace( vector_2, 4, 120, 180 ) #extend with a couple of downviewing los
#Append( vector_1, vector_2 )
nelemGet( nelem, vector_1 )

# Sensor altitude from earth surface
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


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

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

# Initialization and adding scattering elements.
# ----------------------------------------------------------------
ScatSpeciesInit

# Adding scattering elements.
# ----------------------------------------------------------------
# Here actually both added elements are indentical. however, for testing and for
# demonstration purposed, having 2 elements is better.
ScatElementsPndAndScatAdd(
  scat_data_files=["testdata/scatData/azi-random_f229-231T214-225r100NP-1ar1_5ice.xml",
                   "testdata/scatData/azi-random_f229-231T214-225r100NP-1ar1_5ice.xml"],
  pnd_field_files=["testdata/testdoit_pnd_field_1D.xml",
                   "testdata/testdoit_pnd_field_1D.xml"] )

#scat_dataCheck( scat_data = scat_data_raw )
scat_dataCalc
#scat_dataCheck

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



# -------------
# DOIT settings
# -------------

# 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. 
DOAngularGridsSet( N_za_grid=19, N_aa_grid=37,
                     za_grid_opt_file="testdata/testdoit_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
  Ignore( f_grid )
  # Alternative method (needs less memory):
  #scat_data_monoCalc
  # Perform iterations: 1. scattering integral. 2. RT calculations with 
  # fixed scattering integral field, 3. convergence test 
  doit_i_field_monoIterate
  # Write the radiation field inside the cloudbox:
  #WriteXMLIndexed( in=doit_i_field_mono, 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
}


# 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
  #Print( doit_iteration_counter, 0 )
}

AgendaSet( iy_cloudbox_agenda ){
  iyInterpCloudboxField
}

} # End of Main