1. Download: Go to http://www.cora.nwra.com/~dbraun/hankel download a gzipped tar file there called Hankel.tar.gz 2. Explanation: The az_hank.c code does both an azimuthal and Hankel transform (I have separate code to do both, but this version does both efficiently without needing to read/write the azimuthal transform). The code tempfftflip.c does the temporal transform and outputs the A and B coefficients (see below). All input and output is via FITS files, the input is 3D and the output is 4D (with the 4th dimension of size 2 indicating real and imaginary components). There is additional code which include the routines for fits I/O (fitssubs.c and fitssubs.h) as well as definitions of keywords in the output FITS file (hank_hdr.h). You will not usually need to vary anything in these source codes. The one exception is params.h which must be modified to specify the details of the annulus and the Hankel transforms. The definition of the Hankel function decomposition is given by equation (4) of Braun 1995 (ApJ 451, 859). In these notes, the terms A and B denote the coefficients of the incoming and outgoing modes which as described by the above citation are functions of m (azimuthal order), L (horizontal wavenumber: L = sqrt(l(l+1)), where l is the spherical harmonic degree), and nu (temporal frequency). The argument of the Hankel functions is L * theta, where theta is the angular distance (in radians) from the target point. This code az_hank.c uses Hankel functions as approximations for the appropriate Legendre polynomials, and therefore is best suited for small annuli centered around local regions, like many sunspot seismology applications. To set (and change) the parameters of the annulus, you need to edit the file params.h and then recompile. Currently the file params.h contains the following definitions: -----------------params.h------------------------------------------ #define TSAM 60. #define NTHET_A1 0 #define NTHET_A2 200 #define DTHET 5.93412e-4 #define NPHI 1258 #define MMAX_A 50 #define NTHET_H1 24 #define NTHET_H2 200 #define MMAX_H 50 #define NLMIN 1 #define NLMAX 88 ------------------------------------------------------------------- TSAM is the interval in time between images (frames) in seconds. NTHET_A1, NTHET_A2, and DTHET define the annulus over which the polar coordinate interpolation and azimuthal transform is computed. The first two quantities are integers defining the inner and outer radii over which the polar coordinate system interpolation is performed, with the spacing in theta given by DTHET (in radians). Note that setting DTHET to be the pixel size of the input data is generally a good idea. To maximize the flexibility of the software, the subsequent Hankel components can be computed over the same, or a smaller subset, of the annulus from which the azimuthal transform is computed. The inner and outer annulus radii for the Hankel transform is is defined by NTHET_H1 and NTHET_H2 respectively. In the above example, the outer radius is 200 pixels for both the azimuthal and Hankel transforms, but the inner radius is 0 and 24 for the azimuthal and Hankel transforms respectively. The selection of the inner and outer annulus radii requires some consideration (see Braun 1995) which, for the sake of brevity, I will not get into here. NPHI= number of pixels in the azimuthal direction to do the polar interpolation. To avoid undersampling at the outer annulus radius, make this value at least 2PI*NTHET_A2. Note it is not necessary to set the spacing in the phi direction, since it is assumed that the NPHI values will cover the full range of 2*PI radians and the spacing will be determined accordingly. MMAX_A= maximum value of the azimuthal order m to compute in the azimuthal transform. Again, a smaller set of m's can be subsequently computed in the Hankel transform (given by MMAX_H) NLMIN and NLMAX are integers denoting the limits of the wavenumbers in the Hankel transform to output (i.e. the set of integers denoted by "j" in Eq. 6 of Braun [1995]) Note that j=0 has no real meaning, and usually I use NLMIN = 0. In general, you shouldn't exceed the spatial Nyquist wavenumber. Therefore NLMAX is usually set to (NTHET_H2-NTHET_H1)/2, which corresponds to the Nyquist wavenumber L_nyq = PI/DTHET. 3. Compile: Only 2 executables are needed: az_hank and tempfftflip although all of the *.c and *.h files are needed and should be in the same directory. You will need the fftw libraries (see www.fftw.org if you never heard of this). I use the gcc compiler: gcc az_hank.c -O3 -lm -fno-builtin -o az_hank gcc tempfftflip.c -O3 -lm -lfftw3 -fno-builtin -o tempfftflip 4. Run: To test if the code compile and runs on your system, compile and run the code (without modification of params.h) on the included test file I0_10_in which is a 3D file which describes an m=0 ingoing Hankel function of a single wavenumber. First run az_hankel and then tempfftflip: az_hankel I0_10_in HA_tran 450 250 tempfftflip HA_tran_in AB_coeff If all works well, the files this creates "HA_tran" and "AB_coeff" should be identical (or nearly so) to the output files I've included "HA_tran_0_10_in" and "AB_coeff_0_10_in" To use these codes in general: Usage: az_hankel INPUT Hankel_trans X0 Y0 [nframe_start nframe_end] INPUT = 3D fits file containing time sequence of NAXIS3 input images, each of size NAXIS1 by NAXIS2, where NAXIS[1-3] must be specified in header. Can be either of type BITPIX=16 or -32 Hankel_trans = output file (to be used by tempfftflip). Can be deleted after running tempfftflip. It is of only limited use by itself. If you are really interested in knowing the format, ask! X0,Y0 = float values of the target position in pixels (assuming pixels start with first row or column = 0) optional parameters: nframe_start and nframe_end are integers which define the range of images in the 3D input to use (0=first image). Leaving them out will cause all images to be processed. To maximize efficiency of the subsequent call of tempfftflip, the total number of images, (nframe_end+1-nframe_start), should be factorable but does not need to be a power of 2. Usage: tempfftflip Hankel_trans AB_coeff_file Hankel_trans = temporary file created by az_hankel (see above). AB_coeff_file = output file containing A and B coefficients as a 4D FITS file. The format is described below. 5. format of AB_coeff_file This is a 4 dimensional fits file, where the 4th dimension is always of size 2 (NAXIS4=2) indicating the real (first element) and imaginary (2nd element) of the output. The other 3 dimensions are as follows (with reference to keywords in the header). We use the convention here that all elements are indexed from "0". NAXIS1: the temporal frequency axis, with units given by DAXIS1 in units of Hz, or by the keyword DELTANU in units of mHz. The k-th element, in the first half of the output, corresponds to the A coefficient for the frequency k/TSAM. The second half stores the B coefficient in backwards order (as in the positive and negative frequencies of a typical FFT). Thus, the A and B coefficients corresponding to the frequency k/TSAM is given by elements k and (NAXIS1-k) respectively Note that the 0th (the “DC”) and NAXIS1/2-th (the Nyquist frequency, when n is even) elements of the output should be disregarded, as their meaning in terms of inward and outward propagating functions is ambiguous (I should just reset these to zero someday). NAXIS2: the azimuthal order, with units of 1 (DAXIS2 will always be 1). The order is m=0,1,...MMAX_H, -MMAX_H, -MMAX_H+1, ...-1 (i.e. the 2nd half are the negative "m"s in backwards order). There is a keyword MAXM in the header which is the same as MMAX_H. The dimension of this axis is 2*MAXM+1 NAXIS3: the wavenumbers, as given by the set defined by Eq. 6 of Braun 1995. The spacing in L is given by DAXIS3, as well as the keyword DELTAL. The first element is always L=0 (in spite of the fact that you specify NLSTART > 0!), but is always set to zero as this value has no meaning in terms of a separation of inward and outward wave components.