J and Astronomy

J is an extremely terse programming language. It is interactive, i.e., you type in an expression and it is executed as soon as you hit ``return''. One strength is that whole arrays of numbers can be represented by a single name and can be manipulated without any explicit reference to the indices of individual elements of the array. J is descended from APL, and like it has a rich array of operations. Unlike APL, it uses only the normal ASCII character set -- operators are the usual symbols, plus others formed by appending a period or colon after ASCII symbols, i.e., + , +. and +: are distinct functions. Thus "B % A" where B and A are numbers or arrays represents B divided by A, while "%. A" is the inverse of the matrix A and "B %. A" computes the solution of the linear system A x = B. As another example, "a ^ b" is a to the power b, but "a ^:n" applies the operation "a" n times, while "a ^:_" iterates "a" until convergence, reducing the notational apparatus drastically. At present, J is *free*, and runs under Linux, OS X, and (even) Windows. Get it at http://www.jsoftware.com.

Some definitions referred to (and needed for) the foregoing: local_defs.ijs. Also note that where I have written a file name as "/your_path/..." you should replace this by
the full path name to the folder where you have put the .ijs script in question.
Check back later as this page will grow. (In the code given here, anything following "NB." is a comment.)


Astronomy

Here's a routine to correct coordinates for precession: precess.ijs.

This is for the integration of the Lane-Emden (polytrope) equation: polytrope.ijs.

This code evaluates scattering and absorption light by spheres using Mie theory: mie.ijs.

This Mie theory code evaluates scattering of polarized light by spheres as a function of the scattering angle: angle_scat.ijs.

This code gives the location in the sky of the sun, moon & planets: Sky.ijs. It needs some data files: SS_DAT


Multi-level atoms

The following code deals with the solution of an n-level atom or ion under conditions relevant to the interstellar medium.
The J routines are here: N_pop.ijs.
There must be a file of atomic data for each ion considered. Some are given in this directory: ATOMIC_DATA.
Just a note on the sort of output produced:N_pop.txt.
Here is a discussion of the theory behind these routines: N-level.pdf.


Transfer of radiation in plane-parallel atmospheres

Here are some routines for evaluating the radiation field in plane-parallel atmospheres, based on the integral equations of the problem. We assume the source function and/or Planck function can be well approximated by a cubic spline. A J routine "SpMD.ijs" (see below) is used to construct matrix representations of the integral of the function against against a kernel function, based on this spline representation.

For example the routine IQMsi.ijs creates a matrix which operates on the source function to give the emergent intensity from a semi-infinite, plane parallel atmosphere, I(mu), as a function of the emergent angle mu = cos[theta].

Here is another routine, MTsi.ijs, which creates the matrix representations of the integration of a source function against various exponential integrals. For example, the $\Phi$ transform evaluates the flux at each level in the atmosphere by integration of the source function against the 2nd exponential integral. Perhaps the most important transform is the $\Lambda$-transform, which gives the mean intensity at a given point in the atmosphere by integration over the source function. For a given grid of optical depths "tau", the expression "1 MTsi tau" gives $\Lambda$-transform and "2 MTsi tau" gives the $\Phi$-transform. Other integers on the left of MTsi give transforms involving higher order exponential integrals useful for applications involving polarization (see below). Note that MTsi assumes a semi-infinite atmosphere and thus integrates from the last (largest) tau point to infinity, using a linear extrapolation of the source function.


Polarization

If the polarization of the scattered radiation is taken into account, the equations become more complicated. If the scattering is by free electrons or molecules, the scattering will follow Rayleigh's law. The integral equations in this case can be found in Harrington (1970) Astrophysics & Space Science, 8, 227-242. In addition to the Lambda-transform, we introduce two additional transforms, the "M-transform" and the "N-transform", which involve higher exponential integrals. The J code to create the matrix representation is "3 MTsi tau" for the "M-transform" and "5 MTsi tau" for "N-transform". In addition, to evaluate the flux of polarized radiation, we also need the "$\Lambda^{(4)}$-transform" (the integral against the E4 kernel), given by "4 MTsi tau". The messy details of the equations behind all this are given in these notes.pdf.

