#DEFINITIONS:  -*-sh-*-
#
# The script demonstrates a 3D radio link calculation. The settings mimics a
# satellite-to-ground link, but by changing the receiver position a radio
# occultation geometry can be ac hived. Calculations are done for three
# frequencies (f_centre+-f_bwidth/2), where dispersion is considered.
#
# The basic atmosphere is Fascode tropical. An ionosphere is added. By
# changing an include file, a non-ionosphere case can be obtained (resulting in
# faster calculations).
#
# Run with report level 0: 
#   arts -r0 link_budget.arts
# to obtain formatted screen output with results.
#
# Author: Patrick Eriksson

Arts2{

#============================================================================
# User settings
#============================================================================

#
# Position of receiver 
#
VectorCreate( reciever )      # z [m], lat [deg], lon  [deg]  
#
VectorSet( reciever, [ 0e3, 15, 78 ] )         

#
# Position of transmitter 
#
# (If propagation close to the poles, modifications below can be needed.) 
#
VectorCreate( transmitter )   # z [m], lat [deg], lon  [deg]  
#
VectorSet( transmitter, [ 800e3, 17, 88 ] )

#
# Frequency
#
NumericCreate( f_centre )     # Centre frequency [Hz]
NumericCreate( f_bwidth )     # Bandwidth (full) [Hz]
#
NumericSet( f_centre, 2e9 )  
NumericSet( f_bwidth, 0.1e9 ) 

#
# Transmitted polarisation
#
# Do 'arts -d instrument_pol' for information on the coding of polarisation states.
# For example, V and H are coded as 5 and 5, respectively.
#
ArrayOfIndexSet( instrument_pol, [ 5 ] )


#
# Partition function handling
#
IndexCreate( bad_partition_functions_ok )
# ok to use extrapolation of partition functions (i.e. fitting range of
#  partition function paameterization does not cover the actual atmospheric
#  temperatures)?
# ( 0: not ok, 1: ok )
IndexSet( bad_partition_functions_ok, 0 )




#============================================================================
# Basic and RT settings
#============================================================================

#
# General initialisation
#
INCLUDE "general/general.arts"
INCLUDE "general/continua.arts"
INCLUDE "general/agendas.arts"

# 
# Propgation path agendas
#
Copy( ppath_agenda,            ppath_agenda__TransmitterReceiverPath )
Copy( ppath_step_agenda,       ppath_step_agenda__RefractedPath      )
Copy( refr_index_air_agenda,   refr_index_air_agenda__GasMicrowavesEarth  )

#
# Radiative transfer step lengths
#
NumericSet( ppath_lmax, 10e3 )
NumericSet( ppath_lraytrace, 500 )

#
# Here we consider dispersion (overkill if not ionopshere included)
#
Copy( iy_main_agenda, iy_main_agenda__Freqloop )
Copy( iy_sub_agenda, iy_sub_agenda__Radiolink )


#
# Transmitted signal
#
Copy( iy_transmitter_agenda, iy_transmitter_agenda__UnitPolIntensity )


#
# f_grid
#
NumericCreate( grid_hwidth )
NumericCreate( grid_start  )
NumericCreate( grid_stop   )
#
NumericScale( grid_hwidth, f_bwidth,    0.5         )
NumericAdd(   grid_stop,   f_centre,    grid_hwidth )
NumericScale( grid_hwidth, grid_hwidth, -1          )
NumericAdd(   grid_start,  f_centre,    grid_hwidth )
#
VectorNLinSpace( f_grid, 3, grid_start, grid_stop   )

# Postion and line-of-sight of sensor
# 
VectorSet( rte_los, [] )   # Dummy value
#
Copy( rte_pos, reciever )
Copy( rte_pos2, transmitter )

#
# Basic stuff
#
IndexSet( stokes_dim, 4 )
AtmosphereSet3D
#
jacobianOff
cloudboxOff
sensorOff
#
VectorSet( lat_true, [] )
VectorSet( lon_true, [] )





#============================================================================
# Generate lat_grid and lon_grid
#============================================================================

#
# Both grids: 21 points covering = +-10 degs around receiver
#
NumericCreate( grid0 )
IndexCreate( grid_length )
#
NumericSet( grid_hwidth, 15 )   # Half-width [deg]
IndexSet(   grid_length, 21 )   # Total number of grid points
#
Extract(      grid0,       reciever,    1           )
NumericAdd(   grid_stop,   grid0,       grid_hwidth )
NumericScale( grid_hwidth, grid_hwidth, -1          )
NumericAdd(   grid_start,  grid0,       grid_hwidth )
#
VectorNLinSpace( lat_grid, grid_length, grid_start, grid_stop )
#
Extract(      grid0,       reciever,    2           )
NumericAdd(   grid_start,  grid0,       grid_hwidth )
NumericScale( grid_hwidth, grid_hwidth, -1          )
NumericAdd(   grid_stop,   grid0,       grid_hwidth )
#
VectorNLinSpace( lon_grid, grid_length, grid_start, grid_stop )




#============================================================================
# Load atmosphere and absorption lines
#============================================================================

#INCLUDE "earth_tropical_false3d_0to64km_lowfreq.arts"
INCLUDE "earth_tropical_false3d_0to1000km_lowfreq.arts"





#============================================================================
# Calculate and output
#============================================================================

#
# Run data checks
#
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc( bad_partition_functions_ok=bad_partition_functions_ok )
atmgeom_checkedCalc
cloudbox_checkedCalc


#
# Auxilary variables
#
ArrayOfStringSet( iy_aux_vars, [ "Faraday rotation", "Extra path delay" ] ) 

#
# Determine transmitted signal (just for report purpose)
#
MatrixCreate( iy_trans )
#
VectorSet( rtp_pos, [] )  # Dummy values
VectorSet( rtp_los, [] )
AgendaExecute( iy_transmitter_agenda )
Copy( iy_trans, iy )


#
# Perform RT
#
iyCalc

#
# Extract total Faraday rotation
#
Tensor4Create( farrot_total )
Extract( farrot_total, iy_aux, 0 )


#
# Report
#
StringCreate( message )
StringCreate( seperator )
StringCreate( empty_line )
StringSet( seperator, "-" )
StringSet( empty_line, "" )
#
StringSet( message, "---- Result of radio link calculation ---" )
Print( empty_line, 0 )
Print( message, 0 )
#
StringSet( message, "Frequencies [GHz]:" )
VectorCreate( f_ghz )
VectorScale( f_ghz, f_grid, 1e-9 )
Print( empty_line, 0 )
Print( message, 0 )
Print( seperator, 0 )
Print( f_ghz, 0 )
#
StringSet( message, "Transmitted Stokes vector (each row a frequency):" )
Print( empty_line, 0 )
Print( message, 0 )
Print( seperator, 0 )
Print( iy_trans, 0 )
#
StringSet( message, "Received Stokes vector:" )
Print( empty_line, 0 )
Print( message, 0 )
Print( seperator, 0 )
Print( iy, 0 )
#
StringSet( message, "Total Faraday rotation [deg]:" )
Print( empty_line, 0 )
Print( message, 0 )
Print( seperator, 0 )
Print( farrot_total, 0 )
#
Print( empty_line, 0 )



}