#DEFINITIONS:  -*-sh-*-
#
# filename: TestRT4.arts
#
# Demonstration of a RT4 scattering calculation
#
# The setup follows the one used by TestDisort (in disort/ folder) and
# TestDOITFromIndividualFields.arts (in doit/ folder).
#
# Author: Jana Mendrok
# 

Arts2 {

# for test case, first, run the proprietary RT4 test.
# REMOVE (or outcomment) when you use this as template for your RT4-based
# calculations.
INCLUDE "TestRT4Test.arts"


# here starts the the actual ARTS RT4 interface test/example
INCLUDE "artscomponents/disort/indivfieldsatmo_setup.arts"

# the "old way": use RT4's own Lambertian surface methods (uncomment the
# following three call blocks)
#   first, touch all explicit-input surface reflection variables, such that they
#   exist. fill the ones required for your choice of ground type below (see
#   online doc for what choices are available and what input is required for
#   each of them).
#Touch(surface_complex_refr_index)
#Touch(surface_reflectivity)
#Touch(surface_scalar_reflectivity)

# Set surface skin temperature based on t_field
#Extract( t_surface, t_field, 0 )
#VectorCreate( tskin_vector )
#VectorExtractFromMatrix( tskin_vector, t_surface, 0, "column" )
#Extract( surface_skin_t, tskin_vector, 0 )

# Set reflection type
#StringCreate( ground_type )
#StringSet( ground_type, "L" )


# Set cloudbox limits 
# ---------------------------------------------------------------------
#cloudboxSetFullAtm
cloudboxSetAutomatically( particle_field=scat_species_mass_density_field )
cloudboxSetAutomatically( particle_field=scat_species_mass_flux_field,
                          cloudbox_limits_old=cloudbox_limits )
cloudboxSetAutomatically( particle_field=scat_species_number_density_field,
                          cloudbox_limits_old=cloudbox_limits )
cloudboxSetAutomatically( particle_field=scat_species_mean_mass_field,
                          cloudbox_limits_old=cloudbox_limits )

# Particle Number Density field calculation
# -----------------------------------------
pnd_fieldCalcFromscat_speciesFields



# 7. RT4 settings-----------------------------------------------------
#------------------------------------------------------------------------

# For RT4 we use the same angular grids as output.
#DOAngularGridsSet( N_za_grid=38 )

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

AgendaSet( iy_cloudbox_agenda ){
  iyInterpCloudboxField
  #iyInterpCloudboxField( za_inter_order=3, za_restrict=1 )
}


# 8. The RT calculation--------------------------------------------------
#------------------------------------------------------------------------

# Consistency checks
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc  
sensor_checkedCalc
scat_data_checkedCalc

# Perform scattering calculation
RT4Calc #( nstreams=64 )
# alternatively, the old way using one of RT4's own surface methods:
#RT4CalcWithRT4Surface( ground_type=ground_type ) #( nstreams=32 )
#WriteXML( output_file_format, doit_i_field )

Tensor7Create( ifield_old_optprop )
Copy( ifield_old_optprop, doit_i_field )
RT4Calc( new_optprop=1 ) #, nstreams=64 )
CompareRelative( ifield_old_optprop, doit_i_field, 1e-6 )
Compare( ifield_old_optprop, doit_i_field, 1e-30 )

# Calculate RT from cloudbox boundary to the sensor
yCalc

#WriteXML( output_file_format, y )

# Verify results
VectorCreate(yREFERENCE)
ReadXML( yREFERENCE, "artscomponents/disort/yREFERENCEfromDOIT.y.xml" )
Compare( y, yREFERENCE, 1. )

} # End of Main