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. J is *free*. It runs under Linux, OS X, and Windows. The source is available under GPL 3. Get it at http://www.jsoftware.com.

AND... You can run J on your iPhone! Just go to the app store and search on "jsoftware inc". It's a free download. This app is the full J language; all the code on this page (except a few cases where I call a FORTRAN routine) will run on an iPhone. (On my iPhone 4S, J could invert a 500 x 500 matrix in 5.5 sec - and that's in double precision. On my iPhone 6, it only takes 0.5 sec.)

Here are some definitions needed for the following code: 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.)

What Aulus Gellius (ca 160 CE) wrote about the study of logic applies nicely to the study of J:

"... and I need only add by way of advice, that the study and knowledge of this science in its rudiments does indeed, as a rule, seem forbidding and contemptible, as well as disagreeable and useless. But when you have made some progress, then finally its advantages will become clear to you, and a kind of insatiable desire for acquiring it will arise; so much so, that if you do not set bounds to it, there will be great danger lest, as many others have done, you should reach a second childhood amid those mazes and meanders of logic, as if among the rocks of the Sirens."

"Attic Nights" Bk XVI, viii, 15-17.


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

This is for the integration of the Lane-Emden (polytrope) equation: polytrope.ijs.
and this is for the case of the isothermal sphere: iso_sphere.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 evaluates the Voigt function, which describes the shape of spectral lines: voigt.ijs.

(Keplers_Equation.ijs) - This J code solves Kepler's Equation: KE=: ,~ (p. sin)^:_ ] (Yes, that is the entire program!)

Here's a routine to evaluate the Chandrasekhar H(mu) function: H_funct.ijs. (With this function we may evaluate the intensity of light scattered in any direction from a thick, isotropically scattering atmosphere illuminated by radiation incident at an arbitrary angle. For Rayleigh scattering of a beam with polarization see "ARay.ijs" below.)

This code gives the location in the sky of the sun, moon & planets. It also displays the moon's shape
and orientation: Sky.ijs. It needs some data files: SS_DAT.
Here's a new version that plots the location of the sun, moon and planets against the background of all stars brighter than 5.2 magnitudes: Skys.ijs.

This little routine gives the times at which the Milky Way crosses the observer's meridian and it's angle (azimuth): Milky_Way.ijs. Why? A tool for investigating the possible alignments of ancient structures, e.g., The Great Hopewell Road.

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.


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. (See also Harrington (1969).) 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.

The generation of the matrix transforms takes less a second for a fine grid of 80 optical depths. Once they have been computed the solution of the equations is very fast and the transforms need not be recomputed unless the optical depth grid is changed.

Cool Stellar Atmospheres

As a practical example, let's calculate the polarization from a cool stellar atmosphere. From the MARCS website (http://marcs.astro.uu.se/), we get the parameters for a 2500 K, log g = 3.0, solar abundance model. Taking a wavelength of 5000A, we compute the run of the Planck function B_nu(T) with optical depth (the red line is log B_nu). From the MARCS opacity tables we find the values of absorption/(scattering+absorption)= $\lambda$ (blue line). Here are the model_parameters plotted against log(tau). Using the code given in "pol_example2" above, we then find the following results for the source terms s(tau)[blue line] and p(tau)[red line, scaled up by 100], plotted here against log(tau): model_s_and_p. From s and p, we easily obtain the intensity I(mu)[blue line] and polarization fraction [red line, scaled up by 100]: model_I_and_pol, plotted as a function of mu = cos(theta). We see that the polarization only exceeds 1% for mu<0.08, quite close to the limb. This is because scattering only becomes important at small optical depths (< ~ 0.01) in this model. Here is the I(mu) and percent polarization for 4200A, where the Rayleigh scattering is stronger: model_I_and_pol_4200. The polarization now exceeds 1% for mu < ~0.18. (The upturn of I(mu) at small mu is real -- the thermal emission at the surface is very small, but at glancing angles we see the scattered radiation.)

With such code it easy to compute the emergent polarization from any of the MARCS atmosphere models. (E.g., with this J code: extract_pol.ijs). Here are files which contain I(mu), Q(mu) and p(mu)=Q(mu)/I(mu) extracted from 100 of the MARCS model atmospheres in the temperature range 2500K - 5000K and with log(g)= 3, 3.5, 4, 4.5 and 5: MARCS I(mu) Q(mu). The polarization is higher for metal-poor atmospheres. Here are results for a set of MARCS models with 1/10 solar abundances: 1/10 solar MARCS I(mu) Q(mu). It should be pointed out, however, that only the continuum absorption has been included in these calculations. Since the cooler stars are blanketed by line absorption such that there is little or no free continuum, the polarization given here will be overestimated.

