#
#  cloud_rad_field.arts
#
# This file demonstrates how a 1D scattering calculation for one specified 
# frequency can be performed. The radiation field inside the cloudbox is
# calculated and all iterations are written into files which can be plotted
# using the MATLAB plotting scripts (set_plot.m, plot_rad.m).

Main
{
##############################
#     General settings       #
##############################

# The frequency grid has to contain at least two frequencies.
# For a monochromatic calculation you have to set the first frequency value
# to the value you want to calculate. For the second frequency, the
# scattering calculations will not be performed.  

  VectorNLinSpace( f_grid ) {
      start = 3.18e11
      stop  = 3.19015e11
      n     = 2	
  }

# You do the scalar calculation ( stokes_dim = 1 ) or teh vector 
# calculation for a arbitrary number of components ( stokes_dim = 2, 3 or 4).

  IndexSet( stokes_dim ) {
      value = 2
  }

AntennaSet1D{}

output_file_formatSetAscii{}

##############################
#    Atmosphere              #
##############################
	
# Set the atmospheric dimension.

  AtmosphereSet1D {
  }

# Definition of the pressure grid.

  VectorNLogSpace( p_grid ) {
     start = 1000e2
     stop  = 100e2
     n     = 100
  }

WriteXML(p_grid){""}
	

# Latitude for a 1D observation.
	
  NumericSet( lat_1d ) {
     value = 45
  }
  
#The meridian angle for a 1D observation.

  NumericSet( meridian_angle_1d ) {
    value = 0
  }

# Define a set of gas species for the atmosphere.
 
 gas_speciesSet{
	species = [ "H2O",
	 	"O3",
		"O2", 
		"N2" ]
	 }		
	
# Read atmospheric profiles (temperature, pressure, altitude, vmrs).
# You have to specify the path for the scenario you want to have. The 
# data has to be in XML-format.

  AtmRawRead{
	basename = "@ac_arts_xml_data@/atmosphere/fascod/tropical"
	}

# Create workspace variables, which are used in the program. 
 
 AtmFieldsCalc{}

# Write the created workspace variables. 
 
 WriteXML(t_field){""}
 WriteXML(vmr_field){""}

####################################
#  ground setup                    #
####################################

   r_geoidWGS84 {
  }

  MatrixSet( z_ground ) {
     nrows = 1
     ncols = 1
     value = 0.5e3
  }

  AgendaSet( ground_refl_agenda ) {
    Ignore( a_los ) {
    }
    Ignore( r_geoid ) {
    }
    Ignore( z_ground ) {
    }	
    GroundTreatAsBlackbody {
    }
  }


################################
#    Gaseous absorption        #
################################

# Initialize the lookup table workspace variable.

 gas_abs_lookupInit{}

# Read a pre-calculated lookup table from a file.
	
ReadXML(gas_abs_lookup){"/home/home02/sbuehler/arts_calc/abs_lookup/lookup_01_table.xml"}


# Adapt the lookup table. In the varible *gas_abs_lookup* holds only the 
# specified frequencies and species.

gas_abs_lookupAdapt{}

# Write the lookup table.

WriteXML(gas_abs_lookup){""}

# Agenda for scalar gas absorption calculation.

AgendaSet( scalar_gas_absorption_agenda ){
	abs_scalar_gasExtractFromLookup{}	
	}

# Agenda for calculating optical properties of gaseous species
# (extinction matrix and absorption vector)

  AgendaSet( opt_prop_gas_agenda ) {
	ext_matInit{}
	abs_vecInit{}
      ext_matAddGas{}		 
      abs_vecAddGas{}
   }


###################################
#     Scattering                  #
###################################

# Angular grids for scattering calculation.

  VectorNLinSpace( scat_za_grid ) {
     start = 0
     stop  = 180	
     n     = 30
  } 

 VectorNLinSpace( scat_aa_grid ) {
    start = 0
   stop  = 360
  n     = 30
}
 
# Define the scattering region (cloudbox)
 
  CloudboxSetManually {
      p1   = 335e2
      p2   = 265e2
      lat1 = 0
      lat2 = 0
      lon1 = 0
      lon2 = 0
  }

  ArrayOfIndexPrint( cloudbox_limits ) {
   }

# Specification of particle properties.

ParticleTypeInit{}

# Read amplitude matrix from database

#ParticleTypeAdd{
#  filename_amp_mat = "@ac_arts_xml_data@/scattering/cyl_s200_f318-324_ar0.5_o0.0_ampmat.xml"
##filename_amp_mat = "@ac_arts_xml_data@/scattering/sph_s200_f318-324_ampmat.xml" 
#filename_pnd_field = "  "
#}

#WriteXML( amp_mat_raw ){""}
#WriteXML( pnd_field_raw ){""}
    
ReadXML( amp_mat_raw ){""}
#ReadXML( pnd_field_raw ){""}   
 

# Preliminary pnd_field (instead of data_file)

  Tensor4Set (pnd_field){
	nbooks = 1
	npages = 100
	nrows = 1
	ncols = 1
	value = 1297
  }

# So far we always calculate only for 1 particle type. 

  VectorSet (part_types){
	length = 1
	value = 1
  }

# Calculate single particle optical properties
 
 AgendaSet(spt_calc_agenda){
      pha_mat_sptCalc{}
      ext_mat_sptCalc{}
      abs_vec_sptCalc{}
  }

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

 AgendaSet( opt_prop_part_agenda ){
	ext_matAddPart{}
	abs_vecAddPart{}
  }



#####################################
#    Radiative Transfer             #
#####################################


#  Cosmic background 

  AgendaSet( i_space_agenda ) {
    Ignore( a_pos ) {
    }
    Ignore( a_los ) {
    }
    MatrixCBR( i_space, f_grid ) {
    }
 #   MatrixPrint( sensor_pos ) {
 #   }
  }


#  Standard radiative transfer function.
AgendaSet( rte_agenda ){
	RteEmissionStd{}
}
 	
# Calculation of propagation path

  AgendaSet( ppath_step_agenda ) { 
	ppath_stepGeometric{
	       lmax = -1
		}
#	PpathPrint(ppath_step){}
	 } 


# Define the convergence test method. This agenda is executed in 
# *i_fieldIterate* which is part of scat_mono_agenda.
  AgendaSet( convergence_test_agenda) {
	iteration_counterIncrease{}
	Tensor6WriteIteration( i_field ) {
	     iterations = [0]
  	}
	convergence_flagAbs{
#	epsilon = [1e-17, 1e-19, 1e-21, 1e-21] 
	epsilon = [1e-17, 1e-19 ]
#	epsilon = [1e-17]
	}
  }


#-----------------------------------------------------------------------------
#  scat_mono_agenda 
#----------------------------------------------------------------------------- 
#
# This agenda is executed in the method *ScatteringMain*. It does the whole 
# scattering calculation for 1 frequency.
#
# In this example it includes the following methods:
#   1. amp_matCalc:
#      Here the amplitude matrix data is extracted for one frequency 
# 	(specified in ScatteringMain) and interpolated on the 
#	scattering angle grids. 
#   2. i_fieldSetClearsky: 	
#      Set the boundary condition. The clearsky field is interpolated on all 
#      points inside the cloudbox.
#   3. i_fieldIterate:
#      This is the main function which performs the iterative solution method.
#      It uses the following methods and agendas:
#	  1. scat_fieldCalc: Calculate the scattered field unsing 
#	     	trapzoidal integration method
#	  2. i_fieldUpdate1D: Performs a radiative transfer step and 
# 		updates the radiation field.
#		The following agendas are required:
# 			scalar_gas_absorption_agenda
#				(Calculate scalar gas absorption)
#			spt_calc_agenda
#				(Calculate single particle properties)
#			opt_prop_part_agenda
#				(Calculate particle optical properties)
#			opt_prop_gas_agenda
#				(Calculate gaseous optical properties)
# 			ppath_agenda 
#				(Calculate geometrical propagation path)
#	  3. Execute convergence_test_agenda.
#   4. scat_i_put:
# 	Put the obtained solution in boundary variables (for further 
#       use in the clearsky part)

AgendaSet(scat_mono_agenda){
  amp_matCalc{}
  i_fieldSetClearsky{}
  output_file_formatSetAscii{}
  WriteXML(i_field){"initial_field.xml"}	
  i_fieldIterate{}
  scat_iPut{} 
 }		

# Initialize variables for the scattering calculation.
ScatteringInit{}

# Calculate radiation field on the cloudbox boundary.
CloudboxGetIncoming{}

# Start scattering calculation.
# The ScatteringMain function executes *scat_mono_agenda* for each frequency
# defined in f_grid.
ScatteringMain{}
	
}