LOC

LOC - Line radiative transfer with OpenCL

LOC is a program, or rather a family of programs, for the radiative transfer modelling of molecular lines. The main novelty is in the use of OpenCL libraries for the parallelisation of the computations so that they can be run on multi-core CPUs and GPUs. The basic concept is the same as in the case of SOC, which is a corresponding radiative transfer code for calculations of dust emission and scattering.

At this moment there exist versions of LOC for 1D spherically symmetric grids and for hierarchical 3D octree grids (see the page on GitHub ). In addition to “normal” molecular line modelling, one can model spectra with hyperfine structure by assuming either local thermodynamic equilibrium between the components or (in the 1D program) by doing the full calculation, including the effects of line overlap.

Topics:

Input files ^

Model cloud description ^

Cartesian 3D clouds

The 3D model is described by a separate binary file that defines (in the case of the Cartesian grid) the model dimensions and specifies for each cell seven quantities: the number density, the kinetic temperature, the amount of microturbulence (sigma), three components of the macroscopic velocity (vx, vy, and vz), and the fractional abundance. The file starts with the dimensions NX, NY, and NZ (three 4 byte integers), following by the data for the NX×NY×NZ cells, each consisting of the aforementioned seven numbers (4 byte floats; the data for a single cell being consecutive elements in the file).

Of the three spatial coordinates, it is the X coordinate that should change fastest in the file. This corresponds to the “C-order” of elements. One must remember that in Python arrays the last index runs fastest. In other words, if one creates the cloud as an python array before writing it to disk, that array must be indexed as CLOUD[k, j, i, ifield], where k corresponds to the LOC Z axis and i to the LOC X axis. See Full example of a LOC run or the file test_loc_velocity_axis.py in GitHub for practical examples that hopefully clarify the coordinate system (LOC vs. Python).

The values in the file can be in physical units. For example, the density would be in units cm-3 and the microturbulence and velocities in km/s. Usually the density would be the density of H or of H2, the abundance being the abundance of the examined species relative to that density. However, the raw values read from the file may be rescaled using appropriate keywords in the parameter file.

Octree clouds

The octree cloud file starts with five integer (int32) values: NX, NY, NZ, OTL, CELLS. Here (NX, NY, NZ) are the dimensions of the root grid, OTL is the number of hierarchy levels in the octree structure, and CELLS is the total number of cells in the model. This is followed by the data for the root grid, followed by data for the first refinement level and so forth. On each hierarchy level, the file contains data for the same seven quantities as for the regular 3D clouds (density, Tkin, sigma, vx, vy, vz, abundance). Each parameter vector (containing one float32 number per each cell on that hierarchy level) is preceded by the number of cells on the hierarchy level. LOC follows the c-convention where the fastest changing index is called x. If one is working with Python arrays, there the last array index corresponds to the x axis and the first index to the z axis: array[z, y, x]. The overall file structure is

  NX, NY, NZ, OTL, CELLS
  cells_0, rho_0, cells_0, Tkin_0, cells_0, sigma_0, ...
  cells_1, rho_1, cells_1, Tkin_1, ...
  cells_N, rho_N, ..., cells_N, vz_n, cells_N, abundance_N

In the above cells_* are the number of cells per hierarchy level (int32 scalars) and the others (rho_, Tkin_, …) are vectors with those dimensions (e.g. rho_1[cells_1], with float32 values).

The actual grid hierarchy is contained in the rho_* vectors. When a cell is refined (divided into eight subcells) the density value rho_x[i] is replaced by value -(float )&ind, when ind is the index of the first subcell in the vector rho_(x+1), i.e. the density vector of the next hierarchy level. In other words, one takes the subcell index ind, casts that as a float32 number, and stores it multiplied by -1 (to enable distinction between actual density values and rho_ entries that are links to sub-cells). For the other quantities (Tkin, sigma, etc.) the values in the parent cells (i.e. refined cells) are usually undefined and are not used by LOC.

Initialisation file ^