If the run of the Planck function with optical depth, B(tau), and the fractional scattering, lambda(tau), is given, we can use the matrix representations for an iterative solution of the source functions and resultant polarization of the emergent radiation. This code obtains the solution by iterating to convergence: s_and_p0.ijs. Of course, the set of linear equations can be solved directly: s_and_p.ijs. Once the source terms s(tau) and p(tau) have been determined, the total flux at each tau point is given by "Flx=. (F mm s)+(F4 mm p)", where "F=. 2 MTsi tau" and "F4=. 4 MTsi tau". Here, mm is matrix multiplication, defined by the the J expression "mm =: +/ . *"

A series of commands given in the file pol_example.ijs (or pol_example2.ijs) demonstrates how these routines can be used to find the polarization of the emergent radiation.

For the case of frequency-independent absorption and scattering - the "grey atmosphere" case - we may use the condition of radiative equilibrium to show that B(tau) = s(tau) for the integrated radiation. This problem can be solved by forming the matrix equations for the unknown vector {s,p}. The routine in Grey_si.ijs demonstrates how this can be done in J.

While the generation of the matrix transforms is not very quick for a fine grid of optical depths, once they have been computed the solution of the equations is fast and the transforms need not be recomputed unless the optical depth grid is changed.


Slab Geometry

The foregoing code has assumed a semi-infinite atmosphere. We might wish, however, to consider a slab of finite optical depth. In that case, the integration is not extrapolated to infinity, but simply stops at the lower boundary. The routines to generate the matrix transforms are then a bit simpler: MTs.ijs. Also, the integration of the emergent radiation needs a different matrix: IQMs.ijs. And example of the use of these routines is given here: pol_example3.ijs.

In many cases, the slab problem will be symmetric about the central plane. If the slab has an optical thickness of 2T, all the functions will be symmetric about tau = T. So we need only compute functions over the range [0,T], resulting in matrices 1/4 the size of the method considered above. Due to the symmetry, the source function, etc. should have a zero derivative at tau=T, and we can require this of the spline fit. The code to construct the matrices for this symmetric slab case are MTss.ijs and IQMss.ijs . The same example using these transforms is here: pol_example4.ijs.


Monte Carlo Methods for Radiative Transfer Problems: Isotropic Scattering

Another approach to transfer problems is with Monte Carlo methods. Here is the code for the radiation emerging from a uniform source in a slab of finite thickness T: Monte_Uslab.ijs. The same problem can be solved using the matrix methods above: Uniform.ijs. This comparison is useful to study the statistics of the convergence of the Monte Carlo results to the exact solution as a function of the number of photons, the number of times scattered, etc.

We can also examine the case where the sources are not uniformly distributed, but rather are all emitted from the mid-plane of the slab. The Monte Carlo code is MC_mid-slab.ijs, while the matrix approach is Mid_plane.ijs.

Another problem of interest is the scattering of an external beam incident at some angle cosine mu_i. The Monte Carlo code in J is MC_beam.ijs, while the matrix transform J code for the same problem is Beam.ijs.

The equations underlying this code are discussed in these notes2.pdf.

Of course with the Monte Carlo runs we can follow the diffusion of the photons in x,y,z coordinates;
here is an example for a point source in the mid-plane of a slab: Point_in_Slab.ijs.
(Or, this slightly more general code, where we may place the source at any distance below the surface: Point_and_Slab.ijs.)

Here is a note on the Henyey-Greenstein and Rayleigh phase functions, which are often used in Monte Carlo calculations to explore non-isotropic scattering.
The codes HG_slab.ijs and Rayleigh.ijs implement these phase functions for the point-in-slab problem.

