home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdfft.c < prev    next >
C/C++ Source or Header  |  1995-03-13  |  43KB  |  1,301 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdfft.c    2.24  3/13/95
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  *  grdfft.c
  9.  *
  10.  * A program to do various operations
  11.  * on grdfiles in the frequency domain.
  12.  *
  13.  * W.H.F. Smith, 7 March 1990.
  14.  *
  15.  * Version 2.0.2 has option to do power spectral estimates.
  16.  * added by WHFSmith, 4 Feb 1992.
  17.  *
  18.  */
  19.  
  20. #include "gmt.h"
  21.  
  22. #ifndef FSIGNIF
  23. #define FSIGNIF 24
  24. #endif
  25. #define UP_DOWN_CONTINUE 0
  26. #define AZIMUTHAL_DERIVATIVE 1
  27. #define DIFFERENTIATE 2
  28. #define INTEGRATE 3
  29. #define ISOSTASY 4
  30. #define FILTER 5
  31. #define SPECTRUM 6
  32.  
  33. /* Macro definition ij_data(i,j) finds the array index to an element
  34.     containing the real data(i,j) in the padded complex array:  */
  35.  
  36. #define    ij_data(i,j) (2*(nx2*((j)+j_data_start)+(i)+i_data_start))
  37.  
  38.  
  39. int    narray[2];
  40. int    ndim = 2, ksign, iform = 1, i, j, k, n, nx2 = 0, ny2 = 0, ndatac, i_data_start, j_data_start;
  41. int    pad[4];
  42.  
  43. BOOLEAN    map_units = FALSE, force_narray = FALSE, suggest_narray = FALSE, leave_trend_alone = FALSE, n_user_set = FALSE;
  44.  
  45. double    data_var, data2_var, delta_kx, delta_ky;
  46. double    a[3];    /* Plane fitting coefficients  */
  47. double    mGal_at_45 = 980619.9203; /* Moritz's 1980 IGF value for gravity in mGal at 45 degrees latitude */
  48.  
  49. char    *infile = NULL, *outfile = NULL;
  50.  
  51. struct GRD_HEADER h;
  52. struct F_INFO {
  53.     double    lc[3];        /* Low-cut frequency for r, x, and y    */
  54.     double    lp[3];        /* Low-pass frequency for r, x, and y    */
  55.     double    hp[3];        /* High-pass frequency for r, x, and y    */
  56.     double    hc[3];        /* High-cut frequency for r, x, and y    */
  57.     double    ltaper[3];    /* Low taper width for r, x, and y    */
  58.     double    htaper[3];    /* High taper width for r, x, and y    */
  59.     int    do_this[3];    /* T/F this filter wanted for r, x, and y    */
  60.     int    set_already;
  61. } f_info;
  62.  
  63. struct FFT_SUGGESTION {
  64.     int    nx;
  65.     int    ny;
  66.     int    worksize;    /* # single-complex elements needed in work array  */
  67.     int    totalbytes;    /* (8*(nx*ny + worksize))  */
  68.     double    run_time;
  69.     double    rms_rel_err;
  70. }    fft_sug[3];    /* [0] holds fastest, [1] most accurate, [2] least storage  */
  71.  
  72.  
  73. float    *datac, *workc;
  74. double    scale_out = 1.0;
  75.  
  76. main (argc, argv)
  77. int    argc;
  78. char    **argv;
  79. {
  80.  
  81.     int    op_count = 0, n_op_count = 0, par_count = 0, n_par = 0, *operation = NULL;
  82.     BOOLEAN error = FALSE, stop, chose_spectrum = FALSE, give_wavelength;
  83.     double    *par = NULL;
  84.  
  85.     for (i = 0; i < 3; i++) {
  86.         f_info.lc[i] = f_info.lp[i] = -1.0;    /* Set negative, below valid freqency range  */
  87.         f_info.hp[i] = f_info.hc[i] = MAXDOUBLE;    /* Set huge positive, above valid freqency range  */
  88.         f_info.ltaper[i] = f_info.htaper[i] = 0.0;    /* 1/width of taper, zero when not used  */
  89.         f_info.do_this[i] = FALSE;
  90.     }
  91.     f_info.set_already = FALSE;
  92.     
  93.     argc = gmt_begin (argc, argv);
  94.  
  95.         for (i = 1; i < argc; i++) {
  96.                 if (argv[i][0] == '-') {
  97.                         switch (argv[i][1]) {
  98.                         
  99.                 /* Common parameters */
  100.             
  101.                 case 'V':
  102.                 case '\0':
  103.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  104.                     break;
  105.                 
  106.                 /* Supplemental parameters */
  107.             
  108.                                 case 'A':
  109.                     n_op_count++;
  110.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  111.                                         operation[op_count] = AZIMUTHAL_DERIVATIVE;
  112.                     op_count++;
  113.                     n_par++;
  114.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  115.                     if ((sscanf(&argv[i][2], "%lf", &par[par_count])) != 1) {
  116.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option:  Cannot read azimuth\n", gmt_program);
  117.                         error++;
  118.                     }
  119.                     par_count++;
  120.                                         break;
  121.                                 case 'C':
  122.                     n_op_count++;
  123.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  124.                                         operation[op_count] = UP_DOWN_CONTINUE;
  125.                     op_count++;
  126.                     n_par++;
  127.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  128.                     if ((sscanf(&argv[i][2], "%lf", &par[par_count])) != 1) {
  129.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Cannot read zlevel\n", gmt_program);
  130.                         error++;
  131.                     }
  132.                     par_count++;
  133.                                         break;
  134.                                 case 'D':
  135.                     n_op_count++;
  136.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  137.                                         operation[op_count] = DIFFERENTIATE;
  138.                     op_count++;
  139.                     n_par++;
  140.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  141.                                         if (argv[i][2]) {
  142.                         par[par_count] = (argv[i][2] == 'g' || argv[i][2] == 'G') ? mGal_at_45 : atof (&argv[i][2]);
  143.                         if (par[par_count] == 0.0) {
  144.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  scale must be nonzero\n", gmt_program);
  145.                             error++;
  146.                         }
  147.                     }
  148.                     par_count++;
  149.                                         break;
  150.                                 case 'E':
  151.                     n_op_count++;
  152.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  153.                                         operation[op_count] = SPECTRUM;
  154.                     op_count++;
  155.                     n_par++;
  156.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  157.                     par[par_count] = 0.0;
  158.                     give_wavelength = FALSE;
  159.                     j = 2;
  160.                     while(argv[i][j]) {
  161.                          if (argv[i][2] == 'x' || argv[i][2] == 'X')
  162.                             par[par_count] = 1.0;
  163.                         else if (argv[i][2] == 'y' || argv[i][2] == 'Y')
  164.                             par[par_count] = -1.0;
  165.                         else if (argv[i][j] == 'w' || argv[i][j] == 'W')
  166.                             give_wavelength = TRUE;
  167.                         j++;
  168.                     }
  169.                     par_count++;
  170.                     chose_spectrum = TRUE;
  171.                                         break;
  172.                                 case 'F':
  173.                                     if (!(f_info.set_already)) {
  174.                         n_op_count++;
  175.                         operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  176.                                             operation[op_count] = FILTER;
  177.                                             op_count++;
  178.                                             f_info.set_already = TRUE;
  179.                                         }
  180.                                         if (parse_f_string(&argv[i][2])) {
  181.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -F option: ", gmt_program);
  182.                                             error++;
  183.                                         }
  184.                                         break;
  185.                                 case 'G':
  186.                                         outfile = &argv[i][2];
  187.                                         break;
  188.                                 case 'I':
  189.                     n_op_count++;
  190.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  191.                                         operation[op_count] = INTEGRATE;
  192.                     op_count++;
  193.                     n_par++;
  194.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  195.                                         if (argv[i][2]) {
  196.                         par[par_count] = (argv[i][2] == 'g' || argv[i][2] == 'G') ? mGal_at_45 : atof (&argv[i][2]);
  197.                         if (par[par_count] == 0.0) {
  198.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  scale must be nonzero\n", gmt_program);
  199.                             error++;
  200.                         }
  201.                     }
  202.                     par_count++;
  203.                                         break;
  204.                                 case 'L':
  205.                                         leave_trend_alone = TRUE;
  206.                                         break;
  207.                                 case 'M':
  208.                                         map_units = TRUE;
  209.                                         break;
  210.                                 case 'N':
  211.                     if (argv[i][2] == 'f' || argv[i][2] == 'F')
  212.                         force_narray = TRUE;
  213.                     else if (argv[i][2] == 'q' || argv[i][2] == 'Q')
  214.                         suggest_narray = TRUE;
  215.                     else {
  216.                                             sscanf(&argv[i][2], "%d/%d", &nx2, &ny2);
  217.                         n_user_set = TRUE;
  218.                     }
  219.                                         break;
  220.                                 case 'S':
  221.                                         scale_out = (argv[i][2] == 'd' || argv[i][2] == 'D') ? 1.0e6: atof (&argv[i][2]);
  222.                                         break;
  223.                                 case 'T':
  224.                     n_op_count++;
  225.                     operation = (int *) memory((char *)operation, n_op_count, sizeof(int), "grdfft");
  226.                                         operation[op_count] = ISOSTASY;
  227.                     op_count++;
  228.                      n_par += 5;
  229.                     par = (double *) memory((char *)par, n_par, sizeof(double), "grdfft");
  230.                                         n = sscanf (&argv[i][2], "%lf/%lf/%lf/%lf/%lf", &par[par_count],
  231.                         &par[par_count+1], &par[par_count+2], &par[par_count+3], &par[par_count+4]);
  232.                     for (j = 1, k = 0; j < 5; j++) if (par[par_count+j] < 0.0) k++;
  233.                     if (n != 5 || k > 0) {
  234.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -T option.  Correct syntax:\n", gmt_program);
  235.                         fprintf (stderr, "\t-T<te>/<rhol>/<rhom>/<rhow>/<rhoi>, all densities >= 0\n");
  236.                         error++;
  237.                     }
  238.                     par_count += 5;
  239.                                         leave_trend_alone = TRUE;
  240.                                         break;
  241.                                default:
  242.                                         error = TRUE;
  243.                     gmt_default_error (argv[i][1]);
  244.                                         break;
  245.                         }
  246.                 }
  247.                 else
  248.                         infile = argv[i];
  249.         }
  250.  
  251.     if (argc == 1 || gmt_quick) {
  252.         fprintf (stderr, "grdfft %s - Perform mathematical operations on grdfiles in the frequency domain\n\n", GMT_VERSION);
  253.         fprintf (stderr,"usage: grdfft <in_grdfile> [-G<out_grdfile>] [-A<azimuth>] [-C<zlevel>]\n");
  254.         fprintf (stderr,"\t[-D[<scale>]] [-E[x_or_y][w]] [-F[x_or_y]<lc>/<lp>/<hp>/<hc>] [-I[<scale>]] [-L] [-M]\n");
  255.         fprintf (stderr,"\t[-N<stuff>] [-S<scale>] [-T<te/rl/rm/rw/ri>] [-V]\n\n");
  256.         
  257.         if (gmt_quick) exit (-1);
  258.         
  259.         fprintf (stderr,"\tin_grdfile is the input netCDF grdfile\n");
  260.         fprintf (stderr, "\tOPTIONS:\n");
  261.         fprintf (stderr,"\t-G filename for output netCDF grdfile\n");
  262.         fprintf (stderr,"\t-A<azimuth> Take azimuthal derivative along line <azimuth> degrees CW from North.\n");
  263.         fprintf (stderr,"\t-C<zlevel>  Continue field upward (+) or downward (-) to zlevel (meters).\n");
  264.         fprintf (stderr,"\t-D Differentiate, i.e., multiply by kr [ * scale]  [Default scale gives mGal from m]\n");
  265.         fprintf (stderr,"\t-E Estimate spEctrum of r [x] [y].  Write f, power[f], 1 std dev(power[f]) to stdout.\n");
  266.         fprintf (stderr,"\t\tAppend w to write wavelength instead of frequency.\n");
  267.         fprintf (stderr,"\t-F Filter r [x] [y] freq according to four wavelengths <lc>/<lp>/<hp>/<hc>.\n");
  268.         fprintf (stderr,"\t   freq outside <lc>/<hc> are cut; inside <lp>/<hp> are passed, rest are tapered.\n");
  269.         fprintf (stderr,"\t   replace wavelength by - to skip, e.g.  -F-/-/500/100 is a low-pass filter.\n");
  270.         fprintf (stderr,"\t-I Integrate, i.e., devide by kr [ * scale]  [Default scale gives m from mGal]\n");
  271.         fprintf (stderr,"\t-L Leave trend alone.  Do not remove least squares plane from data.  [Default removes plane].\n");
  272.         fprintf (stderr,"\t-M Map units used.  Convert grid dimensions from degrees to meters.\n");
  273.         fprintf (stderr,"\t-N<stuff>  Choose or inquire about suitable grid dimensions for FFT.\n");
  274.         fprintf (stderr,"\t\t-Nf will force the FFT to use the dimensions of the data.\n");
  275.         fprintf (stderr,"\t\t-Nq will inQuire about more suitable dimensions.\n");
  276.         fprintf (stderr,"\t\t-N<nx>/<ny> will do FFT on array size <nx>/<ny> (Must be >= grdfile size).\n");
  277.         fprintf (stderr,"\t\tDefault chooses dimensions >= data which optimize speed, accuracy of FFT.\n");
  278.         fprintf (stderr,"\t\tIf FFT dimensions > grdfile dimensions, data are extended and tapered to zero.\n");
  279.         fprintf (stderr,"\t-S multiply field by scale after inverse FFT [1.0]\n");
  280.         fprintf (stderr,"\t   Give -Sd to convert deflection of vertical to micro-radians\n");
  281.         fprintf (stderr,"\t-T Compute isostatic response.  Input file is topo load. Append elastic thickness,\n");
  282.         fprintf (stderr,"\t   and densities of load, mantle, water, and infill, all in SI units\n");
  283.         fprintf (stderr,"\t   It also implicitly sets -L\n");
  284.         fprintf (stderr,"\tList operations in the order desired for execution.\n");
  285.         exit (-1);
  286.     }
  287.  
  288.     if (!(n_op_count)) {
  289.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one operation\n", gmt_program);
  290.         error++;
  291.     }
  292.     if (n_user_set && (nx2 <= 0 || ny2 <= 0)) {
  293.         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  nx2 and/or ny2 <= 0\n", gmt_program);
  294.         error++;
  295.     }
  296.     if (scale_out == 0.0) {
  297.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  scale must be nonzero\n", gmt_program);
  298.         error++;
  299.     }
  300.     if (!infile) {
  301.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  302.         error++;
  303.     }
  304.     if (chose_spectrum && outfile) {
  305.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  -G ignored (stdout used for ASCII output)\n", gmt_program);
  306.         error++;
  307.     }
  308.     if (!chose_spectrum && !outfile) {
  309.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", gmt_program);
  310.         error++;
  311.     }
  312.     if (error) exit (-1);
  313.  
  314.     if (read_data(argc, argv) ) {
  315.         fprintf (stderr,"grdfft: Fatal memory error.  Exiting.\n");
  316.         exit (-1);
  317.     }
  318.     
  319.     /* Check that no NaNs are present */
  320.     
  321.     stop = FALSE;
  322.     for (j = 0; !stop && j < h.ny; j++) for (i = 0; !stop && i < h.nx; i++) stop = bad_float ((double)datac[ij_data(i,j)]);
  323.     if (stop) {
  324.         fprintf (stderr, "grdfft: Input grid cannot have NaNs!\n");
  325.         exit (-1);
  326.     }
  327.  
  328.     if (gmtdefs.verbose) fprintf (stderr, "grdfft: Forward FFT\n");
  329.     
  330.     ksign = -1;
  331.  
  332.     fourt_(datac, narray, &ndim, &ksign, &iform, workc);
  333.  
  334.     for (op_count=0, par_count = 0; op_count < n_op_count; op_count++) {
  335.         switch(operation[op_count]) {
  336.             case UP_DOWN_CONTINUE:
  337.                 if (gmtdefs.verbose) ((par[par_count] < 0.0) ? fprintf (stderr, "grdfft: Downward continuation\n") : fprintf (stderr, "grdfft: Upward continuation\n"));
  338.                 do_continuation (&par[par_count]);
  339.                 par_count++;
  340.                 break;
  341.             case AZIMUTHAL_DERIVATIVE:
  342.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: azimuthal derivative\n");
  343.                 do_azimuthal_derivative (&par[par_count]);
  344.                 par_count++;
  345.                 break;
  346.             case DIFFERENTIATE:
  347.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: differentiate\n");
  348.                 do_differentiate (&par[par_count]);
  349.                 par_count++;
  350.                 break;
  351.             case INTEGRATE:
  352.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: integrate\n");
  353.                 do_integrate (&par[par_count]);
  354.                 par_count++;
  355.                 break;
  356.             case ISOSTASY:
  357.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: isostasy\n");
  358.                 do_isostasy (&par[par_count]);
  359.                 par_count += 5;
  360.                 break;
  361.             case FILTER:
  362.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: filter\n");
  363.                 do_filter ();
  364.                 break;
  365.             case SPECTRUM:
  366.                 if (gmtdefs.verbose) fprintf (stderr, "grdfft: spectrum\n");
  367.                 do_spectrum (&par[par_count], give_wavelength);
  368.                 par_count += 2;
  369.                 break;
  370.         }
  371.     }
  372.  
  373.     if (outfile) {
  374.  
  375.         if (gmtdefs.verbose) fprintf (stderr, "grdfft: Inverse FFT\n");
  376.  
  377.         ksign = 1;
  378.         fourt_(datac, narray, &ndim, &ksign, &iform, workc);
  379.  
  380.         scale_out *= (2.0 / ndatac);
  381.         for (i = 0; i < ndatac; i+= 2) datac[i] *= scale_out;
  382.  
  383.         /* The data are in the middle of the padded array:  */
  384.  
  385.         if (write_grd (outfile, &h, datac, h.x_min, h.x_max, h.y_min, h.y_max, pad, TRUE)) {
  386.             fprintf (stderr, "grdfft: Error writing file %s\n", outfile);
  387.             exit (-1);
  388.         }
  389.     }
  390.     
  391.     free((char *)workc);
  392.     free((char *)datac);
  393.  
  394.     if (gmtdefs.verbose) fprintf (stderr, "grdfft: Done\n");
  395.  
  396.     gmt_end (argc, argv);
  397.  
  398.     exit(0);
  399. }
  400.  
  401. read_data(argc, argv)
  402. int argc;
  403. char **argv;
  404. {
  405.     int    worksize, factors[32];
  406.     double    tdummy, edummy;
  407.     void    suggest_fft();
  408.     void    fourt_stats();
  409.  
  410.     if (read_grd_info (infile, &h)) {
  411.         fprintf (stderr, "grdfft: Error opening file %s\n", infile);
  412.         return (-1);
  413.     }
  414.     
  415.     grd_init (&h, argc, argv, TRUE);
  416.  
  417.     /* Get dimensions as may be appropriate:  */
  418.     if (n_user_set) {
  419.         if (nx2 < h.nx || ny2 < h.ny) {
  420.             fprintf(stderr,"grdfft: Error: You specified -Nnx/ny smaller than input grdfile.  Ignored.\n");
  421.             n_user_set = FALSE;
  422.         }
  423.     }
  424.     if (!(n_user_set) ) {
  425.         if (force_narray) {
  426.             nx2 = h.nx;
  427.             ny2 = h.ny;
  428.         }
  429.         else {
  430.             suggest_fft(h.nx, h.ny, fft_sug, (gmtdefs.verbose || suggest_narray));
  431.             if (fft_sug[1].totalbytes < fft_sug[0].totalbytes) {
  432.                 /* The most accurate solution needs same or less storage
  433.                  * as the fastest solution; use the most accurate's dimensions:  */
  434.                 nx2 = fft_sug[1].nx;
  435.                 ny2 = fft_sug[1].ny;
  436.             }
  437.             else {
  438.                 /* Use the sizes of the fastest solution  */
  439.                 nx2 = fft_sug[0].nx;
  440.                 ny2 = fft_sug[0].ny;
  441.             }
  442.         }
  443.     }
  444.  
  445.     /* Get here when nx2 and ny2 are set to the vals we will use.  */
  446.     narray[0] = nx2;
  447.     narray[1] = ny2;
  448.     fourt_stats(nx2, ny2, factors, &edummy, &worksize, &tdummy);
  449.     if (gmtdefs.verbose) fprintf(stderr,"grdfft:  Data dimension %d %d\tFFT dimension %d %d\n",
  450.         h.nx, h.ny, nx2, ny2);
  451.  
  452.     /* Make an array of floats 2 * nx2 * ny2 for complex data:  */
  453.     ndatac = 2 * nx2 * ny2;
  454.         datac = (float *) memory (CNULL, ndatac, sizeof(float), "grdfft");
  455.     memset ((char *)datac, 0, ndatac*sizeof(float)); 
  456.     if (worksize) {
  457.                 if (worksize < nx2) worksize = nx2;
  458.                 if (worksize < ny2) worksize = ny2;
  459.         worksize *= 2;
  460.             workc = (float *) memory (CNULL, worksize, sizeof(float), "grdfft");
  461.         memset ((char *)workc, 0, worksize*sizeof(float)); 
  462.     }
  463.     else {
  464.         workc = (float *) memory (CNULL, 4, sizeof(float), "grdfft");
  465.         memset ((char *)workc, 0, 4*sizeof(float)); 
  466.     }
  467.  
  468.     /* Put the data in the middle of the padded array:  */
  469.  
  470.     i_data_start = pad[0] = (nx2 - h.nx)/2;    /* zero if nx2 < h.nx+1  */
  471.     j_data_start = pad[3] = (ny2 - h.ny)/2;
  472.     pad[1] = nx2 - h.nx - pad[0];
  473.     pad[2] = ny2 - h.ny - pad[3];
  474.  
  475.     if (read_grd (infile, &h, datac, h.x_min, h.x_max, h.y_min, h.y_max, pad, TRUE)) {
  476.         fprintf (stderr, "grdfft: Error reading file %s\n", infile);
  477.         return (-1);
  478.     }
  479.  
  480.     if (!(leave_trend_alone) )remove_plane();
  481.     if (!(force_narray) )taper_edges();
  482.  
  483.     delta_kx = 2 * M_PI / (nx2 * h.x_inc);
  484.     delta_ky = 2 * M_PI / (ny2 * h.y_inc);
  485.     if (map_units) {
  486.         /* Give delta_kx, delta_ky units of 2pi/meters  */
  487.         double tmp;
  488.         tmp = 2.0 * M_PI * gmtdefs.ellipse[N_ELLIPSOIDS-1].eq_radius / 360.0;    /* GRS-80 sphere m/degree */
  489.         delta_kx /= (tmp * cos(0.5 * (h.y_min + h.y_max) * D2R) );
  490.         delta_ky /= tmp;
  491.     }
  492.  
  493.     return(0);
  494. }
  495.  
  496. int    remove_plane()
  497. {
  498.     /* Remove the best-fitting plane by least squares.
  499.  
  500.     Let plane be z = a0 + a1 * x + a2 * y.  Choose the
  501.     center of x,y coordinate system at the center of 
  502.     the array.  This will make the Normal equations 
  503.     matrix G'G diagonal, so solution is trivial.  Also,
  504.     spend some multiplications on normalizing the 
  505.     range of x,y into [-1,1], to avoid roundoff error.  */
  506.  
  507.     int    one_or_zero;
  508.     double    x_half_length, one_on_xhl, y_half_length, one_on_yhl;
  509.     double    sumx2, sumy2, x, y, z;
  510.  
  511.     one_or_zero = (h.node_offset) ? 0 : 1;
  512.     x_half_length = 0.5 * (h.nx - one_or_zero);
  513.     one_on_xhl = 1.0 / x_half_length;
  514.     y_half_length = 0.5 * (h.ny - one_or_zero);
  515.     one_on_yhl = 1.0 / y_half_length;
  516.  
  517.     sumx2 = sumy2 = data_var = 0.0;
  518.     a[2] = a[1] = a[0] = 0.0;
  519.  
  520.     for (j = 0; j < h.ny; j++) {
  521.         y = one_on_yhl * (j - y_half_length);
  522.         for (i = 0; i < h.nx; i++) {
  523.             x = one_on_xhl * (i - x_half_length);
  524.             z = datac[ij_data(i,j)];
  525.             a[0] += z;
  526.             a[1] += z*x;
  527.             a[2] += z*y;
  528.             sumx2 += x*x;
  529.             sumy2 += y*y;
  530.         }
  531.     }
  532.     a[0] /= (h.nx*h.ny);
  533.     a[1] /= sumx2;
  534.     a[2] /= sumy2;
  535.     for (j = 0; j < h.ny; j++) {
  536.         y = one_on_yhl * (j - y_half_length);
  537.         for (i = 0; i < h.nx; i++) {
  538.             x = one_on_xhl * (i - x_half_length);
  539.             datac[ij_data(i,j)] -= (a[0] + a[1]*x + a[2]*y);
  540.             data_var += (datac[ij_data(i,j)] * datac[ij_data(i,j)]);
  541.         }
  542.     }
  543.     data_var = sqrt(data_var / (h.nx*h.ny - 1));
  544.     /* Rescale a1,a2 into user's units, in case useful later:  */
  545.     a[1] *= (2.0/(h.x_max - h.x_min));
  546.     a[2] *= (2.0/(h.y_max - h.y_min));
  547.     if(gmtdefs.verbose)fprintf (stderr,"grdfft: Plane removed.  Mean, S.D., Dx, Dy:  %.8lg\t%.8lg\t%.8lg\t%.8lg\n", a[0],data_var,a[1],a[2]);
  548. }
  549.  
  550. int    taper_edges ()
  551. {
  552.     int    im, jm, il1, ir1, il2, ir2, jb1, jb2, jt1, jt2;
  553.     double    scale, cos_wt;
  554.  
  555.     /* Note that if nx2 = h.nx+1 and ny2 = h.ny + 1, then this routine
  556.         will do nothing; thus a single row/column of zeros may be
  557.         added to the bottom/right of the input array and it cannot
  558.         be tapered.  But when (nx2 - h.nx)%2 == 1 or ditto for y,
  559.         this is zero anyway.  */
  560.  
  561.  
  562.     /* First reflect about xmin and xmax, point symmetric about edge point:  */
  563.  
  564.     for (im = 1; im <= i_data_start; im++) {
  565.         il1 = -im;    /* Outside xmin; left of edge 1  */
  566.         ir1 = im;    /* Inside xmin; right of edge 1  */
  567.         il2 = il1 + h.nx - 1;    /* Inside xmax; left of edge 2  */
  568.         ir2 = ir1 + h.nx - 1;    /* Outside xmax; right of edge 2  */
  569.         for (j = 0; j < h.ny; j++) {
  570.             datac[ij_data(il1,j)] = 2.0*datac[ij_data(0,j)] - datac[ij_data(ir1,j)];
  571.             datac[ij_data(ir2,j)] = 2.0*datac[ij_data((h.nx-1),j)] - datac[ij_data(il2,j)];
  572.         }
  573.     }
  574.  
  575.     /* Next, reflect about ymin and ymax.
  576.         At the same time, since x has been reflected,
  577.         we can use these vals and taper on y edges:  */
  578.  
  579.     scale = M_PI / (j_data_start + 1);
  580.  
  581.     for (jm = 1; jm <= j_data_start; jm++) {
  582.         jb1 = -jm;    /* Outside ymin; bottom side of edge 1  */
  583.         jt1 = jm;    /* Inside ymin; top side of edge 1  */
  584.         jb2 = jb1 + h.ny - 1;    /* Inside ymax; bottom side of edge 2  */
  585.         jt2 = jt1 + h.ny - 1;    /* Outside ymax; bottom side of edge 2  */
  586.         cos_wt = 0.5 * (1.0 + cos(jm * scale) );
  587.         for (i = -i_data_start; i < nx2 - i_data_start; i++) {
  588.             datac[ij_data(i,jb1)] = cos_wt * (2.0*datac[ij_data(i,0)] - datac[ij_data(i,jt1)]);
  589.             datac[ij_data(i,jt2)] = cos_wt * (2.0*datac[ij_data(i,(h.ny-1))] - datac[ij_data(i,jb2)]);
  590.         }
  591.     }
  592.  
  593.     /* Now, cos taper the x edges:  */
  594.  
  595.     scale = M_PI / (i_data_start + 1);
  596.     for (im = 1; im <= i_data_start; im++) {
  597.         il1 = -im;
  598.         ir1 = im;
  599.         il2 = il1 + h.nx - 1;
  600.         ir2 = ir1 + h.nx - 1;
  601.         cos_wt = 0.5 * (1.0 + cos(im * scale) );
  602.         for (j = -j_data_start; j < ny2 - j_data_start; j++) {
  603.             datac[ij_data(il1,j)] *= cos_wt;
  604.             datac[ij_data(ir2,j)] *= cos_wt;
  605.         }
  606.     }
  607.  
  608.     /* Now, renormalize the variance:  COMMENTED OUT--STUPID IDEA
  609.  
  610.     data2_var = 0.0;
  611.     for (i = 0; i < ndatac; i+= 2) data2_var += (datac[i] * datac[i]);
  612.  
  613.     data2_var = sqrt(data2_var / (ndatac/2 - 1));
  614.     scale = data_var / data2_var;
  615.     for (i = 0; i < ndatac; i+= 2) datac[i] *= scale;
  616.  
  617.     if (gmtdefs.verbose) fprintf (stderr, "grdfft:  Cosine tapering rescaled data by %.8lg\n", scale);  */
  618.     return(0);
  619. }
  620.  
  621. double    kx(k)
  622. int    k;
  623. {
  624.     /* Return the value of kx given k,
  625.         where kx = 2 pi / lambda x,
  626.         and k refers to the position
  627.         in the datac array, datac[k].  */
  628.  
  629.     int    ii;
  630.  
  631.     ii = ((int)(k/2))%nx2;
  632.     if (ii > nx2/2) ii -= nx2;
  633.     return(ii * delta_kx);
  634. }
  635.  
  636. double    ky(k)
  637. int k;
  638. {
  639.     /* Return the value of ky given k,
  640.         where ky = 2 pi / lambda y,
  641.         and k refers to the position
  642.         in the datac array, datac[k].  */
  643.  
  644.     int    jj;
  645.  
  646.     jj = ((int)(k/2))/nx2;
  647.     if (jj > ny2/2) jj -= ny2;
  648.     return(jj * delta_ky);
  649. }
  650.  
  651. double    modk(k)
  652. int    k;
  653. {
  654.     /* Return the value of sqrt(kx*kx + ky*ky),
  655.         given k, where k is array position.  */
  656.  
  657.     return (hypot (kx(k), ky(k)));
  658. }
  659.  
  660. int    do_differentiate (par)
  661. double *par;
  662. {
  663.     int    k;
  664.     double scale, fact;
  665.  
  666.     /* Differentiate in frequency domain by multiplying by kr [scale optional] */
  667.     
  668.     scale = (*par != 0.0) ? *par : 1.0;
  669.     datac[0] = datac[1] = 0.0;
  670.     for (k = 2; k < ndatac; k+= 2) {
  671.         fact = scale * modk(k);
  672.         datac[k] *= fact;
  673.         datac[k+1] *= fact;
  674.     }
  675. }
  676.  
  677. int    do_integrate (par)
  678. double *par;
  679. {
  680.     /* Integrate in frequency domain by deviding by kr [scale optional] */
  681.     int    k;
  682.     double fact, scale;
  683.  
  684.     scale = (*par != 0.0) ? *par : 1.0;
  685.     datac[0] = datac[1] = 0.0;
  686.     for (k = 2; k < ndatac; k+= 2) {
  687.         fact = 1.0 / (scale * modk(k) );
  688.         datac[k] *= fact;
  689.         datac[k+1] *= fact;
  690.     }
  691. }
  692.  
  693. int    do_continuation (zlevel)
  694. double *zlevel;
  695. {
  696.     int    k;
  697.     double tmp;
  698.  
  699.     /* If z is positive, the field will be upward continued using exp[- k z].  */
  700.  
  701.     for (k = 2; k < ndatac; k+= 2) {
  702.         tmp = exp (-(*zlevel) * modk(k) );
  703.         datac[k] *= tmp;
  704.         datac[k+1] *= tmp;
  705.     }
  706. }
  707.  
  708. int    do_azimuthal_derivative (azim)
  709. double *azim;
  710. {
  711.     int    k;
  712.     float    tempr, tempi, fact;
  713.     double cos_azim, sin_azim;
  714.  
  715.     cos_azim = cos ((*azim) * D2R);
  716.     sin_azim = sin ((*azim) * D2R);
  717.  
  718.     datac[0] = datac[1] = 0.0;
  719.     for (k = 2; k < ndatac; k+= 2) {
  720.         fact = sin_azim * kx(k) + cos_azim * ky(k);
  721.         tempr = -(datac[k+1] * fact);
  722.         tempi = (datac[k] * fact);
  723.         datac[k] = tempr;
  724.         datac[k+1] = tempi;
  725.     }
  726. }
  727.  
  728. int    do_isostasy (par)
  729. double *par;
  730. {
  731.  
  732.     /* Do the isostatic response function convolution in the Freq domain.
  733.     All units assumed to be in SI (that is kx, ky, modk wavenumbers in m**-1,
  734.     densities in kg/m**3, Te in m, etc.
  735.     rw, the water density, is used to set the Airy ratio and the restoring
  736.     force on the plate (rm - ri)*gravity if ri = rw; so use zero for topo in air.  */
  737.     int    k;
  738.     double    airy_ratio, rigidity_d, d_over_restoring_force, mk, k2, k4, transfer_fn;
  739.  
  740.     double    te;    /* Elastic thickness, SI units (m)  */
  741.     double    rl;    /* Load density, SI units  */
  742.     double    rm;    /* Mantle density, SI units  */
  743.     double    rw;    /* Water density, SI units  */
  744.     double    ri;    /* Infill density, SI units  */
  745.     
  746.     /* Kernighan and Ritchie 2nd edition ANSI C says these should be const double, but
  747.         my lint won't recognise that type:  */
  748.         
  749.     double    youngs_modulus = 1.0e11;    /* Pascal = Nt/m**2  */
  750.     double    normal_gravity = 9.80619203;    /* m/s**2  */
  751.     double    poissons_ratio = 0.25;
  752.  
  753.     te = par[0];    rl = par[1];    rm = par[2];    rw = par[3];    ri = par[4];
  754.     rigidity_d = (youngs_modulus * te * te * te) / (12.0 * (1.0 - poissons_ratio * poissons_ratio));
  755.     d_over_restoring_force = rigidity_d / ( (rm - ri) * normal_gravity);
  756.     airy_ratio = -(rl - rw)/(rm - ri);
  757.  
  758.     if (te == 0.0) {    /* Airy isostasy; scale global variable scale_out and return */
  759.         scale_out *= airy_ratio;
  760.         return;
  761.     }
  762.     
  763.     for (k = 0; k < ndatac; k+= 2) {
  764.         mk = modk(k);
  765.         k2 = mk * mk;
  766.         k4 = k2 * k2;
  767.         transfer_fn = airy_ratio / ( (d_over_restoring_force * k4) + 1.0);
  768.         datac[k] *= transfer_fn;
  769.         datac[k+1] *= transfer_fn;
  770.     }
  771. }
  772.  
  773. int    do_filter()
  774. {
  775.     int    k;
  776.     double    weight, f_wt();
  777.     for (k = 0; k < ndatac; k += 2) {
  778.         weight = f_wt(k);
  779.         datac[k] *= weight;
  780.         datac[k+1] *= weight;
  781.     }
  782. }
  783.  
  784. double    f_wt(k)
  785. int    k;
  786. {
  787.     double    freq, return_value = 1.0;
  788.     int    j;
  789.     
  790.     /* Always return immediately when zero weight encountered.  */
  791.     
  792.     for (j = 0; j < 3; j++) {
  793.         if (!(f_info.do_this[j])) continue;
  794.         switch (j) {
  795.             case 0:
  796.                 freq = modk(k);
  797.                 break;
  798.             case 1:
  799.                 freq = kx(k);
  800.                 break;
  801.             case 2:
  802.                 freq = ky(k);
  803.                 break;
  804.         }
  805.         if (freq <= f_info.lc[j] || freq >= f_info.hc[j])
  806.             /* In fully cut range.  Weight is zero.  */
  807.             return(0.0);
  808.         else if (freq > f_info.lc[j] && freq < f_info.lp[j])
  809.             return_value *= 0.5 * (1.0 + cos (M_PI * (freq - f_info.lp[j]) * f_info.ltaper[j]));
  810.         else if (freq > f_info.hp[j] && freq < f_info.hc[j])
  811.             return_value *= 0.5 * (1.0 + cos (M_PI * (freq - f_info.hp[j]) * f_info.htaper[j]));
  812.         /* Else freq is in the fully passed range, so weight is multiplied by 1.0  */
  813.     }
  814.  
  815.     /* Get here when all weights have been multiplied on return_value:  */
  816.     
  817.     return(return_value);
  818. }
  819.  
  820. int    parse_f_string(c)
  821. char    *c;
  822. {
  823.     int    i, j, n_tokens, descending;
  824.     double    fourvals[4];
  825.     char    line[80], *p;
  826.     
  827.     strcpy(line, c);
  828.     i = j = 0;    /* j is Filter type:  r=0, x=1, y=2  */
  829.     
  830.     if (line[i] == 'x') {
  831.         j = 1;
  832.         i++;
  833.     }
  834.     else if (line[i] == 'y') {
  835.         j = 2;
  836.         i++;
  837.     }
  838.     
  839.     f_info.do_this[j] = TRUE;
  840.     
  841.     n_tokens = 0;
  842.     p = strtok(&line[i],"/");
  843.     while (p) {
  844.         if (n_tokens > 3) {
  845.             fprintf(stderr,"Too many slashes.\n");
  846.             return(TRUE);
  847.         }
  848.         if(p[0] == '-') {
  849.             fourvals[n_tokens] = -1.0;
  850.         }
  851.         else {
  852.             if ((sscanf(p, "%lf", &fourvals[n_tokens])) != 1) {
  853.                 fprintf(stderr,"grdfft:  Cannot read token %d.\n", n_tokens);
  854.                 return(TRUE);
  855.             }
  856.         }
  857.         n_tokens++;
  858.         p = strtok(CNULL, "/");
  859.     }
  860.     
  861.     if (n_tokens != 4) {
  862.         fprintf(stderr,"Cannot find 4 tokens separated by slashes.\n");
  863.         return(TRUE);
  864.     }
  865.  
  866.     descending = TRUE;
  867.     for (i = 1; i < 4; i++) {
  868.         if (fourvals[i] == -1.0 || fourvals[i-1] == -1.0) continue;
  869.         if (fourvals[i] > fourvals[i-1]) descending = FALSE;
  870.     }
  871.     if (!(descending)) {
  872.         fprintf(stderr,"Wavelengths are not in descending order.\n");
  873.         return(TRUE);
  874.     }
  875.     if ( (fourvals[0] * fourvals[1]) < 0.0 || (fourvals[2] * fourvals[3]) < 0.0) {
  876.         fprintf(stderr,"Pass/Cut specification error.\n");
  877.         return(TRUE);
  878.     }
  879.     
  880.     /* Now everything is OK  */
  881.     
  882.     if (fourvals[0] >= 0.0 || fourvals[1] >= 0.0) {    /* Lower end values are set  */
  883.         f_info.lc[j] = (2.0 * M_PI)/fourvals[0];
  884.         f_info.lp[j] = (2.0 * M_PI)/fourvals[1];
  885.         if (fourvals[0] != fourvals[1]) f_info.ltaper[j] = 1.0/(f_info.lc[j] - f_info.lp[j]);
  886.     }
  887.     
  888.     if (fourvals[2] >= 0.0 || fourvals[3] >= 0.0) {    /* Higher end values are set  */
  889.         f_info.hp[j] = (2.0 * M_PI)/fourvals[2];
  890.         f_info.hc[j] = (2.0 * M_PI)/fourvals[3];
  891.         if (fourvals[2] != fourvals[3]) f_info.htaper[j] = 1.0/(f_info.hc[j] - f_info.hp[j]);
  892.     }
  893.     
  894.     return(FALSE);
  895. }
  896.  
  897. int    do_spectrum(par, give_wavelength)
  898. double *par;
  899. int    give_wavelength;
  900. {
  901.     /* This is modeled on the 1-D case, using the following ideas:
  902.      *    In 1-D, we ensemble average over samples of length L = 
  903.      *    n * dt.  This gives n/2 power spectral estimates, at
  904.      *    frequencies i/L, where i = 1, n/2.  If we have a total
  905.      *    data set of ndata, we can make nest=ndata/n independent
  906.      *    estimates of the power.  Our standard error is then
  907.      *    1/sqrt(nest).
  908.      *    In making 1-D estimates from 2-D data, we may choose
  909.      *    n and L from nx2 or ny2 and delta_kx, delta_ky as 
  910.      *    appropriate.  In this routine, we are giving the sum over
  911.      *     all frequencies in the other dimension; that is, an
  912.      *    approximation of the integral.
  913.      */
  914.  
  915.     int    k, nk, nused, ifreq;
  916.     double    delta_k, r_delta_k, freq, *power, eps_pow;
  917.     PFD    get_k;
  918.     char    format[64];
  919.  
  920.     if (*par > 0.0) {
  921.         /* X spectrum desired  */
  922.         delta_k = delta_kx;
  923.         nk = nx2/2;
  924.         get_k = (PFD)kx;
  925.     }
  926.     else if (*par < 0.0) {
  927.         /* Y spectrum desired  */
  928.         delta_k = delta_ky;
  929.         nk = ny2/2;
  930.         get_k = (PFD)ky;
  931.     }
  932.     else {
  933.         /* R spectrum desired  */
  934.         if (delta_kx < delta_ky) {
  935.             delta_k = delta_kx;
  936.             nk = nx2/2;
  937.         }
  938.         else {
  939.             delta_k = delta_ky;
  940.             nk = ny2/2;
  941.         }
  942.         get_k = (PFD)modk;
  943.     }
  944.  
  945.     /* Get an array for summing stuff:  */
  946.     power = (double *) memory(CNULL, nk, sizeof(double), "grdfft");
  947.     for (k = 0; k < nk; k++) power[k] = 0.0;
  948.  
  949.     /* Loop over it all, summing and storing, checking range for r:  */
  950.  
  951.     r_delta_k = 1.0 / delta_k;
  952.     
  953.     for (nused = 0, k = 2; k < ndatac; k+= 2) {
  954.         freq = (*get_k)(k);
  955.         ifreq = irint(fabs(freq)*r_delta_k) - 1;
  956.         if (ifreq < 0) ifreq = 0;    /* Might happen when doing r spectrum  */
  957.         if (ifreq >= nk) continue;    /* Might happen when doing r spectrum  */
  958.         power[ifreq] += (datac[k]*datac[k] + datac[k+1]*datac[k+1]);
  959.         nused++;
  960.     }
  961.  
  962.     /* Now get here when array is summed.  */
  963.     eps_pow = 1.0/sqrt((double)nused/(double)nk);
  964.     delta_k /= (2.0*M_PI);    /* Write out frequency, not wavenumber  */
  965.     sprintf (format, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  966.     for (k = 0; k < nk; k++) {
  967.         freq = (k + 1) * delta_k;
  968.         if (give_wavelength) freq = 1.0/freq;
  969.         printf(format, freq, power[k], eps_pow * power[k]);
  970.     }
  971.     return;
  972. }
  973.  
  974. void    fourt_stats(nx, ny, f, r, s, t)
  975. int    nx, ny, f[], *s;
  976. double    *r, *t;
  977. {
  978.     /* Find the proportional run time, t, and rms relative error, r,
  979.      * of a Fourier transform of size nx,ny.  Also gives s, the size
  980.      * of the workspace that will be needed by the transform.
  981.      * To use this routine for a 1-D transform, set ny = 1.
  982.      * 
  983.      * This is all based on the comments in Norman Brenner's code
  984.      * FOURT, from which our C codes are translated.
  985.      * Brenner says:
  986.      * r = 3 * pow(2, -FSIGNIF) * sum{ pow(prime_factors, 1.5) }
  987.      * where FSIGNIF is the smallest bit in the floating point fraction.
  988.      * 
  989.      * Let m = largest prime factor in the list of factors.
  990.      * Let p = product of all primes which appear an odd number of
  991.      * times in the list of prime factors.  Then the worksize needed
  992.      * s = max(m,p).  However, we know that if n is radix 2, then no
  993.      * work is required; yet this formula would say we need a worksize
  994.      * of at least 2.  So I will return s = 0 when max(m,p) = 2.
  995.      *
  996.      * I have two different versions of the comments in FOURT, with
  997.      * different formulae for t.  The simple formula says 
  998.      *     t = n * (sum of prime factors of n).
  999.      * The more complicated formula gives coefficients in microsecs
  1000.      * on a cdc3300 (ancient history, but perhaps proportional):
  1001.      *    t = 3000 + n*(500 + 43*s2 + 68*sf + 320*nf),
  1002.      * where s2 is the sum of all factors of 2, sf is the sum of all
  1003.      * factors greater than 2, and nf is the number of factors != 2.
  1004.      * We know that factors of 2 are very good to have, and indeed,
  1005.      * Brenner's code calls different routines depending on whether
  1006.      * the transform is of size 2 or not, so I think that the second
  1007.      * formula is more correct, using proportions of 43:68 for 2 and
  1008.      * non-2 factors.  So I will use the more complicated formula.
  1009.      * However, I realize that the actual numbers are wrong for today's
  1010.      * architectures, and the relative proportions may be wrong as well.
  1011.      * 
  1012.      * W. H. F. Smith, 26 February 1992.
  1013.      *  */
  1014.  
  1015.     int    n_factors, i, ntotal, sum2, sumnot2, nnot2;
  1016.     int    nonsymx, nonsymy, nonsym, storage;
  1017.     double    err_scale, sig_bits = FSIGNIF;
  1018.  
  1019.  
  1020.     /* Find workspace needed.  First find non_symmetric factors in nx, ny  */
  1021.     n_factors = get_prime_factors(nx, f);
  1022.     nonsymx = get_non_symmetric_f(f, n_factors);
  1023.     n_factors = get_prime_factors(ny, f);
  1024.     nonsymy = get_non_symmetric_f(f, n_factors);
  1025.     nonsym = MAX(nonsymx, nonsymy);
  1026.  
  1027.     /* Now get factors of ntotal  */
  1028.     ntotal = nx * ny;
  1029.         n_factors = get_prime_factors(ntotal, f);
  1030.     storage = MAX(nonsym, f[n_factors - 1]);
  1031.     if (storage == 2)
  1032.         *s = 0;
  1033.     else
  1034.         *s = storage;
  1035.  
  1036.     /* Now find time and error estimates:  */
  1037.  
  1038.     err_scale = 0.0;
  1039.     sum2 = 0;
  1040.     sumnot2 = 0;
  1041.     nnot2 = 0;
  1042.     for(i = 0; i < n_factors; i++) {
  1043.         if (f[i] == 2)
  1044.             sum2 += f[i];
  1045.         else {
  1046.             sumnot2 += f[i];
  1047.             nnot2++;
  1048.         }
  1049.         err_scale += pow((double)f[i], 1.5);
  1050.     }
  1051.     *t = 1.0e-06*(3000.0 + ntotal * (500.0 + 43.0*sum2 + 68.0*sumnot2 + 320.0*nnot2));
  1052.     *r = err_scale * 3.0 * pow(2.0, -sig_bits);
  1053.     return;
  1054. } 
  1055.  
  1056. int    get_non_symmetric_f(f, n)
  1057. int    f[], n;
  1058. {
  1059.     /* Return the product of the non-symmetric factors in f[]  */
  1060.     int    i = 0, j = 1, retval = 1;
  1061.  
  1062.     if (n == 1) return(f[0]);
  1063.  
  1064.     while(i < n) {
  1065.         while(j < n && f[j] == f[i]) j++;
  1066.         if ((j-i)%2) retval *= f[i];
  1067.         i = j;
  1068.         j = i + 1;
  1069.     }
  1070.     if (retval == 1) retval = 0;    /* There are no non-sym factors  */
  1071.     return(retval);
  1072. }
  1073.  
  1074. void    suggest_fft(nx, ny, fft_sug, do_print)
  1075. int    nx, ny, do_print;
  1076. struct FFT_SUGGESTION fft_sug[];
  1077. {
  1078.     int    f[64], xstop, ystop;
  1079.     int    nx_best_t, ny_best_t;
  1080.     int    nx_best_e, ny_best_e;
  1081.     int    nx_best_s, ny_best_s;
  1082.         int     nxg, nyg;       /* Guessed by this routine  */
  1083.         int     nx2, ny2, nx3, ny3, nx5, ny5;   /* For powers  */
  1084.     double    current_time, best_time, given_time, s_time, e_time;
  1085.     int    current_space, best_space, given_space, e_space, t_space;
  1086.     double    current_err, best_err, given_err, s_err, t_err;
  1087.     void    fourt_stats();
  1088.  
  1089.  
  1090.     fourt_stats(nx, ny, f, &given_err, &given_space, &given_time);
  1091.     given_space += nx*ny;
  1092.     given_space *= 8;
  1093.     if (do_print) fprintf(stderr,"grdfft:  Data dimension\t%d %d\ttime factor %.8lg\trms error %.8le\tbytes %d\n",
  1094.         nx, ny, given_time, given_err, given_space);
  1095.  
  1096.     best_err = s_err = t_err = given_err;
  1097.     best_time = s_time = e_time = given_time;
  1098.     best_space = t_space = e_space = given_space;
  1099.     nx_best_e = nx_best_t = nx_best_s = nx;
  1100.     ny_best_e = ny_best_t = ny_best_s = ny;
  1101.  
  1102.     xstop = 2 * nx;
  1103.     ystop = 2 * ny;
  1104.  
  1105.         for (nx2 = 2; nx2 <= xstop; nx2 *= 2) {
  1106.           for (nx3 = 1; nx3 <= xstop; nx3 *= 3) {
  1107.             for (nx5 = 1; nx5 <= xstop; nx5 *= 5) {
  1108.                 nxg = nx2 * nx3 * nx5;
  1109.                 if (nxg < nx || nxg > xstop) continue;
  1110.  
  1111.                 for (ny2 = 2; ny2 <= ystop; ny2 *= 2) {
  1112.                   for (ny3 = 1; ny3 <= ystop; ny3 *= 3) {
  1113.                     for (ny5 = 1; ny5 <= ystop; ny5 *= 5) {
  1114.                         nyg = ny2 * ny3 * ny5;
  1115.                         if (nyg < ny || nyg > ystop) continue;
  1116.  
  1117.             fourt_stats(nxg, nyg, f, ¤t_err, ¤t_space, ¤t_time);
  1118.             current_space += nxg*nyg;
  1119.             current_space *= 8;
  1120.             if (current_err < best_err) {
  1121.                 best_err = current_err;
  1122.                 nx_best_e = nxg;
  1123.                 ny_best_e = nyg;
  1124.                 e_time = current_time;
  1125.                 e_space = current_space;
  1126.             }
  1127.             if (current_time < best_time) {
  1128.                 best_time = current_time;
  1129.                 nx_best_t = nxg;
  1130.                 ny_best_t = nyg;
  1131.                 t_err = current_err;
  1132.                 t_space = current_space;
  1133.             }
  1134.             if (current_space < best_space) {
  1135.                 best_space = current_space;
  1136.                 nx_best_s = nxg;
  1137.                 ny_best_s = nyg;
  1138.                 s_time = current_time;
  1139.                 s_err = current_err;
  1140.             }
  1141.  
  1142.             }
  1143.           }
  1144.         }
  1145.  
  1146.         }
  1147.       }
  1148.     }
  1149.  
  1150.     if (do_print) {
  1151.         fprintf(stderr,"grdfft:  Highest speed\t%d %d\ttime factor %.8lg\trms error %.8le\tbytes %d\n",
  1152.             nx_best_t, ny_best_t, best_time, t_err, t_space);
  1153.         fprintf(stderr,"grdfft:  Most accurate\t%d %d\ttime factor %.8lg\trms error %.8le\tbytes %d\n",
  1154.             nx_best_e, ny_best_e, e_time, best_err, e_space);
  1155.         fprintf(stderr,"grdfft:  Least storage\t%d %d\ttime factor %.8lg\trms error %.8le\tbytes %d\n",
  1156.             nx_best_s, ny_best_s, s_time, s_err, best_space);
  1157.     }
  1158.     /* Fastest solution:  */
  1159.     fft_sug[0].nx = nx_best_t;
  1160.     fft_sug[0].ny = ny_best_t;
  1161.     fft_sug[0].worksize = (t_space/8) - (nx_best_t * ny_best_t);
  1162.     fft_sug[0].totalbytes = t_space;
  1163.     fft_sug[0].run_time = best_time;
  1164.     fft_sug[0].rms_rel_err = t_err;
  1165.     /* Most accurate solution:  */
  1166.     fft_sug[1].nx = nx_best_e;
  1167.     fft_sug[1].ny = ny_best_e;
  1168.     fft_sug[1].worksize = (e_space/8) - (nx_best_e * ny_best_e);
  1169.     fft_sug[1].totalbytes = e_space;
  1170.     fft_sug[1].run_time = e_time;
  1171.     fft_sug[1].rms_rel_err = best_err;
  1172.     /* Least storage solution:  */
  1173.     fft_sug[2].nx = nx_best_s;
  1174.     fft_sug[2].ny = ny_best_s;
  1175.     fft_sug[2].worksize = (best_space/8) - (nx_best_s * ny_best_s);
  1176.     fft_sug[2].totalbytes = best_space;
  1177.     fft_sug[2].run_time = s_time;
  1178.     fft_sug[2].rms_rel_err = s_err;
  1179.  
  1180.     return;
  1181. }
  1182.  
  1183. int    get_prime_factors(n, f)
  1184. int    n, f[];
  1185. {
  1186.     /* Fills the integer array f with the prime factors of n.
  1187.      * Returns the number of locations filled in f, which is
  1188.      * one if n is prime.
  1189.      *
  1190.      * f[] should have been malloc'ed to enough space before
  1191.      * calling prime_factors().  We can be certain that f[32]
  1192.      * is enough space, for if n fits in a long, then n < 2**32,
  1193.      * and so it must have fewer than 32 prime factors.  I think
  1194.      * that in general, ceil(log2((double)n)) is enough storage
  1195.      * space for f[].
  1196.      *
  1197.      * Tries 2,3,5 explicitly; then alternately adds 2 or 4
  1198.      * to the previously tried factor to obtain the next trial
  1199.      * factor.  This is done with the variable two_four_toggle.
  1200.      * With this method we try 7,11,13,17,19,23,25,29,31,35,...
  1201.      * up to a maximum of sqrt(n).  This shortened list results
  1202.      * in 1/3 fewer divisions than if we simply tried all integers
  1203.      * between 5 and sqrt(n).  We can reduce the size of the list
  1204.      * of trials by an additional 20% by removing the multiples
  1205.      * of 5, which are equal to 30m +/- 5, where m >= 1.  Starting
  1206.      * from 25, these are found by alternately adding 10 or 20.
  1207.      * To do this, we use the variable ten_twenty_toggle.
  1208.      *
  1209.      * W. H. F. Smith, 26 Feb 1992, after D.E. Knuth, vol. II  */
  1210.  
  1211.     int    current_factor;    /* The factor currently being tried  */
  1212.     int    max_factor;    /* Don't try any factors bigger than this  */
  1213.     int    n_factors = 0;    /* Returned; one if n is prime  */
  1214.     int    two_four_toggle = 0;    /* Used to add 2 or 4 to get next trial factor  */
  1215.     int    ten_twenty_toggle = 0;    /* Used to add 10 or 20 to skip_five  */
  1216.     int    skip_five = 25;    /* Used to skip multiples of 5 in the list  */
  1217.     int    m;    /* Used to keep a working copy of n  */
  1218.  
  1219.  
  1220.     /* Initialize m and max_factor  */
  1221.     m = abs(n);
  1222.     if (m < 2) return(0);
  1223.     max_factor = floor(sqrt((double)m));
  1224.  
  1225.     /* First find the 2s  */
  1226.     current_factor = 2;
  1227.     while(!(m%current_factor)) {
  1228.         m /= current_factor;
  1229.         f[n_factors] = current_factor;
  1230.         n_factors++;
  1231.     }
  1232.     if (m == 1) return(n_factors);
  1233.  
  1234.     /* Next find the 3s  */
  1235.     current_factor = 3;
  1236.     while(!(m%current_factor)) {
  1237.         m /= current_factor;
  1238.         f[n_factors] = current_factor;
  1239.         n_factors++;
  1240.     }
  1241.     if (m == 1) return(n_factors);
  1242.  
  1243.     /* Next find the 5s  */
  1244.     current_factor = 5;
  1245.     while(!(m%current_factor)) {
  1246.         m /= current_factor;
  1247.         f[n_factors] = current_factor;
  1248.         n_factors++;
  1249.     }
  1250.     if (m == 1) return(n_factors);
  1251.  
  1252.     /* Now try all the rest  */
  1253.  
  1254.     while (m > 1 && current_factor <= max_factor) {
  1255.  
  1256.         /* Current factor is either 2 or 4 more than previous value  */
  1257.  
  1258.         if (two_four_toggle) {
  1259.             current_factor += 4;
  1260.             two_four_toggle = 0;
  1261.         }
  1262.         else {
  1263.             current_factor += 2;
  1264.             two_four_toggle = 1;
  1265.         }
  1266.  
  1267.         /* If current factor is a multiple of 5, skip it.  But first,
  1268.             set next value of skip_five according to 10/20 toggle:  */
  1269.  
  1270.         if (current_factor == skip_five) {
  1271.             if (ten_twenty_toggle) {
  1272.                 skip_five += 20;
  1273.                 ten_twenty_toggle = 0;
  1274.             }
  1275.             else {
  1276.                 skip_five += 10;
  1277.                 ten_twenty_toggle = 1;
  1278.             }
  1279.             continue;
  1280.         }
  1281.  
  1282.         /* Get here when current_factor is not a multiple of 2,3 or 5:  */
  1283.  
  1284.         while(!(m%current_factor)) {
  1285.             m /= current_factor;
  1286.             f[n_factors] = current_factor;
  1287.             n_factors++;
  1288.         }
  1289.     }
  1290.  
  1291.     /* Get here when all factors up to floor(sqrt(n)) have been tried.  */
  1292.  
  1293.     if (m > 1) {
  1294.         /* m is an additional prime factor of n  */
  1295.         f[n_factors] = m;
  1296.         n_factors++;
  1297.     }
  1298.     return (n_factors);
  1299. }
  1300.  
  1301.