; STFUN
;---------------------------------------------------------------
;! Task to calculate a structure function image
;# TASK ANALYSIS
;-----------------------------------------------------------------------
;;  Copyright (C) 1995, 2009
;;  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
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:
  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)<A(1), and
                           abs(Y)<A(3) ) , for A(1) and A(3) > 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 <F> can be obtained, and an autocorrelation function R(s)
for the irregularities dF(r) = F(r) - <F> 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 <F> makes the calculation of R(s) uncertain, while the structure
function D(s), which is independent of <F>, 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

