; STFUN ;--------------------------------------------------------------- ;! Task to calculate a structure function image ;# TASK ANALYSIS ;----------------------------------------------------------------------- ;; Copyright (C) 1995 ;; Associated Universities, Inc. Washington DC, USA. ;; ;; This program is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation; either version 2 of ;; the License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public ;; License along with this program; if not, write to the Free ;; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, ;; MA 02139, USA. ;; ;; Correspondence concerning AIPS should be addressed as follows: ;; Internet email: aipsmail@nrao.edu. ;; Postal address: AIPS Project Office ;; National Radio Astronomy Observatory ;; 520 Edgemont Road ;; Charlottesville, VA 22903-2475 USA ;----------------------------------------------------------------------- ;--------------------------------------------------------------- STFUN LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC STFUN: Task to calculate a structure function image USERID -32000.0 32000.0 User ID. 0=current user 32000=>all users INNAME Input image name (name) INCLASS Input image name (class) INSEQ 0.0 9999.0 Input image name (seq. #) INDISK 0.0 9.0 Input image disk drive # OUTNAME Output image name (name) OUTCLASS Output image name (class) Default is 'STFUN' OUTSEQ 0.0 9999.0 Output image name (seq. #) OUTDISK 0.0 9.0 Output image disk drive # BLC BLC of area to be used TRC TRC of area to be used IMSIZE 0.0 1024.0 Output image size APARM User supplied parameters: (1) min x-lag (pixels) 0=>0 (2) max x-lag (pixels) 0=> largest possible (3) min y-lag (pixels) 0=>0 (4) max y-lag (pixels) 0=> largest possible (5) >=0 => blank zero-lag pixel, else => zero (6) <=0 =>zero pixels of lag smaller than min x-lag and min y-lag, else => blank (7) <=0 =>zero pixels of lag larger than max x-lag and max y-lag, else => blank see HELP STFUN ----------------------------------------------------------------- STFUN Type: Task Use: STFUN is a task which calculates a structure function image. Output image pixel (X,Y) has the value < [MAP(I,J) - MAP(I+X,J+Y)]**2 > where MAP(I,J) are values in the input image and the angular brackets denote an average over all I,J. The pixel coordinates X,Y are the "lags" in the x and y directions. Currently, the maximum allowed input subimage size is 512 by 512. This may change. The program will complain if the input subimage is too large. Adverbs: USERID......The ID of the owner of the images. 0 => current user. 3200 0=> all users. INNAME......First image name (name ). blank => any. INCLASS.....First image name (class). blank => any. INSEQ.......First image name (seq. #). 0 => any. INDISK......Disk drive # for the first image. 0 => any. OUTNAME.....Output image name (name). Standard defaults. OUTCLASS....Output image name (class). Default is 'STFUN' OUTSEQ......Output image name (seq. #). 0 => highest unique. OUTDISK.....Output image disk drive #. BLC.........BLC of area to be used. Standard default. TRC.........TRC of area to be used. Standard default. IMSIZE......Output image size. 0,0 => 2 * min[ (TRC-BLC), (A(2),A(4)) ] + (1,1). APARM.......User supplied parameters. A(1) = x-lag minimum (in pixels). 0 => 0. A(2) = x-lag maximum (in pixels). 0 => TRC(1) - BLC(1). A(3) = y-lag minimum (in pixels). 0 => 0. A(4) = y-lag maximum (in pixels). 0 => TRC(2) - BLC(2). A(5) >= 0 => blank central pixel (X=0,Y=0), else zero. A(6) <= 0 => zero central square ( abs(X) 0, else blank. A(7) <= 0 => zero outer padded frame ( abs(X)>A(2), or abs(Y)>A(4) ), else blank. ----------------------------------------------------------------------- STFUN: Task to compute a structure function image DOCUMENTOR: John H. Simonetti NRAO RELATED PROGRAMS: COMB PURPOSE ------------- A structure function is useful in analyzing the random irregularities present in an image. It can be used to discover the angular correlation scale of the irregularities, and learn something about their power spectrum. It is most useful when the correlation scale of the irregularities is larger than the image size (as it may be, for instance, when one looks at a Rotation Measure map where the irregularities are due to a foreground screen). For a map of the quantity F(r), where the angular correlation scale for the irregularities dF is smaller than the image size, an accurate mean can be obtained, and an autocorrelation function R(s) for the irregularities dF(r) = F(r) - can calculated as R(s) = < dF(r) dF(r+s) >. The variables r and s are vector angular separations within the map, and the angular brackets denote an average over r. Under these circumstances the structure function D(s) = < [F(r) - F(r+s)]**2 > and the autocorrelation function are simply related: D(s) = 2 v [ 1 - R(s)/R(0) ], where v = R(0) = < dF**2 > is the variance of the irregularities. For a map F(r) where the angular correlation scale of the irregularities is larger than the image or subimage, having no handle on makes the calculation of R(s) uncertain, while the structure function D(s), which is independent of , is still calculable, and to within the random uncertainties due to estimating each point of D(s) from a limited number of map pairs, the calculated structure function will approximate that found from a larger image with the same dF statistics. The power spectrum of the map irregularities is usually related to the structure function in a simple way (just as the power spectrum and autocorrelation function are fourier transforms of each other). The exact relationship will depend upon the how the variable you are measuring in the map has been produced, and is something you should think about. Generally, a gaussian power spectrum would yield D(s) proportional to s**2 for s < the correlation angular scale of the irregularities. You should expect D(s) to be proportional to s**2 for s smaller than the HPBW of the resolving beam for the input map. The behavior of D(s) is most easily seen in a slice of the D(s) map from the center pixel outward along an arbitrary radius. Indeed, it is recommended that slices be made of D(s) maps where the central zero-lag pixel (X=0,Y=0) has been blanked, and these slices be plotted on a log-log scale, if possible (which would have to be done outside AIPS in its current form). Different behaviors for different radial directions would indicate the irregularities have anisotropic statistics. COMMENTS ------------ GENERAL: STFUN uses quite a bit of cpu time for large input subimages and where the structure function is calculated over a large range of lags (see EXECUTION TIMES). It is recommended that a structure function analysis of an image be approached conservatively by suitable blanking of unimportant parts of the input image, choice of BLC and TRC, and the APARMS. In particular, blanking of unimportant parts of the image to leave out regions that should not be included in the statistical analysis, will speed up the process considerably. The program also runs faster if structure function values are calculated only for a lags much smaller than the maximum possible lag within the input subimage. This should be adequate when the large lag values would not contain any useful information as when the correlation scale of the irregularities is much smaller than the input subimage size. A useful procedure might be to use a sizeable input subimage (specified through BLC and TRC), with blanking, and first produce a structure function map with values calculated for lags from zero up to some small maximum in x and y (APARM(2) and APARM(4)) that may or may not be comparable to the largest angular scale of the irregularities present in the input map. The result will be a structure function map of limited lag, which uses all the small lag information present in the input subimage. The rest of the output image (controlled by IMSIZE) out to the maximum lags possible in the x and y directions could be filled with zeros. Then if it seems useful to find values for lags larger than APARM(2) and APARM(4) from the first run, a new output map can be made with zeros filling the pixels for lags smaller than the limits set in the first run, and calculated values for larger lag limits (set using APARM(1), (2), (3) and (4)). Again the very largest lag pixels can be given zero values filling the map out IMSIZE. The resulting map from the second run can now be overlayed (using COMB) on the map from the first run giving a complete structure function map out to the new maximum lags. The process can be repeated, or stopped at any time without having taken up too much time to obtain a useful final result. (See EXAMPLE below.) BLC, TRC: Should be used to limit the region of the input map used in the statistical calculation. Currently, the input subimage must be smaller than 512 by 512. Your machine may allow larger dimensions, if so, change the program as directed in the precursor comments in the code. The program will complain if the subimage is too large. IMSIZE: Optionally used to increase the size of the output map beyond the size dictated by the maximum x and y lags specified either by the input subimage size, or by APARM(2) and APARM(4). This could be useful when making a structure function image with small lags only (smaller than the input subimage size). If the resulting image indicates a new image with larger lags might be useful, the final two images can be pieced together easily only if they are of the same size. IMSIZE = 0,0 and APARM(2) & APARM(4) = 0 In this case, the output image size is set by the maximum lag in the input window. IMSIZE =0,0 and the APARMS are non-zero. In this case, the output image size will (2*APARM(2) + 1) by (2*APARM(4) +1) You can of course default just IMSIZE(1) and specify IMSIZE(2), or vice-versa. If IMSIZE is set too small given the APARM(2) and APARM(4) values, it will be automatically increased to allow for the APARM values. IMSIZE should really be an odd number, if you give it an even number STFUN will increase it by one and tell you about it. APARMS: The APARM(1) through APARM(4) parameters are used to select the lag space (in pixels) for which structure function values will be calculated. The defaults are to calculate values for all lags possible from 0 to a maximum x-lag and y-lag, equal to the x any y width - 1 of the input subimage. In general, the APARMS can be set to calculate structure function values for pixels within a rectangle. APARM(5) is used to blank the central, zero-lag pixel of the output map independent of APARM(6). This is the default and is useful if you plan to make radial slices outward from the central pixel and plot the results on a log-lag scale (otherwise the central pixel has a zero value). APARM(6) and APARM(7) are used to zero the pixels in the output map for which structure function values have not been calculated. The default zeroing of these pixels allows for later overlaying of output maps as described above under GENERAL comments, and in the following EXAMPLE. EXAMPLE OF THE USE OF IMSIZE, AND APARMS (1) THROUGH (4): If the input subimage is 100 by 100 pixels in size, the maximum lag possible is 99 in the x and y directions. If we choose to calculate the structure function for x and y lags from 0 to 30 (in pixels) we set APARMS(1-4) as 0,30,0,30. If we wish the output map to be only large enough to show these calculated values, use IMSIZE 0,0. Then the output map will have size 61 by 61 (61 = 2 * 30 + 1, which takes into account the positive and negative lags, and zero lag). If we wish to extend the output map to the maximum size it would be if we were to calculate the structure function for all possible lags, we would use IMSIZE 199,199 (199 = 2 * 99 + 1). In this latter case we might zero the pixel values for lags beyond 30 using APARM(7) = 0. If we then wished to calculate the structure function for lags larger than 30 but no more than 50 in the x and y directions, we would repeat the task with APARMS(1-4) set to 31,50,31,50. If IMSIZE was still set to 199,199, and APARM(6) = 0, the new result can be added to the old result using COMB to obtain a complete structure function map for lags from 0 to 50. UNCERTAINTIES IN THE CALCULATED VALUES: If we take the map value at any point to be one possible sample from an ensemble of samples that would be normally distributed, and the map values are uncertain due to a measurement error of sigma, then the variance var(x,y) for an output structure function value D(x,y) (about the ensemble mean structure function) can be approximated by / [2/N(x,y)] D(x,y)**2, for D(x,y) > 2 sigma**2, var(x,y) = | \ [2/N(x,y)] 4 sigma**4, for D(x,y) < 2 sigma**2, where N(x,y) is the number of pixel pairs in the imput subimage used to calculate D(x,y). The second line is included since D(x,y) may be less than 2 sigma**2, due to N(x,y) being small or where s = sqrt[x**2 + y**2] is less than the resolving beamwidth. An average D(x,y) over a complete ensemble of maps with similar statistics will be greater than 2 sigma**2, as long as s is greater than the resolving beamwidth. These expressions are only valid if lags much greater than the angular correlation length were sampled in the image. If this is not the case, you should expect the uncertainties will be larger than given above. In using the above expressions, it may be necessary to use an estimate for sigma when sigma varies over the imput subimage. EXECUTION TIMES: This program does not use Q routines (AP code), although it probably could. Thus, its a bit of a clunker on the VAX. However, it has been optimized on the Convex C-1 and it is crucial that the program is installed with the vectorizing option turned on. Roughly a factor of 50 gain in speed is obtained between the local scalar (O0) and vector (O2) optimization levels on the C-1. On the VAX 11/780, a factor of 2 in speed was picked up by turning on the optimizer. Note that the standard AIPS installation uses /NOOPT on the VAX and the O0 level on the Convex. The following guide is for a square input image of width N pixels (none blank) and for calculating values at all possible lags. The CPU time is T seconds. For small images, all the overhead is in the AIPS machinery of opening images and writing them etc. However, because the number of operations goes roughly as N**4 in the main work routine, the structure function computation rapidly becomes the dominant factor. VAX 11/780 / NOOPT ------------------ N <= 25 T <= 15 N > 25 log T = 3.8 * log N - 4.1 until page swapping kills it VAX 11/780 / OPT ----------------- N <= 25 T <= 7 N > 25 log T = 3.9 * log N - 4.8 until page swapping kills it CONVEX C-1 level O0 optimization -------------------------------- N <= 25 T <= 7 N > 25 log T = 3.9 * log N - 4.8 Almost identical to the optimized VAX 11/780 CONVEX C-1 level O1 optimization -------------------------------- N <= 35 T <= 10 N > 35 log T = 3.9 * log N - 5.1 CONVEX C-1 level O2 optimization -------------------------------- N <= 60 T <= 10 N > 60 log T = 3.6 * log N - 5.7