If you want to use this method, but are (inexplicably :-) adverse to using "J", here is a directory with a bit of Fortran code and some necessary data (i.e., the pre-computed matrix transforms): For_Pol.d.

Such polarization has been discussed recently in connection with the detection of exoplanets (e.g., Davidson et al. 2010, AAS meeting No. 215 [http://adsabs.harvard.edu/abs/2010AAS...21542303D]; Carciofi & Magalhaes 2005, ApJ, 635, 570). Here is a bit of code that computes the polarization from a transiting planet, given the limb darkening I(mu) and Stokes Q(mu): transit. Here are some results:T2500_4600A.jpg. The scattering and hence polarization will be largest for early type (electron scattering) and late type (Rayleigh scattering) atmospheres, but not very significant for solar type stars. Here are the results we obtain for the transit of a Jupiter size planet across the MARCS solar atmosphere. The light curves are thus:sun_light. Here is the polarization: sun_pol. The maximum polarization is 1.8e-6 (0.00018 %), largely in agreement with the results of Carciofi & Magalhaes (2005).

The polarization from cool atmospheres increases strongly towards shorter wavelengths. This is because (1) Rayleigh scattering increases as lambda^(-4), and (2) the Planck function has a steeper gradient at short wavelengths, resulting in radiation which is more strongly peaked perpendicular to the atmosphere's surface, which in turn leads to higher polarization of radiation scattered near the surface. Here are the polarization curves for the transit of a Jupiter radius planet across a K5 V star with atmosphere parameters T_eff=4500K and log g=4.5. The curves are for wavelengths of 6000A, 5500A, 5200A, 4800A, 4600A, 4400A, 4200A, 4000A and 3800A; the polarization is increasing with decreasing wavelength: K5V_Jupiter.

The discovery of the Kepler-16 system, where a planet orbits a binary system provides a timely example. We take the larger star "A" (R=0.65 R_sun) to have a temperature of 4500K and surface gravity of log g=4.5. Here are the polarization curves for transits of "A" by both star "B" (R=0.22 R_sun) and by planet "b" (R=0.0775 R_sun): Kepler-16. The transit by the companion star produces nearly 5 times the polarization as the planet's transit. (Here are the theoretical light curves: Kepler-16 light curves.)

Hot Stellar Atmospheres

Polarization can also arise in hot stellar atmospheres as a result of electron scattering. Unfortunately, for very hot stars where electron scattering is most important, the gradient of the Plank function with depth is shallow in the visible region so that the radiation field is not strongly peaked outwards. Thus we do not approach the 11.7% polarization Chandrasekhar found for a pure scattering atmosphere. Many extensive grids of hot stellar atmospheres can be found on the web, but for computation of the polarization at a given wavelength, we need the temperature, absorption and scattering as a function of the monochromatic optical depth. For example, the atmospheres tabulated at http://star.arm.ac.uk/~csj/models/Grid.html give the run of the needed variables with depth. Unfortunately, the absorption and optical depth are given only for a wavelength of 4000A. Here is a modification of the "extract_pol.ijs" code given above for use with these atmospheres: extract_hot.ijs. (for cooler stars extract_kurucz)
Let us consider a star of spectral type ~O7, T_eff = 40,000K and log(g)=4.0. At a wavelength of 4000A, the polarization of the emerging radiation is as shown here: T40K_lg4_pol. The maximum polarization is ~3.5% at the limb, but drops quickly, reversing sign at ~83 deg, and reaching a secondary maximum of ~ 0.36% at 69 deg. This reversal in sign corresponds to a switch from polarization parallel to the surface of the atmosphere (-) to perpendicular to the surface (+). As explained in the reference above [Harrington (1970)], this is due to the shallow gradient of the Plank function in the visible for very hot stars.
We might be able to observe such polarization during eclipses of hot stars, though the effects are small. Here are predicted polarization curves for two identical O7 stars as given above: O-star_eclipse.pdf. Three cases are shown for different impact parameters of the eclipse track. We see here that the polarization reversal of the emergent radiation leads to complicated reversals in the eclipse polarization.

Applications to Rotating Stars and Binaries

One application is to evaluate the net polarization integrated over the surface of a rapidly rotating star which will be flattened by rotation and also have a non-uniform surface brightness due to gravity darkening. The (non-realistic) case of a pure electron scattering atmosphere is treated by this code: pure scattering with the gravity darkening following von Zeipel's law. This is an alternate version using a gravity darkening law formulated by Espinosa Lara & Rieutord (2011) pure scattering 2. A more realistic approach is to interpolate in temperature and gravity in the hot star atmospheres above; this code does that: rotating hot star. There is also an effect seen even in slowly rotating stars if we can resolve the narrow spectral lines. Here is an example of code to evaluate this phenomenon with a Milne-Eddington line model: The Öhman effect. This version includes macroturbulence: The Öhman effect No. 2. Another program uses the output of the "2nd_stellar_spectrum" program: Öhman effect No. 3. The theory and results of these programs are discussed in the ``Stellar Polarization'' section of this web page.

Stars in a close binary system will also depart from spherical symmetry, and this could give rise to net polarization. If the stellar shape can be described by a Roche model (mass concentrated at the center), we may calculate the amount of polarization which results. The equations are given on the "Polarization from Tidal Distortion" notes linked on the main web page. Here is the code to do the calculation for an electron scattering atmosphere: Binary.ijs.
And this is the code to use emergent I & Q calculated from MARCS model atmospheres: BinA.ijs.

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.

Polarization of an Incident Beam Scattered by a Plane Parallel Atmosphere

Another problem of interest is the scattering by a slab of an incoming beam of radiation, including polarization. (For example, consider an exoplanet with a scattering atmosphere illuminated by a nearby star. We would like to find the polarization of the scattered radiation.) In this case, the intensity and polarization is a function of two angles that specify the emerging radiation: (1) the cosine of the angle wrt the normal to the slab (mu), and (2) the azimuthal angle (phi) measured from the plane which contains the normal and the incoming beam. (The results will be symmetrical about this plane.) This is the J code to solve this problem: Pol_beam.ijs. This routine also keeps the radiation emerging after each scattering, so that one can see the details of multiple scattering. This also lets you obtain the results for any albedo by weighting the nth scattering by albedo^n.
As an example, here are the results for a slab of total optical thickness 9, illuminated by a beam that enters at mu_0 = 0.7 (~ 45 degrees). We set the albedo to unity and follow 60 scatterings per photon (splitting off the escaped fractional intensity). This plot shows the fractional polarization of all the radiation which emerges from the upper surface of the slab, plotted against the azimuthal angle in degrees and against mu (the cosine of the emergent angle with respect to the slab normal): pol_beam.pdf. And here is the emergent intensity: PB_intensity.
In this case, 82.65% of the radiation emerges from the top of the slab and 9.04% from the bottom. With such a thick slab and albedo=1, the remaining 8.31% is still trapped in the slab after 60 scatterings. Of course even a small amount of absorption would extinguish the trapped radiation.
If we just consider the *first* scattering, we can easily compute the emergent radiation explicitly. Here is the program: Pure_RS. The results from this routine agree well with the Monte Carlo results for radiation emerging after the first scattering. While the polarization from the first scattering reaches 100%, in the case of a thick slab like the tau=9 case above, only 23% of the radiation scattered back is from a single scattering. If we look at the sum of the remaining scatterings (2nd, 3rd, ...) we obtain the polarization from multiple scatterings, which never exceeds 30%. Here, for this case, is the angular distribution of intensity for multiply scattered photons.

These Monte Carlo results can be checked by comparison with the analytic formulae for an illuminated, semi-infinite atmosphere, as derived by Abhyankar & Fymat(1970. A&Ap 4, 101.) and presented by Madhusudhan & Burrows (2012. Ap.J 747, 25.) Here is the J code to evaluate these expressions, which involves iterative solution of several integral equations: ARay.ijs. (In J, iterating expressions to convergence is simple with the power operator ^:_ )
Here is a comparison of the Monte Carlo results with the analytic evaluation: Analytic vs Monte Carlo. The figure shows the fractional polarization as a function of azimuthal angle for 30 values of mu=cos(theta). The curves are indistinguishable except for the extreme values of mu: for mu near zero, relatively few photons emerge in the Monte Carlo runs. (A more general bit of code allows the incoming beam to have arbitrary polarization: ARay_f.ijs.) Of course the analytic expressions take only seconds to evaluate, but they cannot be generalized to the cases of depth varying albedo, non-Rayleigh scattering laws, etc.

Here are some examples for a planet with a thick Rayleigh scattering atmosphere: An illuminated planet:.

All the foregoing can be generalized from Rayleigh scattering to Mie scattering cases. Here we discuss the Monte Carlo code for problems with such more general scattering functions: Code for Mie scattering by particle distributions: And here is one example of a generalization of the Pol_beam code for a case where we have a Rayleigh scattering atmosphere with a layer of haze/dust with a Mie scattering phase function: PB_Mie3.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 in a form suitable to the J "readB" verb and as a plain text file:log_alpha_table_161.ascii) 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 set of linear equations (ER06 eq 27) for the S_i:

Here, delta tau 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 63 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.

The code for the more general problem of an arbitrary number of levels (rather than only 3) is actually somewhat simpler (at least in J!). This is my J code for the general n-level problem implemented for the first 12 rotational levels (and 11 lines) of the CO molecule: CO12.ijs. It is easily applied to other cases. The core subroutines (in J we call them ``verbs'') ABF, ABF0, Iter, p_i, MMs, Alphas, and Esc, as well as most of RunCO are completely general, while SetCO gets the data and sets up the arrays for the CO rotational levels. In this code we call the LAPACK linear equation solver ``dgesv'' (rather than the native J ``%.'') for speed as we are inverting matrices of dimension (# levels)x(# slabs) for each Newton-Raphson iteration (576 x 576 for the case shown below). Here, we first solve the problem for a single zone, and use this as a first guess for the multi-slab problem.
(Here is this code modified for the OI problem above: OI_3_X -- it gives exactly the same results as "OI_3_level" above.)
The rate equations for the level populations are compactly expressed in the central three lines of the function "ABF":

GM=. C_ij +"2 pp*"2 AA
OUT=. II*"1/~ +/"2 ggt*"2 GM
GM=. g 0}"2 (gg*"2 GM)- OUT