The following list will describe briefly the keywords used in the initialisation file (ini file). This is the file that is given as a command line parameter for the LOC programme. Each line in the parameter file starts with a keyword, possibly followed by one or more parameters separated by spaces. Comment lines in inline comments are indicated with a ‘#’ character.

  • abundance float scaling of the abundance values read from the model file

  • alpha float for 1D models only, parameter affecting the radial discretisation: the default linear discretisation ri, where i is the shell index. Radii will be placed at ri**(1/(1+alpha)). Values alpha<0 will thus make the 1D grid more dense towards the centre of the model, and with values alpha>0 the grid is more dense towards the outer parts. The value of alpha must be greater than -1.

  • angle float apparent size of (root) grid cell in units of arcsec. For 1D models, this is corresponds to the radius of the spherical model.

  • bandwidth float total bandwidth used in the simulations and for the output spectra, in units if km/s

  • cabfile filename optional file giving the fractional abundance of each collisional partner (abundances relative to the density values used to define the model). In 3D models (LOC_OT.py), the file starts with NX, NY, NZ, and the number of collisional partners (four 4 byte integers). The rest of the file lists the abundances for every model cell, the values of a given cell being consecutive (float32 numbers, with the index over partners running faster than the index over cells: if created in Python, the array would have the shape CAB[CELLS, PARTNERS]). In 1D models (LOC1D.py) the file content is just the CELLS*PARTNERS abundance values (float32 numbers), without any leading integers.

  • channels int number of velocity channels (in the simulation and the output spectra)

  • cloud filename name of the model file (with the density etc. values for every model cell)

  • colden 1 (with any integer argument >0) instructs LOC to save information of the column density. In the case of 1D models (LOC1D.py) column densities are saved to a binary file, as described in section on column densities . In the case of 3D models (LOC_OT.py) results are saved to FITS files, for the same pixels as for which the spectra are saved. The file name will be prefix.NH2.000.fits. If maps are saved towards multiple directions, ‘000’ is replaced with a running index over the directions. LOC_OT.py saves at the same time also some other files, as described here .

  • cooling if given, the total cooling rate of the current species is calculated and saved to a file (one value per cell)

  • cores int asks LOC to only use the specified number of cores. It depends on the OpenCL implementation, whether this works (has been tested with Intel CPU driver, should in principle also work with POCL and with AMD GPU drivers).

  • crtemit specifies the name for a file with additional continuum emission in each cell. This binary file should start with two int32 numbers that are the number of cells and the number of transitions (matching the current run). This is followed by the main array emission[cells, transitions] of float32 numbers, the continuum emission in units of photons / s / Hz / H2 (assuming that the model is defined in terms of H2 density). In the file, the transition index runs faster.

  • crttau specifies the name for a file with additional continuum opacity in each cell. This binary file should start with two int32 numbers that are the number of cells and the number of transitions (matching the current run). This is followed by the main array tau[cells, transitions] of float32 numbers, the optical depth per pc. In the file, the transition index runs faster.

  • density float factor by which the density values read from the model file are multiplied before calculations

  • device string select the OpenCL device for the calculations. The string can be ‘c’ or ‘g’, to select the first available CPU or GPU, respectively. Otherwise, the string can be any unique string that appears as part of the device name. One can check the available OpenCL devices and their names using the clinfo utility, which is available in most repositories or can be obtained from github.

  • direction float float [string] specifies the direction of the observer as angles (in degrees) from the positive Z axis and rotation from positive X towards Y. The optional third parameter is an additional string included in the name of the output files. For direction 90 0, the observer is in the direction of the positive X axis and the maps are oriented so that cloud Y axis is right and the Z axis up. For direction 0 0, observer is in the direction of positive Z axis and in the maps Y is right and X is down. For direction 90 90, the observer is in the direction of the positive Y axis and in the maps the cloud X axis points left and the cloud Z axis up.

  • distance float specifies the distance between the model and the observer, in units of pc

  • fits [RA DEC] makes LOC_OT.py to write results into FITS files instead of custom binary files. This includes the files for the spectra, the optical depths, and the column densities. The optional parameters are the centre coordinates in degrees that are used in the written FITS files. If there is only one argument, that is interpreted as a flag whether to use FITS files (‘fits 1’) or not (‘fits 0’)

  • fraction float scaling applied to the abundance values read from the model file (synonym for keyword abundance)

  • gpu int for any value of int>0, LOC tries to select for the calculations a GPU instead of a CPU device. See also keyword platform. The keyword device is an alternative way to select the exact platform and device to use.

  • grid float specifies the step between spectra in the output maps (~pixel size). The unit is arcsec.

  • hfsfile filename gives the name of the file that specifies the relative velocities and intensities of hyperfine components, when these are assumed to be in LTE.

  • isotropic float specifies the temperature (in degrees of Kelvin) for isotropic background radiation, assumed to have the spectrum of a blackbody.

  • iterations int specifies the number of iterations (consisting of the simulation of the radiation field and the updating of the level populations)

  • load filename at the beginning of the run, load previously saved level populations from the given binary file.

  • local int will change the default local work group size with the size given as the argument (see OpenCL documentation: defaults are probably 1 for CPU and 32 for GPU)

  • levels int specifies the number of lowest energy levels included in the computations (less or equal to the number of levels for which the molecule file contains data for)

  • lowmem turns on some optimisations to minimise the memory usage, at the expense of some increase in the run time

  • mapview theta phi nx ny [ cx, cy, cz] is in LOC_OT.py an alternative way to specify the mapping parameters (instead of using keywords points and direction). Theta is the angle from the positive z-axis, and phi is the rotation around the z axis from the positive x axis, both given in degrees. Parameters nx and ny are the number of pixels in the output map. The optional parameters cx, cy, and ‘cz give a position in cell units that will be towards the centre of the map.

  • molecule filename specifies the name for the file containing the molecular data. The file can be in lamda format.

  • nray int applies only to 1D models (LOC1D.py), where it is the number of rays (number of impact parameters) used during the calculation of the radiation field. Depending on the model (number of shells) values are typically some hundreds or thousands.

  • nside int specifies the angular grid for the simulated rays, which follows the healpix scheme. The total number of rays per model surface element is 12 times this value squared.

  • octree filaname is in LOC_OT.py similar to the cloud keyword and specifies the name of the cloud file. However, this was used mainly to test different ray-tracing implementations, which are indicated by an optional integer attached to the keyword. One should probably not use this keyword, although for example octree4 should be safe (cloud would correspond to octree0).

  • overlap filename instructs LOC to take into account the possible overlap line profiles of some transitions (not assumed to be in LTE). The argument is the name of a file that lists these transitions (see molecular parameters )

  • platform int (int) specifies the OpenCL platform to be used in the calculations. The optional second integer parameter selects a device within the platform (running index over all devices within the platform that are of the requested type, CPU or GPU). By default all available platforms are tried in order, choosing the first platform and the first device on the platform. Instead of the gpu and platform keywords, one can use just the device keyword.

  • points int int specifies the dimensions of the output spectrum files (the number of spectra along horizontal and vertical directions)

  • prefix string prefix for the output files

  • tausave int for int>0, LOC1D.py saves optical depths to a file with suffix “.tau”, for the same lines of sight as for the save spectra. In LOC_OT.py, optical depths can be saved to FITS files (see section on optical depths )

  • save filename gives the name for a binary file where the level populations are saved at the end of a run

  • sigma float scaling to be applied to the microturbulence value read from the file for the model description

  • spectra int int … specifies the transitions for which the spectra are to be computed and saved. The arguments are pairs “upper lower”, two integers that specify the upper and lower energy level of the transition. The energy levels are numbered sequentially started from zero (in the order these appear in the molecule file). For example, in the case of linear molecules, the level index is likely to be the same as the rotational quantum number.

  • speray int, for 1D models, this specifies that spectra will be written for the number of offsets given here as the argument, instead of the default number of (equidistant) rays given by the keywork nray.

  • stop float iterations are automatically stopped when the largest relative change in level populations fall below this limit. Only levels up to the level specified with the keyword uppermost are tested.

  • temperature float scalar scaling applied to the kinetic temperature values read from the model description file.

  • transitions int int … list of transitions (using running index of the transitions) for which the excitation temperatures are to be saved

  • uppermost int specifies the uppermost energy level that is checked for level population changes (see keyword stop)

  • velocity float scaling applied to the velocity values (three components if macroscopic velocity) that are read from the model description file

