#DEFINITIONS:  -*-sh-*-
#This control file performs an ARTS-MC radiative transfer simulation
#For a single line of sight and a Gaussian antenna response


Arts2 {

INCLUDE "general/general.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/planet_earth.arts"

# Agenda for scalar gas absorption calculation
Copy(abs_xsec_agenda, abs_xsec_agenda__noCIA)

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# Planck as blackbody radiation
Copy( blackbody_radiation_agenda, blackbody_radiation_agenda__Planck )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

# blackbody surface with skin temperature interpolated from t_surface field
Copy( surface_rtprop_agenda, surface_rtprop_agenda__Blackbody_SurfTFromt_field )


#### LOAD DATA: these files were created with MCDataPrepare.arts ######

ReadXML( f_grid, "TestMonteCarloDataPrepare.f_grid.xml" )

IndexSet( f_index, 0 )

ReadXML( p_grid, "p_grid.xml" )

AtmosphereSet3D

ReadXML( lat_grid, "lat_grid.xml" )

ReadXML( lon_grid, "lon_grid.xml" )

ReadXML( t_field, "TestMonteCarloDataPrepare.t_field.xml" )

ReadXML( z_field, "TestMonteCarloDataPrepare.z_field.xml" )

ReadXML( vmr_field, "TestMonteCarloDataPrepare.vmr_field.xml" )

ReadXML( z_surface, "TestMonteCarloDataPrepare.z_surface.xml" )

ReadXML( abs_lookup, "TestMonteCarloDataPrepare.abs_lookup.xml" )

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

abs_lookupAdapt

FlagOn( cloudbox_on )
ReadXML( cloudbox_limits, "TestMonteCarloDataPrepare.cloudbox_limits.xml" )

ReadXML( pnd_field, "TestMonteCarloDataPrepare.pnd_field.xml" )

ReadXML( scat_data_array_mono, "TestMonteCarloDataPrepare.scat_data_array_mono.xml" )


#### Define Agendas #################################################

# absorption from LUT
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__LookUpTable )


#### Define viewing position and line of sight #########################

rte_losSet( rte_los, atmosphere_dim, 99.7841941981, 180 )

rte_posSet( rte_pos, atmosphere_dim, 95000.1, 7.61968838781, 0 )

Matrix1RowFromVector( sensor_pos, rte_pos )

Matrix1RowFromVector( sensor_los, rte_los )


#### Set some Monte Carlo parameters ###################################

IndexSet( stokes_dim, 4 )

StringSet( iy_unit, "RJBT" )

NumericSet( ppath_lmax, 3e3 )

MCSetSeedFromTime

Print( mc_seed, 1 )

mc_antennaSetGaussianByFWHM( mc_antenna, 0.1137, 0.239 )


#### Check atmosphere ##################################################

atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc( scat_data_array=scat_data_array_mono )


#### Perform Monte Carlo RT Calculation #################################

NumericSet( mc_std_err, -1 )
IndexSet( mc_max_time, 20 )
IndexSet( mc_max_iter, -1 )

abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
MCGeneral

#### Save calculated Stokes vector and std. err.#########################

output_file_formatSetAscii
WriteXML( output_file_format, y )

WriteXML( output_file_format, mc_error )

#### Print number of photons and radiance units ########################

Print( mc_iteration_count, 1 )

Print( iy_unit, 1 )

#### Tests ########################

# Radiance test
NumericCreate( mc_error_0 )
Extract( mc_error_0, mc_error, 0 )
NumericScale( mc_error_0, mc_error_0, 4. )

VectorCreate( y_0 )
Select( y_0, y, [ 0 ] )

VectorCreate( y_ref )
VectorSetConstant( y_ref, 1, 198.6 )

Compare( y_0, y_ref, mc_error_0,
         "Total radiance should be close to 198.6 K" )

# Polarization test
NumericCreate( mc_error_1 )
Extract( mc_error_1, mc_error, 1 )
NumericScale( mc_error_1, mc_error_1, 4. )

VectorCreate( y_1 )
Select( y_1, y, [ 1 ] )

VectorSetConstant( y_ref, 1, 7.6 )

Compare( y_1, y_ref, mc_error_1,
         "Polarization difference should be close to 7.6 K" )
}