- J tends to expose the underlying structure and symmetries of a problem.

The needed data are the energy levels, the statistical weights (just 2J+1), the Einstein A transition rates, and the collision rates with the most important species -- usually hydrogen (atomic or molecular) and/or electrons. Here is the Einstein A matrix: AA. Only the diagonal above the principle diagonal have non-zero values; for other problems, the entire upper triangular matrix could be filled. (Each non-zero entry corresponds to a spectral line.) The collision rates used here consider only (para-) hydrogen molecules. The de-excitation rates are interpolated from Table 1 of Flower and Launay (MNRAS (1985) 214, 271) while the excitation rates follow from detailed balance. Here is the collision rate matrix C_ij for a temperature of 40K and an H_2 density of 10^4 :C_ij. The entries above the diagonal are de-excitation rates while those below the diagonal are excitation rates (all per sub-level). The level energies in eV are just
Elev = 0 0.000476723 0.00143017 0.00286034 0.00476723 0.00715084 0.0100112 0.0133482 0.017162 0.0214525 0.0262198 0.0314637
These data in ASCII are:energy levels, Einstein A values and collision rates with para- H_2.
Here we show the results for a run with a temperature of 40K, H_2 density 1e4 /cm^3, CO column density of 3e17, and the slab divided logarithmically into 48 sub-slabs. Thus the J command is
'xrs N tau SS Ec'=. RunCO fsb;3e17;40;1e4
The results displayed on the screen are CO_output, while the variables xrs, N, etc. hold other output.
Note that if we run our code with just one zone (set slab boundaries 'fsb' to [0 1]), we generate the usual escape probability solution. The populations of the J=0,1,2,3 sub-levels are here plotted as a function of position in the slab (the straight lines are the single zone solutions): pop0123. Here are J=4,5,6,7 :pop4567, and here the log of J=8,9,10,11:pop_8_9_10_11. The populations per sub-level decrease uniformly (J=0>J=1>2>3>4>5... etc) although the total population (2J+1)(sub-level pop) peaks at J=2.
Here pop-over-boltz we plot the populations across the slab divided by the Boltzmann population for 40K. The upper levels are far below the Boltzmann values since the Einstein A's are large, though in the interior of the slab radiative trapping reduces the effective A value. The lower three levels are overpopulated compared to the Boltzmann distribution.
We can look a bit deeper into the influence of the line radiation on the populations by plotting the log of the "net radiative bracket" ("pp" in the code) as a function of the log of "fsb", where fsb is distance from the surface to each slab midpoint divided by the thickness of the whole medium (and thus proportional to the line optical depth for each transition). Here is that plot: CO_NRB. We see that in the interior of the slab the effective radiative decay rates (pp*AA) in the J=2->1, 3->2, 4->3 and 5->4 lines are reduced by over a factor of 100.
Here we plot the emergent flux for the lines:CO_fluxes. The blue bars are this 48 slab solution while the red bars are the single slab approximation.
Here are the emergent line profiles: CO_profiles. (They are plotted in units of the Doppler width of the line. Since this differs from line to line, the areas do not map directly to line fluxes, i.e., the 7->6 line carries more flux than the 6->5 line.)
The line source functions show an interesting variety: CO_SS.
The sum of the line fluxes (+/Ec) is 0.0001481 for the 48 zone model, compared to 0.0001705 for the single zone solution: the single zone approximation overestimates the cooling by 15%.
If the population of an upper sub-level should exceed the population of a lower sub-level, this inversion will lead to maser activity. For CO, this may be a possibility for the J=1 -> J=0 transition. If we look at the stimulated emission correction factor SE = [(pop_lower - pop_upper)/pop_lower] for each line, we see how close that transition is to inversion -- when SE passes 0 and goes negative, the line optical depth becomes negative. Here is ``SE'' for this run: CO_SE. We see that stimulated emission is important for most of the lines, but in particular, the bottom curve, which is the J=1->0 line, ranges from 0.13 at the slab center to 0.10 at the edges.