Molecular parameters ^

LOC reads directly files in the format of the Lamda database .

If the abundances of the collisional partners are assumed to be constant over the model volume, these abundances (with default value 1.0) can be added to the molecule file as a floating point number after the number of temperatures (same line, after the number of temperatures). If the abundances change from cell to cell, one needs to use the cabfile keywords of the ini files.

Possible hyperfine structure of the lines can be handled assuming LTE conditions between the hyperfine components (instead of the general line overlap in frequency space; see the overlap keyword). In this case, the molecule file should specify only the main energy levels, not the individual hyperfine components. Calculations require the use of the keyword hfsfile and a corresponding text file that lists the relative velocities and intensities of the components. For each hyperfine group, the file lists on one line the indices of the upper and lower level (running indices starting with zero) and the number of hyperfine components. The following lines then list the relative velocity (km/s) and the relative intensity of each component. Below is a fictious example of the file format.

 1 0 5       # first transition (=upper and lower level) and the number of components
 -2.4  0.05  # velocity (km/s), relative intensity (sum equal to one)
 -1.2  0.25  # ...
 0.0   0.40
 1.54  0.20
 2.92  0.10

 3 2 3       # a second transition with hyperfine structure, this with three components
 -1.1  0.3
 0.0   0.5
 2.2   0.2