Monte Carlo Methods for Radiative Transfer Problems: Scattering with Polarization

Including polarization in the Monte Carlo routines makes things significantly more complex, as we must follow the transformation of the Stokes parameters at each scattering. Some of the equations that are relevant are presented in these notes3.pdf. Here we give implementations in J of the problem of emission from a plane embedded in a slab which scatters with a Rayleigh phase function (slabP.ijs) and emission from sources distributed uniformly in the slab (slabU.ijs).
The problem of uniform emission can also be solved by the matrix methods discussed above. Here is the code for that approach: (U_sp.ijs). Here is a comparison between the two methods for the polarization from a slab of optical half-depth 2 (MC_vs_matrix.pdf). The solid curve is the (essentially exact) matrix solution, the dots the Monte-Carlo results. The Monte-Carlo code followed groups of 300,000 photons through 50 scatterings; this was repeated ncycle=100 times for a total of 3e7 photons. Note that the maximum polarization only reaches 0.041. The total intensity is easy, but the statistics must be much better to get the polarization right. Our Monte-Carlo results only start to deviate for the smallest values of $\mu$ -- rays just grazing the surface -- where the relative number of emerging photons becomes very small.

More interesting is the problem of a point source within the slab. In this case, we cannot use the matrix methods at all. Here is a generalization of the "Point_and_Slab" code given above to the case of Rayleigh scattering with polarization: PPaS.ijs. Here's a screenshot of the intensity (colors) and polarization (line segments) of radiation from a slab of total optical thickness 5 with a point source at the midplane. This is the view looking in at 58 deg from the normal (mu=0.525) downward along the y-axis: inten_and_pol.jpg.

If we average the escaping radiation generated by PPaS.ijs over the face of the slab, we should get the same result as generated by slabP.ijs. Here is such a test, where the total optical thickness of the slab is 4 and the point source is located at z=1 (i.e., 1 unit from the top face and 3 from the bottom). The polarization from the bottom is not too far from the semi-infinite case (nearly all negative), while that from the top face, where the radiation field is mostly horizontal, is mostly positive. Here is the intensity Ave_Inten.pdf and the polarization Ave_Pol.pdf.

A variation of the PPaS.ijs code which uses the rejection technique (see notes3 above) to find the distribution of scattering angles is given here: PPaSrj.ijs.

Line Transfer Problems: The Coupled Escape Probability (CEP) Method