CEP Method for Spherical Geometries

The Coupled Escape Probability (CEP) method can be applied to cases with spherical geometry. I have written some J code for this purpose. The basic equations are given in these notes:sp_notes. Here is the J code for the two-level atom:2_Lev_sph.ijs. The program uses this pre-computed table: Log_eta_table.dat.

As an example, we present results for a spherical shell with an inner radius of 1 and outer radius of 2. The optical depth along a radius from the center is 5000 and the value of epsilon = 0.0001. The thermal source B is unity throughout the shell. The shell is divided into 63 zones. We compare the resulting source function with that from a plane-parallel slab with an optical half-thickness of 5000 and epsilon = 0.0001 EF="SP_CO12.ijs">SP_CO12.ijs.ijs above). Here is the comparison: T5000eps0001. To see the behavior at the edge of the sphere (slab), we plot the source function against the log of (2 - radius): logT5000eps0001.

Here is another example comparing a sphere with a slab for optical depths of 10 and epsilon=0.01: T10eps01. As expected, the effective optical depth of the sphere is less so the source function is lower. But that is not all: when we raise the radial optical depth of the sphere to 13.85 so the central value of the source function matches the slab, we see that the slopes differ.

Some limited testing indicates that the coupling matrix must be computed very carefully. For these examples, we used 200 mu angles and carried out the mu integrations for 10 positions within each shell, for a total of 126,000 rays.