On the other hand, if the molecule file already describes transitions for which the profiles may overlap, one can use the keyword overlap. Its argument is also a text file, which this time lists the transitions for which the frequency overlap should be taken into account. Below is an example of the structure of that text file.

  3     # number transitions with overlapping line profiles
  5 4   # each transition specified using the upper and lower level index (>=0)
  6 5 
  9 8

Output files ^

Spectrum files ^

The spectra are saved to binary files. Currently the file names are similar to prefix.molecule.01-00.bin, where prefix is given in the ini-file, molecule is replaced by the actual name of the molecule (species), and the number give the transition (based on the 0-offset numbering of energy levels). The file contains:

  • three 4-byte integers that give the number of pixels along the two map pixels (NX, NY), and the number of velocity channels (NCHN)
  • two 4-byte floating point numbers giving the velocity of the first channel and the channel width (km/s)
  • a cube with dimensions (NX, NY, NCHN+2), consisting of 4-byte floating point numbers. Each vector of NCHN+2 elements starts with the x- and y-offsets, followed by the channel values

Excitation temperature ^

Excitation temperatures are saved to binary files, named combining the prefix (as given in the ini-file), the name of the molecule and the transition, and with an ending of tex. The file

  • starts with four 4-byte integers, the dimensions of the model cloud (NX, NY, NZ), and the number of excitation levels included in the calculation (LEVELS)
  • ends with a (NX, NY, NZ) cube of excitation temperature values for the chosen transition, saved as 4-byte floating point numbers

Optical depths ^

In 1D case, the optical depths (saved using the keyword “tausave 1”) is a binary file. The beginning of the file is: number-of-offsets (int32), number-of-channels (int3), velocity-of-first-channel (float32, km/s), channel-width (float32, km/s). These are followed by number-of-offsets times number-of-channels values (float32) giving the optical depths. The optical depths are saved for all the same transitions for which the normal spectra are saved.

