#DEFINITIONS:  -*-sh-*-
#
# ARTS control file for testing 3D propagation path calculations
# 
# It is also demonstration of how to use the ForLopp agenda.
#
# 2012-02-17 Patrick Eriksson

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)

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )


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

# A pressure grid rougly matching 0 to 80 km.
#
VectorNLogSpace( p_grid, 41, 1000e2, 1 )


# Atmospheric dimensionality and lat/lon grids
#
VectorNLinSpace( lat_grid, 21, 35, 55 )
VectorNLinSpace( lon_grid, 21, -40, 40 )
AtmosphereSet3D


# Water vapour needed if refraction will be calculated
#
abs_speciesSet( species=["H2O"] )


# Read a 1D atmospheric case and expand to *atmosphere_dim*
#
AtmRawRead( t_field_raw, z_field_raw, vmr_field_raw, abs_species, 
            "testdata/tropical" )
AtmFieldsCalcExpand1D


# Ground altitude (z_surface)
#
IndexCreate( nlat )
IndexCreate( nlon )
nelemGet( nlat, lat_grid )
nelemGet( nlon, lon_grid )
#
MatrixSetConstant( z_surface, nlat, nlon, 500 )


# Initializing cloudbox: No scattering
#
cloudboxOff

# Activate to make tests with cloudbox on:
#FlagOn( cloudbox_on )
#ArrayOfIndexSet( cloudbox_limits, [ 4, 7, 10, 11, 8, 12 ] )


# A dummy frequency grid
#
VectorSet( f_grid, [10e9] )


# Check if atmosphere OK
#
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc


# Max step length for the representation of the propagation path
#
NumericSet( ppath_lmax, 20e3 )




#
# A single propagation path:
#

# Set a observation position and line-of-sight (LOS)
#
VectorSet( rte_pos, [ 600e3, 70.1, 0 ] )
VectorSet( rte_los, [   113,  180    ] )
VectorSet( rte_pos2, [] )                  # No transmitter involved


#
# no refraction
#
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

# Calculate propagation path
#
ppathCalc                       


# Print complete ppath
#
#Print( ppath, 0 )



# Uncomment to skip batch part below
# 
#Exit()



#
# Run through a number of cases, that should run without error
# 


# Use sensor_pos/los to store each case
#
MatrixSet( sensor_pos, [ 600e3, 20.2,   23;
                         600e3, 70.3,  -12;
                         600e3, 20.1,  -30; 
                         600e3, 70.2,   11;
                         600e3,   41,  1.5;
                         600e3, 37.4,  -20;
                         600e3,   50,   10;
                           6e3, 40.1, 34.6;
                           6e3, 40.1,  -39;
                           6e3,   41,  1.5;
                           6e3, 37.4,  -20;
                           6e3,   50,   10;
                           500,   41, 0.23;
                           500, 37.6, -1.1 ] )
MatrixSet( sensor_los, [   45,    15; 
                           95,   135; 
                          113,    12; 
                          113,  -168;
                          180,     0;
                          180,   156;
                          180,   -29;
                           34,  -160;
                          156,    90;
                            0,     0;
                            0,   156;
                            0,   -29;
                           10,     0; 
                          112,    55 ] )

IndexCreate( ilast )
nrowsGet( ilast, sensor_pos )
IndexStepDown( ilast, ilast )

AgendaSet( forloop_agenda ){
  #Print( forloop_index, 0 )
  VectorExtractFromMatrix( rte_pos, sensor_pos, forloop_index, "row" )
  VectorExtractFromMatrix( rte_los, sensor_los, forloop_index, "row" )
  ppathCalc
}

#
# no refraction
#
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

ForLoop( forloop_agenda, 0, ilast, 1  )





#
# Repeat with refraction
#
Copy( ppath_step_agenda, ppath_step_agenda__RefractedPath )
NumericSet( ppath_lraytrace, 1e3 )
Copy( refr_index_air_agenda, refr_index_air_agenda__GasThayer  )

ForLoop( forloop_agenda, 0, ilast, 1  )





#
# Repeat all with non-spherical ellipsoid
#

refellipsoidEarth( refellipsoid, "WGS84" )

Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )
ForLoop( forloop_agenda, 0, ilast, 1  )

Copy( ppath_step_agenda, ppath_step_agenda__RefractedPath )
Copy( refr_index_air_agenda, refr_index_air_agenda__GasThayer  )
ForLoop( forloop_agenda, 0, ilast, 1  )

}