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 (e.g., inversion of a matrix is a primitive). 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. At present, J is *free*, and runs under Linux, OS X, and (even) Windows. Get it at http://www.jsoftware.com.

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

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

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.

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.

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.