An interesting method for the solution of line transfer problems was introduced by Elitzur and Ramos, MNRAS 365, 779 (2006), hereafter referred to as ``ER06''. Their idea is to divide the medium into zones and compute the probability that a photon emitted is one zone will escape that zone and be absorbed in some other zone. This is simplest in a plane parallel medium divided into slabs. If we assume complete redistribution for the line scattering, the key function needed to compute the CEP coupling matrix is called \alpha(\tau)(ER06 eq 22). This bit of code evaluates it: Alpha. Since it is a smooth function of one variable, it is best to tabulate log(alpha) as a function of log(tau), and interpolate as needed. For our purposes, we have tabulated log(alpha) for 161 values, along with the second derivatives needed for a spline interpolation. (This is that data: log_alpha_table_161.dat[It's in form suitable to the J "readB" verb.]) From this we can easily (and quickly) construct the M_ij matrix given by equation (18) of ER06: MMs.
The use of the CEP matrix is most simply illustrated by application to the classic problem of the scattering line associated with a two level atom. In this case the equation for the source function S is just S = (1-\eps)J + \eps B (ER06 eq 23), where J is the mean intensity, B the Planck function, and eps the ratio of collisional de-excitations to radiative decays. The medium is divided into slabs (i=1,2, ... N) with the S, eps and B assumed constant within each slab. We then obtain a linear equation (ER06 eq 27) for the S_i:

S_i + (\eta_i/dtau_i)* \Sigma_{j=i}^N M_ij*S_j = B_i

Here, \eta=(1-\eps)/\eps and dtau is the difference between tau at the tops and bottoms of the slabs. Here is the J code to solve for the S_i: Two_Level. Once we have the source function we can obtain other quantities, such as the profile of the emerging flux. Here is the J code for the profile (ER06 eq 20): Flux.

This plot shows the source function for B=1 and eps=1e-4. (log S vs log tau). And here is the emergent line profile.
(Which I think bears a strange resemblance to the Minoan symbol that Evans called the ``horns of consecration''!)

More interesting is the problem of multilevel atoms/ions. We implemented the problem used as an example in ER06, the fine structure lines of [O I] in cool neutral clouds where these lines (at 63 and 145 microns) may become optically thick. The method results in non-linear equations which are solved by Newton-Raphson iteration. Here is the J code:OI_3_level. This code takes a couple of seconds to make a model with ~50 zones. It needs some data on the collision cross-sections of oxygen with atomic hydrogen: The data is OI_H_fine_coll.dat. (This data is in a from that can be read by the J function " readB=: (1) 3!:1^:_1 [: 1!:1 < ". It's not plain text, but should be machine independent.)

Here are some results of a run with an oxygen column density of 10^19. The temperature is 300K and the atomic hydrogen density is 2 10^4. Both lines are optically thick (tau(63) = 67 and tau(145)= 7.6). Here are the populations of the upper two levels, here are the source functions of the two lines (the 145 mc source function multiplied by 2), and here are the emergent profiles of the two lines. In some cases, where the first guess values of the populations are not close enough to the solutions, the code will not converge. But one can start with a convergent case and then slowly change the defining parameters, using the last solution as starting guesses. The code can run in this mode, and if it doesn't converge, you can step back to the previous run.
One curious feature of the OI fine structure levels is that the populations may become inverted for temperatures above 66 K at sufficiently low densities. The inversion occurs for n_H < 2e4 at 100K, n_H < 6e4 at 200K, n_H < 9e4 at 300K, etc. This happens because the radiative decay rate from the 2nd level is over 5 times higher than that from the 3rd level. However, if the column density of neutral oxygen is high, trapping of the 67 micron line radiation will counter the radiative drain from the 2nd level, suppressing the inversion. In such cases, inversion will occur first at the slab edge where the escape is greatest. Population inversions will produce laser activity in the 145 micron line; the code given here will not converge under such conditions. In the example given above, n_H = 2e4 at T=300K, inversion would occur at low column density, but for this O column (1e19), the net radiative bracket, even at the slab edges, is < 0.156 and this prevents inversion.


Relevant Maths

This is a Runge-Kutta integration routine for ordinary differential equations: runge_kutta.ijs.

This is an integration routine for stiff systems of differential equations: stiff.ijs.

These are spline interpolation and integration routines: splines.ijs.

If we have a fixed "x" grid but many different functions y(x), it may be useful to precompute a matrix "B" such that (B times y) gives the 2nd derivatives for the spline fit to y(x). This J routine finds such a matrix: SpMD.ijs. SpMD assumes a "natural" cubic spline, i.e., the first and last 2nd derivatives are set to zero. A more general case is given by SpMDD.ijs, where you may specify the first derivative at either or both boundaries. SpMD is useful in constructing matrix representations for the integral of a function against a known kernel function, where the integrals of x^n * kernel(x) are analytic. (See the radiative transfer routines above.)

These are for 2-D bicubic interpolation: bcuint.ijs, and a data file: Bicubic_wts.dat.

The exponential integral function is used frequently in radiative transfer calculations: Ei.ijs.

Here are some utilities to manipulate quaternions: quaternion.ijs. (Quaternions are useful, for example, in rotating coordinates in scattering problems.)

These find the coefficients of the Legendre polynomials and thus find roots and weights for Gaussian quadrature. LegP_coeff.ijs.

You can call LAPACK linear algebraic routines from J.
Here is an example using this to get the roots of a polynomial: poly_root.ijs.