J. Z. 4/20/08

Here's a little tutorial about running the radiative transfer programs.
These programs assume a spherical geometry with arbitrary density,
temperature, and molecular abundance profiles vs radius.
The model includes emission/absorption by both dust and
the molecular lines. 

------------------------------------------------------------------------

I've generated some input files for the radiative transfer
problem based on the basic physical parameters of the
two-component model proposed by Weiss et al 2007 to explain
CO, HCN, and dust in APM 08279 (see their Table 3). The
"cold" component has a density of 1e5 and a temperature
of 65 K, while the "warm" component has a density of 1e4
and a temperature around 220 K. These parameters have been
taken to be constant vs radius. I've arbitrarily set the
physical size (diameter) of both components to 100 pc, which
is within the range given by Weiss et al. (Table 3).

Also, I have arbitrarily set the H20 abundance relative to
H2 to be 1E-6. This is much higher than is found for cool
galactic clouds but probably on the low side for warm
gas, above the ~ 90 K sublimation temperature for H2O
ice mantles to evaporate from dust grains.


There are three basic steps for performing the calculation

1. Generate a radial grid using the "logrid3.f" program
2. Solve the radiative transfer problem using "molnd_7.f".
3. Predict the emergent spectra using "observ_spec_v2.f".

This last step involves choosing a beam size for the 
synthetic observation. I've chosen to match the FHWM
beam size to the source size at a rest frequency of 3 THz.
The beamsize is assumed to scale with wavelength, as
you would expect for a diffraction limited telescope.
Therefore, at lower frequencies observed by Zspec
(1-1.5 THz rest frequency for this source), the beam is 
bigger than the source.



----------------------------------------------------------------------------

***
 weiss_hot_component.param
***

This file contains the parameters defining a radial grid
for simulating the hot gas component of APM08279 according
to the results of Weiss et al 2007
This file is used as input to the program logrid3.f, which
generates the output file:

***
weiss_hot_component.logrid
***

which is how you define the radial shells for the radiative transfer
program molnd_7.f

Here's an example run:

orthanc% ../bin/molnd_7
  Start new (N) or continue old (O) calculation ?
N
  Starting NEW calculation...
  Input file name for radial grid ?
weiss_hot_component.logrid
  Molecular transition data file ?
p-h2o.spec
  Collisional excitation rates file ?
p-h2o.collrates
  Dust parameters file ?
dust.param
  Molecular abundance factor (normally 1.0) ?
1.0
  Dust absorption factor (normally 1.0) ?
1.0
  Output file ?
weiss_hot_component.p-h2o.pop
  Use accelerated (1) or normal (0) lambda-iteration ?
1
  Relative fractional population change
  threshold for convergence (i.e. 1.E-3) ?
1.e-3
  Maximum number of iterations (i.e. 30) ?
10
  Solving rates for iteration # 0
  Solving transfer for iteration # 1
  Solving rates for iteration # 1
  Maximum fractional population change is  8.53544669
  This occurs for level 23( 7( 2, 6) ) in radial shell 50( / 50)
  Old population:   0.00529396742
  New population:   0.000554293065
  Population-weighted r.m.s pop. change is:  0.00939823663
  Solving transfer for iteration # 2
  Solving rates for iteration # 2
  Maximum fractional population change is  1.17859415
  This occurs for level 19( 7( 1, 7) ) in radial shell 50( / 50)
  Old population:   0.00187483721
  New population:   0.000860031051
  Population-weighted r.m.s pop. change is:  0.00477847345
....
this keeps going, converges after around 50 iterations, and
writes the output file:

***
weiss_hot_component.p-h2o.pop
***

Now you are ready to run the program that generates the
predicted spectra. You first need to create a little file

***
weiss_hot_component.observ
***

that specifies the FWHM beam size for the observation. Actually,
it specifies the FWHM beam size * frequency product, i.e.
it assumes that the telescope is diffraction limited.
You also need to specify the velocity range over which to
calculate the spectra. The calculations only need to be done
for positive velocities since the line is (almost) symmetric.
The slight asymmetry from the rising dust spectrum is ignored
(but this isn't hard to change if necessary).

OK, let's calculate the emergent spectrum:

orthanc% ../bin/observ_spec_v2
  Name of file created by molnd ?
weiss_hot_component.p-h2o.out
  File name giving observation info ?
weiss_hot_component.observ
  Abundance factor (normally 1.0) ?
1.
  Dust absorption factor (normally 1.0) ?
1.
  Output file ?
weiss_hot_component.p-h2o.txt
orthanc% 

which produces the output file

***
weiss_hot_component.p-h2o.txt
***

The .txt filename extension is chosen to allow the data
to be easily imported into MATLAB for plotting.

The above stuff needs to be repeated for the ortho-h2o lines;
these are treated separately since there are no radiative
transitions between the ortho and para species.

Then you need to repeat the calculations for for the "cold" component.

--------------------------------------------------------------------

The MATLAB script weiss_both_plot.m does the job of taking the
results in the .txt output files, combining ortho and para results,
and plotting the spectra for both cool and warm components.

--------------------------------------------------------------------