The code for the OI 3-level problem in plane parallel geometry was given in the section above. Here is the code for the same 3-level problem in spherical geometry: Sph_OI_3_level.ijs. We have also adapted the code presented above for a 12-level CO line problem to spherical geometry: SP_CO12.ijs.

Photoionization of Gaseous Nebulae

Here is some code to evaluate the temperature and ionization of gas illuminated by a blackbody:ENEB.
Elements included are H,He,C,N,O and Ne. Diffuse radiation on-the-spot.
nebula.ijs runs the program. SETUP'' solves for the conditions at the inner edge of the spherical shell, RUN'' integrates into the nebula. EVOL'' follows the recombination and cooling of the inner edge if the radiation field is turned off.

Tool for reading MESA files.

MESA is an open source, state-of-the-art stellar evolution program that you can compile and run
on your Mac or Linux box. The project is called "Modules for Experiments in Stellar Astrophysics".
See the: MESA webpage for complete information about this important project.

The MESA program writes ASCII files giving the details of every 50th model, which are stored in the LOGS directory. The files contain 45 global variables and 90 variables that detail the conditions in each zone of the model (usually >1000 zones). Here is a snippet of code to read all this information into J for plotting and analysis: mesa_read.ijs

Relevant Maths

This is a routine to solve for the root of a (maybe nonlinear) equation, using Brent's method: brent.ijs.

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.)

If we have many data points but they have some associated error, we may want to fit a cubic spline using a least squares approach. Here is code to make and evaluate such least squares splines. This is an example of such a fit: Fit to titanium heat data.

Here is a function for 2-D bicubic interpolation: bcuint.ijs.

The exponential integral function is used frequently in radiative transfer calculations: Ei.ijs.
Here is a better version that works on arrays (it's much faster for typical J applications): EIn.ijs.

Here are some utilities to manipulate quaternions: Quaternions.ijs. (Quaternions are useful, for example, in rotating coordinates in scattering problems.) E.g., for fast rotations: qRot.ijs.

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.

Bernstein polynomials are used for Bezier curves: Bezier.ijs.

Here is code for the Incomplete Gamma Function: GamInc.ijs.