In the 3D case, optical depths are saved by default to a binary file with the suffice ‘.tau’, the file containing just the peak optical depth of each spectrum. If one is saving data to FITS files (using ‘fits 1’ in the ini file), the resulting 3D FITS file (…tau#.fits, with # standing for the view direction, the rest of the name including the name of the molecule and the transition) will contain the optical depths for all spectra and all channels.

Column densities ^

In the 1D case, when spectra are written, the program also saves a file (suffix .colden) that lists, among other things, column densities for the radially requidistant set of rays. The contents of this binary file are the following: num (int32, number of radial points), offset[num] (array of float32 values, the offset in cm for each radial point), N[num] (array of float32 values, total column density for each offset, un inuts of cm-2), Nmol[num] (array of float32 values, column density of the simulated species, in units of cm-2), and tau[num] (array of float32 values, peak optical deptha over the line profiles, one value per offset).

In the 3D case, column density maps (if the ini file includes a line ‘colden 1’) are saved to FITS files, where the filename includes the string ‘NH2’. LOC_OT.py saves at the same time also similar maps for the column density of the examined species (file ‘…N…fits’, taking into account the abundance variations), the mass-weighted line-of-sight kinetic temperature (file ‘…Tkin…fits’), and the mass-weighted line-of-sight mean velocity (file ‘…VLOS…fits’).

Cooling rates ^

The cooling rates are saved to a plain binary file that consists of one 4-byte value per cell. The values are the cooling rates in cgs units.

Examples of run times for a 3D model ^

Below are some example of the run times on different computer systems. The model has a 121³ 3d cartesian grid. One is calculating CO spectra, taking into account the ten first energy levels and running the simulations with 128 velocity channels. The OpenCL platforms include Intel, NVidia, and AMD GPUs as well as CPUs that are used either with the open-source POCL or the Intel’s own OpenCL implementations. The run times are given in seconds, for a complete run that includes five iterations and the final writing of the spectrum files. Therefore, the serial Python host code and the IO stand for several seconds of the total run time.

  Run time    Device
  (seconds)
  =====================
  
  Dell XPS laptop
  ------------------------------
  853.87    CPU: Intel i7-8550U 1.80GHz, 4 cores, POCL implementation
  855.89    CPU: Intel i7-8550U 1.80GHz, 4 coresm Intel driver
  279.62    GPU: Intel UHD Graphics 620, Intel driver
  249.54    GPU: Intel UHD Graphics 620, Intel NEO driver
  46.39     GPU: AMD Radeon VII (external GPU enclosure)
        
  Laptop 2 (Tuxedo)
  ------------------------------
  351.70    CPU: Intel i7-8700K CPU @ 3.70GHz, 6 cores, POCL
  31.24     GPU: NVidia GeForce GTX 1080

  Desktop computer
  ------------------------------
  64.68     CPU: AMD Ryzen Threadripper 3960x, 24 cores, POCL
  30.19     GPU: NVidia GeForce GTX 1080 Ti
  14.89     GPU: NVidia GeForce RTX 3090

GPUs are generally more efficient that CPUs but also show a factor of ten differences between the cards. A modern multi-core CPU system does get close. However, the comparison does not include the latest GPU cards, which can be expected to be again a few times faster than the fastest GPUs in the above list.

Full example of a LOC run ^

In the following example we create a 3D model for a cloud, model the CO emission with LOC, and plot some spectra and a cross section of the excitation temperature cube. We use Python to make the model and to plot the results. The scripts can be found in the file LOC_Cartesian_test.zip in our GitHub page .

Cloud model

We discretise the model onto a regular Cartesian grid of 64³ cells. Therefore, we can use the simpler format where the cloud file starts with the dimensions of the root grid (integers) followed by the four dimensional cube, with three spatial coordinates and the fourth index running over the parameters (density, kinetic temperature, three components of the velocity, and the molecular aundance. The following python script make_cloud.py writes the cloud file. Note that when we talk about coordinates X, Y, and Z in LOC, these refer to the “c-order” where it is the “x” coordinate that runs fastest (in the file and in the c-program arrays). In contrast, in python arrays it is the last index that runs fastest. This means that also in the following script a Python array element C[k, j, i, ifield] corresponds to index k along the LOC Z axis, index j along the LOC Y axis, and index i along the LOC X axis!

from numpy import *
from numpy.random import randn

# allocate the cube containing seven parameters for each cell
NX, NY, NZ = 64, 64, 64
# we have same number of cells along each dimension
# however, note the order of dimensions in the Python array !!
C       = ones((NZ, NY, NX, 7), float32)  # start with values 1.0 for all parameters

# calculate cell distances from the centre
k, j, i = indices((NZ, NY, NX), float32)
c       = 0.5*(NX-1.0)                # coordinate value at the centre (since NX=NY=NZ)
i, j, k = (i-c)/c, (j-c)/c, (k-c)/c   # coordinates in the [-1,1] range
r       = sqrt(i*i+j*j+k*k)           # distance from the centre

# Set the values in the cube
# C[k, j, i, ifield] is an element where the LOC x-coordinate is i, y-coordinate j,
#                    and z-coordinate is k !!
C[:,:,:,0]   = exp(-3.0**r*r)         # density values, 3D Gaussian
C[:,:,:,1]  *= 15.0                   # Tkin=15 K in every cell
C[:,:,:,2]  *=  1.0                   # leave microturbulence at 1 km/s
C[:,:,:,3]   = randn(NZ,NY,NX)        # vx, macroscopic velocity
C[:,:,:,4]   = randn(NZ,NY,NX)        # vy  
C[:,:,:,5]   = randn(NZ,NY,NX)        # vz
C[:,:,:,6]   = 1.0                    # abundance values will be scaled later

# Write the actual file
fp = open('my.cloud', 'wb')
asarray([NX, NY, NZ], int32).tofile(fp)  # note again the order of dimensions !!
C.tofile(fp)  # in the file the x-coordinate changes fastest (C-order)
fp.close()

Ini file and LOC run

We define the run-time parameters with a text file my.ini. As described above, this consists of keywords, followed by one or more parameters. The commented file looks like this:

cloud          my.cloud        #  cloud defined on a Cartesian grid
distance       100.0           #  cloud distance [pc]
angle          3.0             #  model cell size  [arcsec]
molecule       co.dat          #  name of the Lamda file
density        1.0e4           #  scaling of densities in the cloud file
temperature    1.0             #  scaling of Tkin values in the cloud file
fraction       1.0e-4          #  scaling o the fraction abundances
velocity       0.5             #  scaling of the velocity
sigma          0.5             #  scaling o the microturbulence
isotropic      2.73            #  Tbg for the background radiation field
levels         10              #  number of energy levels in the calculations
uppermost      3               #  uppermost level tracked for convergence
iterations     5               #  number of iterations
nside          2               #  NSIDE parameter determining the number of rays (angles)
direction      90 0.001        #  theta, phi [deg] for the direction towards the observer
points         64 64           #  number of pixels in the output maps
grid           3.0             #  map pixel size [arcsec]
spectra        1 0  2 1        #  spectra written for transitions == upper lower (pair of integers)
transitions    1 0  2 1        #  Tex saved for transitions (upper and lower levels of each transition)
bandwidth      6.0             #  bandwidth in the calculations [km/s]
channels       128             #  number of velocity channels
prefix         test            #  prefix for output files
load           test.save       #  load level populations
save           test.save       #  save level populations
stop          -1.0             #  stop if level-populations change below a threshold
gpu            1               #  gpu>0 => use the first available gpu

The keyword must be the first string on a line and its arguments are separated by at least by one space (or tab). One can add comments after the arguments or as separate lines starting with ‘#’.

Above the keywords density, temperature, fraction, velocity, and sigma all apply a scaling to the raw values read from the cloud file. We also explicitly choose the number of map pixels (points) and the map pixel size (grid) so that they correspond to the projected size and discretisation of the model cloud. The output files start with the prefix test. We do not specify the OpenCL device so LOC tries to use the first available CPU device. The molecular data are in the file co.dat, which can be downloaded directly from the Leiden Lamda database .

The actual radiative transfer calculations are done from the command line:

LOC_OT.py   my.ini

This should produce a few files, including the spectra in the files test_CO_01-00.spe and test_CO_02-01.spe and the excitation temperatures in test_CO_02-01.tex and test_CO_02-01.tex. The estimated level populations are saved to test.save and can be again loaded (keyword load) as the initial values if the calculations need to be continued.

Plotting spectra

Both the spectra and the excitatio temperatures are saved to binary files. The following Python routines can be used to read the files. These can be used also in the case of octree models (instead of the regular Cartesian grid of this example) - except that the reshape command in LOC_read_Tex_3D must be omitted (there the excitation temperatures are stored as a single, one-dimensional vector that covers all hierarchy levels).

import os, sys
from matplotlib.pylab import *

def LOC_read_spectra_3D(filename):
    """
    Read spectra written by LOC.py (LOC_OT.py; 3D models)
    Usage:
        V, S = LOC_read_spectra_3D(filename)
    Input:
        filename = name of the spectrum file
    Return:
        V  = vector of velocity values, one per channel
        S  = spectra as a cube S[NRA, NDE, NCHN] for NRA
             times NDE points on the sky,
             NCHN is the number of velocity channels
    """
    fp              =  open(filename, 'rb')
    NRA, NDE, NCHN  =  fromfile(fp, int32, 3)
    V0, DV          =  fromfile(fp, float32, 2)
    SPE             =  fromfile(fp, float32).reshape(NRA, NDE,2+NCHN)
    OFF             =  SPE[:,:,0:2].copy()
    SPE             =  SPE[:,:,2:]
    fp.close()
    return V0+arange(NCHN)*DV, SPE
    
    
def LOC_read_Tex_3D(filename):
    """
    Read excitation temperatures written by LOC.py (LOC_OT.py) -- in case of
    a Cartesian grid (no hierarchical discretisation).
    Usage:
        TEX = LOC_read_Tex_3D(filename)
    Input:
        filename = name of the Tex file written by LOC1D.py
    Output:
        TEX = Vector of Tex values [K], one per cell.
    Note:
        In case of octree grids, the returned vector must be
        compared to hierarchy information (e.g. from the density file)
        to know the locations of the cells.
        See the routine OT_GetCoordinatesAllV()
    """
    fp    =  open(filename, 'rb')
    NX, NY, NZ, dummy  =  fromfile(fp, int32, 4)
    TEX                =  fromfile(fp, float32).reshape(NZ, NY, NX)
    fp.close()
    return TEX

We will plot the J=1-0 and J=2-1 spectra towards the cloud centre, cross sections of the J=1-0 and J=2-1 excitation temperature distributions, and a map of the line area W for the J=1-0 transition.

V, T10 = LOC_read_spectra_3D('test_CO_01-00.spe')
V, T21 = LOC_read_spectra_3D('test_CO_02-01.spe')
Tex10  = LOC_read_Tex_3D('test_CO_01-00.tex')
Tex21  = LOC_read_Tex_3D('test_CO_02-01.tex')
NY, NX, NCHN = T10.shape  # number of pixels and channels

# One spectrum from (about) the model centre
subplot(221)
plot(V, T10[NY//2, NX//2, :], 'b-', label='T(1-0)')
plot(V, T21[NY//2, NX//2, :], 'r-', label='T(2-1)')
legend(loc='lower center')
title("Centre spectra (K)")
xlabel(r'$v \/ \/ \rm (km s^{-1})$')
ylabel(r'$T \rm _A \/ \/ (K)$')

# Cross section of the Tex(1-0) cube
subplot(222)
imshow(Tex10[:,:, Tex10.shape[2]//2])
title("J=1-0 excitation temperature (K)")
colorbar()

# Cross section of the Tex(2-1) cube
subplot(223)
imshow(Tex21[:,:, Tex21.shape[2]//2])
title("J=2-1 excitation temperature (K)")
colorbar()

# J=1-0 line area
subplot(224)
W     =  (V[1]-V[0])*sum(T10, axis=2)
imshow(W)
title("Line Area W(1-0) (K km/s)")
colorbar()

show()

The script produces the following figure. The grainy appearance of the Tex and W maps is not caused by numerical noise but by the random velocity field. This can be easily confirmed by recalculating the model without the random velocities (e.g. setting “velocity 0.0” in the ini file or recreating the cloud with some systematic velocity field). The above python programmes (cloud creation and plotting) can be found in LOC_Cartesian_test.zip in the LOC directory of our GitHub page .