#DEFINITIONS:  -*-sh-*-
#
# filename: simpleDISORT.arts
#
# Demonstration of a DISORT scattering calculation
#
# Author: Claudia Emde
# 
# The DISORT calculation is set up in a very similar way to the 
# DOIT calculation. 

Arts {

INCLUDE "general"
INCLUDE "continua.arts"


# Frequency grid 
# --------------
# This example is only monochromatic. Note: The frequency must be contained
# in the gas absorption lookup table.  
IndexSet(nelem){ 1 }
VectorSetConstant(f_grid){
        value = 230e9
}
 

# Pressure grid
ReadXML( p_grid ){
  "./data/p_grid.xml"
}

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

# Atmospheric profiles
AtmRawRead{
 basename = "data/tropical"
}

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

# Lookup tables can be computed using the MATLAB (please contact Stefan)
ReadXML( abs_lookup ){
   "data/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 = -1
}

# Ground altitude (measured from geoid)
nrowsGet(r_geoid)
ncolsGet(r_geoid)
MatrixSetConstant( z_surface){
    value = 500
}

MatrixSetConstant( surface_emissivity_field){
    value = 1
}

basics_checkedCalc



# Definition of sensor position and LOS
#--------------------------------------
# This file holds the viewing angles of the sensor:
#IndexSet(nelem){ 1 }
#VectorSetConstant(vector_1){
#        value = 150
#}

IndexSet(nelem){9}
VectorNLinSpace(vector_1){
	start = 100
	stop = 180
}

# Sensor altitude from earth surface
nelemGet(vector_1){}
VectorSetConstant( vector_2 ) {
   value =  96000.1
}

# Add Earth Radius
VectorAddScalar(vector_2, vector_2){
    value = 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{
      lmax      = -1
  }
}

# Agendas for clearsky radiative transfer calculation
# ----------------------------------------------------
AgendaSet( rte_agenda ){
   RteStd{}
}

# Specification of cloud
# -----------------------
    
# Special DISORT function, because cloudbox needs to be extended over the 
# ehole atmosphere
cloudboxSetDisort
cloudbox_checkedCalc

ParticleTypeInit

# Only one particle type is added in this example
# Note that you must use p20 data, because disort only works with this
# format. 
ParticleTypeAdd{
filename_scat_data="data/scat_data_disort.xml"
filename_pnd_field="data/pnd_field_1D.xml"
}

pnd_fieldCalc
scat_data_rawCheck

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

# For limb calculations is is very important to have a fine resolution     
# about 90°. 
# For DISORT we use the same angular grids as output.
DoitAngularGridsSet{
	N_za_grid = 19
	N_aa_grid = 10
	za_grid_opt_file = ""
	}

# 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
}
	
AgendaSet(iy_cloudbox_agenda){
    Ignore (ppath)
    Ignore (rte_pos)
    Ignore (rte_los)
    iyInterpCloudboxField
    }


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

# Perform scattering calculation 
ScatteringDisort
    
# Calculate RT from cloudbox boundary to the sensor
#
StringSet( y_unit ){ "RJBT" }
#
yCalc


Print(y){
level=1
}

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

} # End of Main