#DEFINITIONS:  -*-sh-*-
#
# A simple, basic test of 1D and 2D antenna patterns. The setting of *y_geo* is
# also tested inside this file.
#
# 2018-12-12, Patrick Eriksson

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 )
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )

# cosmic background radiation
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )

# on-the-fly absorption
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )

# sensor-only path
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )

# no refraction
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )

# geo-positioning
AgendaSet( geo_pos_agenda ){
  geo_posEndOfPpath
}


# Basic settings
#
AtmosphereSet3D
IndexSet( stokes_dim, 2 )
VectorSet( f_grid, [18e9,31e9] )
StringSet( iy_unit, "PlanckBT" )

# no jacobian calculation
jacobianOff

# no scattering
cloudboxOff

# lat and lon true can be left empty for 3D
VectorSet( lat_true, [] ) 
VectorSet( lon_true, [] ) 


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

# No transitions needed
# 
abs_lines_per_speciesSetEmpty

# Atmospheric grids
#
VectorNLogSpace( p_grid, 41, 1050e2, 100e2 )
VectorNLinSpace( lat_grid, 3, -30, 30 )
VectorNLinSpace( lon_grid, 3, -30, 30 )


# Read a 1D atmospheric case and expand to *atmosphere_dim*
#
AtmRawRead( basename = "testdata/tropical" )
AtmFieldsCalcExpand1D


# Surface = Ocean by FASTEM
#
IndexCreate( nlat )
IndexCreate( nlon )
nelemGet( nlat, lat_grid )
nelemGet( nlon, lon_grid )
#
MatrixSetConstant( z_surface, nlat, nlon, 0 )
#
VectorCreate( trv )
nelemGet( v=f_grid )
VectorSetConstant( trv, nelem, 0.9 )
AgendaSet( surface_rtprop_agenda ){
  specular_losCalc
  InterpAtmFieldToPosition( out=surface_skin_t, field=t_field )
  surfaceFastem( wind_speed=5, wind_direction=45, fastem_version=5,
                 transmittance=trv )
}


# Check data and settings (beside sensor)
# 
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc( bad_partition_functions_ok = 1 )
atmgeom_checkedCalc
cloudbox_checkedCalc


# Sensor pos/los
#
MatrixSet( sensor_pos, [ 800e3, 0, 0;800e3, 0, 0 ] )
MatrixSet( sensor_los, [ 130, 20;130, 15 ] )


# First do without antenna to get reference for geo_pos
#
sensorOff
sensor_checkedCalc
#
yCalc
#
MatrixCreate( geo_ref )
Copy( geo_ref, y_geo )


# Define 1D antenna
#
NumericCreate( xwidth )
NumericSet( xwidth, 3 )
VectorCreate( dza )
VectorCreate( zeros )
IndexCreate( ndlos )
#
antenna_responseVaryingGaussian( antenna_response, 1.5, xwidth, 
                                 0.1, 5, 10e9, 40e9, 0 )
IndexSet( sensor_norm, 1 )
IndexSet( antenna_dim, 1 )
MatrixSet( antenna_dlos, [0] )
#
IndexSet( ndlos, 21 )
VectorNLinSpace( dza, ndlos, -2, 2 )
VectorSetConstant( zeros, ndlos, 0 )
Matrix1ColFromVector( mblock_dlos_grid, dza )


# Calculate
#
sensor_responseInit
sensor_responseAntenna
sensor_checkedCalc
#
yCalc
#
VectorCreate( y_ref )
Copy( y_ref, y )
#
Compare( geo_ref, y_geo, 1e-5 )


# Re-do with 2D antenna (but still "1D" mblock_dlos_grid) 
#
antenna_responseVaryingGaussian( antenna_response, 1.5, xwidth, 
                                 0.1, 5, 10e9, 40e9, 1 )
WriteXML("ascii",antenna_response,"R.xml")
IndexSet( antenna_dim, 2 )


# Calculate
#
sensor_responseInit
sensor_responseAntenna
sensor_checkedCalc
#
yCalc
#
Compare( geo_ref, y_geo, 1e-5 )


# Results for *y* should not be identical, but fairly close
Compare( y, y_ref, 0.001 )


# Repeat with circular mblock_dlos_grid
#
mblock_dlos_gridUniformCircular( spacing=0.2, width=2, centre=1 )
WriteXML("ascii",mblock_dlos_grid,"mdg.xml")
sensor_responseInit
sensor_responseAntenna
sensor_checkedCalc
WriteXML("ascii",sensor_response,"H.xml")
#
yCalc
#
Compare( geo_ref, y_geo, 1e-5 )

# Some deviation expected here for *y*
Compare( y, y_ref, 0.01 )

# reset y_ref
Copy( y_ref, y )


# Repeat with rectangular mblock_dlos_grid
#
mblock_dlos_gridUniformRectangular( spacing=0.2, za_width=2, aa_width=2, centre=1 )
sensor_responseInit
sensor_responseAntenna
sensor_checkedCalc
#
yCalc
#
Compare( geo_ref, y_geo, 1e-5 )

# *y*should basically be identical
Compare( y, y_ref, 1e-6 )

}