rings/apps/pspec3.c

/*
 *  pspec3.c                                            ~rick/src/pspec
 *
 *  Calculate the 3-dimensional power spectrum of an input dataset
 *  The module expects a 3-dimensional real dataset
 *
 *  Responsible:
 *      Rick Bogart                             RBogart@solar.stanford.edu
 *
 *  Usage:
 *    pspec3 [-lvx] in= pspec= ...
 *
 *  Arguments: (type            default         description)
 *      in      DataSet none            Input dataset
 *                      A series of records containing 3-d data cubes as one
 *                      of their segments. It is assumed that all records in
 *                      the dataset are from the same dataseries, or at least
 *                      share a common segment structure
 *      segment string  -               Name of the input data segment (only
 *                      required if the input series has multiple 3-dimensional
 *                      segments)
 *      pspec   DataSeries      -               Output data series;  must
 *                      contain at least one 3-dimensional segment
 *    mask_in   Float           0.9375          Inner apodization radius
 *    mask_ex   Float           1.0             Outer apodization radius
 *    apodize   Float           0.96875         Temporal apodization edge
 *    fbin      int             0               output frequency binning
 *
 *  Flags
 *      -l      output direct tpower spectrum rather than scaled log
 *      -n      do not save results (diagnostic only)
 *      -v      run verbose
 *      -x      use double-precision calculation
 *
 *  Bugs:
 *    For large input data cubes (~1 Gpxl) the results may be garbage,
 *      for larger ones (~2 Gpxl) the input may fail altogether. The exact
 *      causes of these problems and consequently the exact limits of
 *      validity are not known, but the code has been verified to work with
 *      "normal" input spectra (Dopplergrams as compressed scaled shorts)
 *      of up to 800 Mpxl.
 *    No checking that input data series corresponds with link in output
 *      series
 *    The value of CDELT3 in the input set overrides the value of Cadence
 *      if propagated, but is not checked for units
 *
 *  Future Updates
 *    Use drms_keyword_lookup()
 *    Create new data series as needed
 *    Parallelize over multiple records
 *    Add recreation option for target output records
 *    Add parameters for overriding essential input keywords
 *    Add parameter for overriding keyword list for propagation
 *    Modify to also work with a series of two-dimensional images - maybe
 *      not such a good idea
 *
 *  Revision history is at end of file.
 */
#include <jsoc_main.h>
#include <fftw3.h>

char *module_name = "pspec3";
char *module_desc = "3-d power spectrum";
char *version_id = "1.1";

ModuleArgs_t module_args[] = {
  {ARG_DATASET, "in", "", "Input data set"},
  {ARG_STRING,  "segment", "Not Specified",
      "input data series segment; ignored if series only has one 3-d segment"}, 
  {ARG_DATASERIES, "pspec", "", "Ouput data series"},
  {ARG_FLOAT, "mask_in", "0.9375", "inner radial apodization edge", "[0.0,)"},
  {ARG_FLOAT, "mask_ex", "1.0", "outer radial apodization edge", "[0.0,)"},
  {ARG_FLOAT,   "apodize", "0.96875", "temporal apodization edge",
        "[0.0,1.0]"},
  {ARG_INT,     "fbin", "0", "output frequency fbins (0-> none)"},
  {ARG_STRING, "copy",  "+",
      "comma separated list of keys to propagate forward"},
  {ARG_FLAG,    "l",    "",
        "output direct power spectrum rather than scaled log"},
  {ARG_FLAG,    "n", "0", "do not save output record (diagnostics only)"},      
  {ARG_FLAG,    "x",    "", "use double-precision calculation"},
  {ARG_FLAG,    "v",    "", "verbose mode"}, 
  {}
};
       /*  list of keywords to propagate (if possible) from input to output  */
char *propagate[] = {"CarrTime", "CarrRot", "CMLon", "LonHG", "LatHG", "LonCM",
    "MidTime", "Duration", "LonSpan", "T_START", "T_STOP", "Coverage",
    "ZonalTrk", "ZonalVel", "MeridTrk", "MeridVel", "MapScale", "MAI", "Ident",
    "Width", "Height", "Size", "MapProj", "Map_PA", "PosAng", "RSunRef"};

#include "keystuff.c"

static int cleanup (int error, DRMS_RecordSet_t *irecs, DRMS_Record_t *orec,
    DRMS_Array_t *orig, DRMS_Array_t *pspec, int dispose) {
  if (orig) drms_free_array (orig);
  if (pspec) drms_free_array (pspec);
  if (irecs) drms_close_records (irecs, DRMS_FREE_RECORD);
  if (orec) {
    if (error) drms_close_record (orec, DRMS_FREE_RECORD);
    else drms_close_record (orec, dispose);
  }
  return error;
}

int read_from_big_cube (DRMS_Segment_t *seg, int dp_calc, void *rdata,
    double *avg, double *var, double apode_edge, double *apodization,
    int xcols) {
  DRMS_Array_t *tmp;
  DRMS_Type_t type;
  void *data;
  double *wdata;
  double dval, dataavg, datavar, weight;
  float *xdata;
  float fval;
  long long n;
  int *start, *stop;
  int m, col, cols, row, rows, area, plane, planes, rank, size;
  int status;

  if (dp_calc) {
    size = sizeof (double);
    type = DRMS_TYPE_DOUBLE;
    wdata = (double *)rdata;
  } else {
    size = sizeof (float);
    type = DRMS_TYPE_FLOAT;
    xdata = (float *)rdata;
  }
  dataavg = datavar = 0.0;

  rank = seg->info->naxis;
  cols = seg->axis[0];
  rows = seg->axis[1];
  planes = seg->axis[rank-1];
  start = (int *)malloc (rank * sizeof (int));
  stop = (int *)malloc (rank * sizeof (int));
  area = 1;
  for (n = 0; n < rank-1; n++) {
    start[n] = 0;
    stop[n] = seg->axis[n] - 1;
    area *= seg->axis[n];
  }
  size *= area;
  n = 0;
  for (plane = 0; plane < planes; plane++) {
    m = 0;
    start[rank-1] = stop[rank-1] = plane;
    tmp = drms_segment_readslice (seg, type, start, stop, &status);
    if (status) {
      fprintf (stderr, "Error reading data plane %d\n", plane);
      return 1;
    }
    data = tmp->data;
    if (apode_edge < 1.0) {
      double t = (2.0 * plane) / (planes - 1.0);
      if (t > 1.0) t = 2.0 - t;
      t = 1.0 - t;
      if (t > apode_edge) {
        t = (t - apode_edge) / (1.0 - apode_edge);
        t *= t;
        weight = (1.0 - t);
        weight *= weight;
      } else weight = 1.0;
    } else weight = 1.0;
    if (dp_calc) {
      for (row = 0; row < rows; row++) {
        for (col = 0; col < cols; col++) {
          dval = ((double *)data)[m];
          if (isnan (dval)) dval = 0.0;
          dval *= weight * apodization[m++];
          dataavg += dval;
          datavar += dval * dval;
          wdata[n++] = dval;
        }
        n += xcols - cols;
      }
    } else {
      for (row = 0; row < rows; row++) {
        for (col = 0; col < cols; col++) {
          fval = ((float *)data)[m];
          if (isnan (fval)) fval = 0.0;
          fval *= weight * apodization[m++];
          dataavg += fval;
          datavar += fval * fval;
          xdata[n++] = fval;
        }
        n += xcols - cols;
      }
    }
    drms_free_array (tmp);
  }
  return 0;
}

DRMS_Array_t *read_big_cube (DRMS_Segment_t *seg, DRMS_Type_t type, int *status) {
  DRMS_Array_t *arr, *tmp;
  int *start, *stop;
  int n, plane, pln, plen, rank, size;
  void *data;

  size = drms_sizeof (type);
  rank = seg->info->naxis;
  plen = seg->axis[rank-1];
  arr = drms_array_create (type, rank, seg->axis, NULL, status);
  if (*status) return NULL;

  start = (int *)malloc (rank * sizeof (int));
  stop = (int *)malloc (rank * sizeof (int));
  plane = 1;
  for (n = 0; n < rank-1; n++) {
    start[n] = 0;
    stop[n] = seg->axis[n] - 1;
    plane *= seg->axis[n];
  }
  size *= plane;
  data = arr->data;
  for (pln = 0; pln < plen; pln++) {
    start[rank-1] = stop[rank-1] = pln;
    tmp = drms_segment_readslice (seg, type, start, stop, status);
    if (*status) return NULL;
    memcpy (data, tmp->data, size);
    drms_free_array (tmp);
    (char *)data += size;
  }
  if (*status) return NULL;
  return arr;
}

int DoIt (void) {
  CmdParams_t *params = &cmdparams;
  DRMS_RecordSet_t *irecs = NULL;
  DRMS_Record_t *irec, *orec = NULL;
  DRMS_Segment_t *iseg, *oseg;
  DRMS_Array_t *orig = NULL, *pspec = NULL;

  double *apodization;
  double *wdata;
  double crpix, crval = 0.0;
  double dnu, domega, dval, r, r2, ra, ra2, rx, ry, r0x, r0y, t, tstep;
  double dataavg, datavar, powrint, normal, weight;
  double bzero, bscale, scale_range;
  float *data, *xdata;
  float ftrb, ftib, fval, vmin, vmax, vmin2;
  long long cube, ntot, l, m, n;
  int axes[3];
  int rank, col, cols, row, rows, plane, planes, area;
  int hcols, hrows, hplanes, opln, xcols, cols_even, rows_even;
  int rgn, rgnct, segs, isegnum, osegnum, found;
  int bin, is, js;
  int key_n, kstat, propct, status;
  int bigcube;
  char **copykeylist;
  char *inser;
  char module_ident[64];
  char pathname[2*DRMS_MAXPATHLEN+1];
  char source[DRMS_MAXQUERYLEN], recid[DRMS_MAXQUERYLEN];
  void *val;

  int keyct = sizeof (propagate) / sizeof (char *);

  char *inds = strdup (params_get_str (params, "in"));
  char *out_series = strdup (params_get_str (params, "pspec"));
  double apode_edge = params_get_double (params, "apodize");
  double edge_inner = params_get_double (params, "mask_in");
  double edge_outer = params_get_double (params, "mask_ex");
  int fbins = params_get_int (params, "fbin");
  char *seg_name = strdup (params_get_str (params, "segment"));
  char *propagate_req = strdup (params_get_str (params, "copy"));
  int log_out = params_isflagset (params, "l") ? 0 : 1;
  int dp_calc = params_isflagset (params, "x");
  int verbose = params_isflagset (params, "v");
  int no_save = params_isflagset (params, "n");
  int dispose = (no_save) ? DRMS_FREE_RECORD : DRMS_INSERT_RECORD;
                                        /* Declaration for FFTW library call */
  fftwf_plan fplan;
  fftw_plan plan;

  snprintf (module_ident, 64, "%s v %s", module_name, version_id);
  if (verbose) {
    printf ("%s:\n", module_ident);
    if (fbins > 1)
      printf ("  output binned by %d frequencies per plane\n", fbins);
    if (no_save)
      printf ("(diagnostic run only, no records will be written to DRMS)\n");
  }
                                                /*  check the output series */
  orec = drms_create_record (drms_env, out_series, DRMS_TRANSIENT, &status);
  if (status) {
    fprintf (stderr,
        "Error: drms_create_record returned %d for data series:\n", status);
    fprintf (stderr, "       %s\n", out_series);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  if ((segs = orec->segments.num_total) < 1) {
    fprintf (stderr, "Error: no data segments in output data series:\n");
    fprintf (stderr, "  %s\n", out_series);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  found = 0;
  for (n = 0; n < segs; n++) {
    oseg = drms_segment_lookupnum (orec, n);
    if (oseg->info->naxis != 3) continue;
    if (!found) osegnum = n;
    found++;
  }
  if (!found) {
    fprintf (stderr, "Error: no segment of rank 3 in output data series:\n");
    fprintf (stderr, "  %s\n", out_series);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  oseg = drms_segment_lookupnum (orec, osegnum);
  if (found > 1) {
    fprintf (stderr,
        "Warning: multiple segments of rank 3 in output data series:\n");
    fprintf (stderr, "  %s\n", out_series);
    fprintf (stderr, "  using %s\n", oseg->info->name);
  }
  switch (oseg->info->type) {
    case DRMS_TYPE_CHAR:
      scale_range = 250.0;
      break;
    case DRMS_TYPE_SHORT:
      scale_range = 65000.0;
      break;
    case DRMS_TYPE_INT:
      scale_range = 4.2e9;
      break;
    case DRMS_TYPE_LONGLONG:
      scale_range = 1.8e19;
      break;
    default:
      scale_range = 1.0;
  }
  drms_close_record (orec, DRMS_FREE_RECORD);
                                                           /*  check input  */
  irecs = drms_open_records (drms_env, inds, &status);
  if (status) {
    fprintf (stderr, "Error (%s): drms_open_records() returned %d for dataset:\n",
        module_ident, status);
    fprintf (stderr, "  %s\n", inds);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }

  rgnct = irecs->n;
  if (rgnct < 1) {
    fprintf (stderr, "No records found in input data set %s\n", inds);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  irec = irecs->records[0];
  inser = strdup (inds);
  if ((segs = drms_record_numsegments (irec)) < 1) {
    fprintf (stderr, "Error: no data segments in input data series:\n");
    fprintf (stderr, "  %s\n", inser);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  found = 0;
  for (n = 0; n < segs; n++) {
    iseg = drms_segment_lookupnum (irec, n);
    if (iseg->info->naxis != 3) continue;
    if (!found) isegnum = n;
    found++;
  }
  if (!found) {
    fprintf (stderr, "Error: no segment of rank 3 in input data series:\n");
    fprintf (stderr, "  %s\n", inser);
    return cleanup (1, irecs, orec, orig, pspec, dispose);
  }
  if (found > 1) {
    if (strcmp (seg_name, "Not Specified")) {
      iseg = drms_segment_lookup (irec, seg_name);
      if (!iseg) {
        fprintf (stderr,
            "Warning: requested segment %s not found in input data series:\n",
            seg_name);
        fprintf (stderr, "  %s\n", inser);
        iseg = drms_segment_lookupnum (irec, isegnum);
        fprintf (stderr, "  using segement %s\n", iseg->info->name);
      } else if (iseg->info->naxis != 3) {
        fprintf (stderr,
            "Warning: requested segment %s in input data series:\n", seg_name);
        fprintf (stderr, "  %s is not 3-dimensional", inser);
        iseg = drms_segment_lookupnum (irec, isegnum);
        fprintf (stderr, " using segment %s\n", iseg->info->name);
      } else isegnum = iseg->info->segnum;
    } else {
      fprintf (stderr,
          "Warning: multiple segments of rank 3 in input data series:\n");
      fprintf (stderr, "  %s\n", inser);
      fprintf (stderr, "  using %s\n", iseg->info->name);
    }
  }
  propct = construct_stringlist (propagate_req, ',', &copykeylist);
  if (verbose) {
    printf ("propagating %d key(s):\n", propct);
    for (key_n = 0; key_n < propct; key_n++) printf ("  %s\n",
        copykeylist[key_n]);
    printf ("\nprocessing %d record(s) in series %s:\n", rgnct, inser);
  }

                                                       /*  process records  */
  for (rgn = 0; rgn < rgnct; rgn++) {
    irec = irecs->records[rgn];
    drms_sprint_rec_query (source, irec);

    if (!(iseg = drms_segment_lookupnum (irec, isegnum))) {
      fprintf (stderr, "Warning: could not find segment \"V\" or \"Vtrack\"\n");
      fprintf (stderr, "       in %s; skipped\n", source);
      continue;
    }
    if (iseg->info->naxis != 3) {
      fprintf (stderr, "Error: rank of data cube (%d) != 3\n", iseg->info->naxis);
      return cleanup (1, irecs, orec, orig, pspec, dispose);
    }
    if (drms_wcs_timestep (irec, 3, &tstep)) {
      dval = (double)drms_getkey_time (irec, "T_STOP", &status);
      tstep = dval - (double)drms_getkey_time (irec, "T_START", &status);
      if (tstep <= 0.0) {
        fprintf (stderr,
            "Error: insufficient key information for frequency determination\n");
        fprintf (stderr,
            "       in record %s; skipped\n", source);
        continue;
      }
      tstep /= planes - 1.0;
    }
                           /*  the input looks okay, open the output record  */
    orec = drms_create_record (drms_env, out_series, DRMS_PERMANENT, &status);
    if (status) {
      fprintf (stderr,
          "Error: drms_create_record returned %d for data series:\n", status);
      fprintf (stderr, "       %s\n", out_series);
      return cleanup (1, irecs, orec, orig, pspec, dispose);
    }
    cols = iseg->axis[0];
    rows = iseg->axis[1];
    planes = iseg->axis[2];
    area = cols * rows;
    cube = planes * area;
    xcols = 2 * ((cols + 2) / 2);
    hcols = cols / 2;
    hrows = rows / 2;
    cols_even = (cols % 2) ? 0 : 1;
    rows_even = (rows % 2) ? 0 : 1;
    hplanes = planes / 2;                    /*  possible bug if planes odd  */
    if (planes % 2) hplanes++;
 
    apodization = (double *)malloc (area * sizeof (double));
    if (edge_outer == 0.0) {
      for (n = 0; n < area; n++) apodization[n] = 1.0;
    } else {
      r0x = hcols + 0.5 * cols_even;
      r0y = hrows + 0.5 * rows_even;
      for (row=0, n=0; row<rows; row++) {
        ry = (row - hrows + 0.5 * rows_even) / r0y;
        for (col=0; col<cols; col++) {
          rx = (col - hcols + 0.5 * cols_even) / r0x;
          r2 = rx * rx + ry * ry;
          r = sqrt (r2);
          if (edge_inner >= edge_outer)
            ra = (r < edge_outer) ? 0.0 : 1.0;
          else
            ra = (r - edge_inner) / (edge_outer - edge_inner);
          ra2 = 1.0 - ra * ra;
          apodization[n++] = (ra >= 1.0) ? 0.0 : (ra <= 0.0) ? 1.0 : ra2 * ra2;
        }
      }
    }
                                          /*  Initialize FFT working space  */
    if (dp_calc) {
      wdata = (double *)malloc (planes * rows * xcols * sizeof (double));
      plan = fftw_plan_dft_r2c_3d (planes, rows, cols, wdata,
          (fftw_complex *)wdata, FFTW_ESTIMATE);
    } else {
      xdata = (float *)malloc (planes * rows * xcols * sizeof (float));
      fplan = fftwf_plan_dft_r2c_3d (planes, rows, cols, xdata,
          (fftwf_complex *)xdata, FFTW_ESTIMATE);
    }
                                     /*  add the argument values to header  */
    if (planes < 2) apode_edge = 2.0;
    bigcube = 0;
    if (cube > 838860800) {
      fprintf (stderr,
          "Warning: call to drms_segment_read()  may fail for %lld values\n",
          cube);
      fprintf (stderr, "         or results may be garbage\n");
      bigcube = 1;
/*
      orig = read_big_cube (iseg, DRMS_TYPE_FLOAT, &status);
*/
      if (dp_calc)
        read_from_big_cube (iseg, dp_calc, wdata, &dataavg, &datavar,
            apode_edge, apodization, xcols);
      else
        read_from_big_cube (iseg, dp_calc, xdata, &dataavg, &datavar,
            apode_edge, apodization, xcols);
    } else {
      orig = drms_segment_read (iseg, DRMS_TYPE_FLOAT, &status);
      if (status) {
        fprintf (stderr, "Error on drms_segment_read in\n");
        fprintf (stderr, "      %s\n", source);
        return cleanup (1, irecs, orec, orig, pspec, dispose);
      }
      l = m = n = 0;
      dataavg = datavar = 0.0;
      data = (float *)orig->data;
      for (plane = 0; plane < planes; plane++) {
        if (apode_edge < 1.0) {
          t = (2.0 * plane) / (planes - 1.0);
          if (t > 1.0) t = 2.0 - t;
          t = 1.0 - t;
          if (t > apode_edge) {
            t = (t - apode_edge) / (1.0 - apode_edge);
            t *= t;
            weight = (1.0 - t);
            weight *= weight;
          } else weight = 1.0;
        } else weight = 1.0;
        if (dp_calc) {
          for (row = 0; row < rows; row++) {
            for (col = 0; col < cols; col++) {
              dval = data[m++];
              if (isnan (dval)) dval = 0.0;
              dval *= weight * apodization[l++];
              dataavg += dval;
              datavar += dval * dval;
              wdata[n++] = dval;
            }
            n += xcols - cols;
          }
        } else {
          for (row = 0; row < rows; row++) {
            for (col = 0; col < cols; col++) {
              fval = data[m++];
              if (isnan (fval)) fval = 0.0;
              fval *= weight * apodization[l++];
              dataavg += fval;
              datavar += fval * fval;
              xdata[n++] = fval;
            }
            n += xcols - cols;
          }
        }
        l = 0;
      }
      drms_free_array (orig);
    }
                                                  /*  necessary for cleanup  */
    orig = NULL;
    dataavg /= m;
    datavar /= m;
    datavar -= dataavg * dataavg;
                                              /*  forward Fourier transform  */
    if (dp_calc) fftw_execute (plan);
    else fftwf_execute (fplan);
                                                 /*  Create output data set  */
    axes[0] = cols;
    axes[1] = rows;
    opln = (fbins > 1) ? hplanes / fbins : hplanes;
    if (fbins > 1 && hplanes % fbins) opln++;
    axes[2] = opln;                          /*  possible bug if planes odd  */
    ntot = cols * rows * opln;
    data = (float *)malloc (ntot * sizeof (float));
                                    /*  copy designated keywords from input  */
    kstat = 0;
    for (key_n = 0; key_n < propct; key_n++) {
      if (strcmp (copykeylist[key_n], "+"))
        kstat += check_and_copy_key (orec, irec, copykeylist[key_n]);
      else kstat += propagate_keys (orec, irec, propagate, keyct);
    }
/*
    kstat = propagate_keys (orec, irec, propagate, keyct);
*/
    dval = drms_getkey_double (irec, "CDELT3", &status);
    if (!status) kstat += check_and_set_key_float (orec, "Cadence", dval);
                                                       /*  and set new ones  */
    crpix = 0.5 * (cols + 1);
    kstat += check_and_set_key_float (orec, "CRPIX1", crpix);
    if (cols != rows) crpix = 0.5 * (rows + 1);
    kstat += check_and_set_key_float (orec, "CRPIX2", crpix);
    crpix = 1.0;
    kstat += check_and_set_key_float (orec, "CRPIX3", crpix);
    kstat += check_and_set_key_float (orec, "CRVAL1", crval);
    kstat += check_and_set_key_float (orec, "CRVAL2", crval);
    kstat += check_and_set_key_float (orec, "CRVAL3", crval);
    kstat += check_and_set_key_str (orec, "CTYPE1", "WAVE-NUM");
    kstat += check_and_set_key_str (orec, "CTYPE2", "WAVE-NUM");
    kstat += check_and_set_key_str (orec, "CTYPE3", "FREQ-ANG");
    kstat += check_and_set_key_float (orec, "Apode_k_min", edge_inner);
    kstat += check_and_set_key_float (orec, "APOD_MIN", edge_inner);
    kstat += check_and_set_key_float (orec, "Apode_k_max", edge_outer);
    kstat += check_and_set_key_float (orec, "APOD_MAX", edge_outer);
    kstat += check_and_set_key_float (orec, "Apode_f", apode_edge);
    if (kstat) {
      fprintf (stderr, "Error writing key value(s) to record in series %s:\n",
          out_series);
      fprintf (stderr, "      series may not have appropriate structure;\n");
      fprintf (stderr, "      record %d skipped\n", rgn);
      continue;
    }

    dval = drms_getkey_double (irec, "MapScale", &status);
    if (!status) {
                                                    /*  assume cols = rows!  */
      dval = 360.0 / dval / cols;                                /*  k (inv radian)  */
      dval /= 696.0;                                         /*  k (inv Mm)  */
      kstat += check_and_set_key_double (orec, "CDELT1", dval);
      if (cols == rows) {
        kstat += check_and_set_key_double (orec, "CDELT2", dval);
        kstat += check_and_set_key_double (orec, "Delta_k", dval);
      } else {
        kstat += check_and_set_key_double (orec, "Delta_kx", dval);
        dval *= (double)cols / (double)rows;
        kstat += check_and_set_key_double (orec, "CDELT2", dval);
        kstat += check_and_set_key_double (orec, "Delta_ky", dval);
      }
      if (kstat) {
        fprintf (stderr, "Error writing key value(s) to record in series %s\n",
          out_series);
        fprintf (stderr, "      series may not have appropriate structure\n");
        fprintf (stderr, "      record %d skipped\n", rgn);
        continue;
      }
    }
    dnu = 1.0e6 / planes / tstep;
    if (fbins > 1) dnu *= fbins;
    kstat += check_and_set_key_double (orec, "Delta_nu", dnu);
    domega = dnu * 2.0e-6 * M_PI;
    kstat += check_and_set_key_double (orec, "D_OMEGA", domega);
    kstat += check_and_set_key_double (orec, "CDELT3", domega);
    if (kstat) {
      fprintf (stderr, "Error writing key value(s) to record in series %s\n",
          out_series);
      fprintf (stderr, "      series may not have appropriate structure\n");
      fprintf (stderr, "      record %d skipped\n", rgn);
      continue;
    }
                /*  convert complex transform to power spectrum and reorder  */
    normal = 1.0 / ntot;
    normal *= normal;
    powrint = 0.0;
    if (fbins > 1) {
      double nfac = normal / fbins / fbins;
      for (plane = 0, n = 0; plane < hplanes; plane += fbins) {
        for (row = 0; row < rows; row++) {
          for (col = 0; col < cols; col++, n++) {
            js = (row == row % hrows) ? row + hrows : row - hrows;
            is = (col == col % hcols) ? col + hcols : col - hcols;
            ftrb = ftib = 0.0;
            for (bin = 0; bin < fbins; bin++) {
                m = (plane + bin) * xcols * rows + js * xcols + 2 * is;
                if (is >= hcols) {
                  m = 2 * (cols - is);
                  if (js != 0) m += (rows - js) * xcols;
                  if ((plane + bin) != 0)
                    m += (planes - plane - bin) * xcols * rows;
                }
                if (dp_calc) {
                  ftrb += wdata[m];
                  ftib += wdata[m+1];
                } else {
                  ftrb += xdata[m];
                  ftib += xdata[m+1];
                }
            }
            data[n] = nfac * (ftrb*ftrb + ftib*ftib);
            if (n == 0) vmax = vmin = vmin2 = data[n];
            if (data[n] > vmax) vmax = data[n];
            if (data[n] < vmin) {
              vmin2 = vmin;
              vmin = data[n];
            }
            powrint += data[n];
          }
        }
      }
      if (n > ntot) {
        fprintf (stderr, "Error: data written beyond output array bounds\n");
        return cleanup (1, irecs, orec, orig, pspec, dispose);
      }
    } else {
      for (plane = 0, n = 0; plane < hplanes; plane++) {
        for (row = 0; row < rows; row++) {
          for (col = 0; col < cols; col++, n++) {
            js = (row == row % hrows) ? row + hrows : row - hrows;
            is = (col == col % hcols) ? col + hcols : col - hcols;
            m = plane * xcols * rows + js * xcols + 2 * is;
            if (is >= hcols) {
                m = 2 * (cols - is);
                if (js != 0) m += (rows - js) * xcols;
                if (plane != 0) m += (planes - plane) * xcols * rows;
            }
            data[n] = (dp_calc) ?
                wdata[m]*wdata[m] + wdata[m+1]*wdata[m+1] :
                xdata[m]*xdata[m] + xdata[m+1]*xdata[m+1];
            data[n] *= normal;
            if (n == 0) vmax = vmin = vmin2 = data[n];
            if (data[n] > vmax) vmax = data[n];
            if (data[n] < vmin) {
              vmin2 = vmin;
              vmin = data[n];
            }
            powrint += data[n];
          }
        }
      }
    }
    if (dp_calc) free (wdata);
    else free (xdata);

    if (verbose) {
      printf ("  P[0] = %11.4e (cf. apodized data mean = %11.4e)\n",
          data[cols * hrows + hcols], dataavg);
      printf ("  S(P) = %11.4e (cf. apodized data var. = %11.4e)\n", powrint, datavar);
    }

    oseg = drms_segment_lookupnum (orec, osegnum);
    if (verbose) printf ("  values range from %.3e to %.3e\n", vmin, vmax);
    if (log_out) {
/*
      double target_min = vmax / exp (scale_range);
printf ("vmin = %.4e <? target_min = %.4e / %.4e = %.4e\n", vmin, vmax,
exp(scale_range), target_min);
*/
      if (vmin < 0.1 * vmin2) {
        fprintf (stderr,
            "** Warning: minimum value = %.2e; value = log (%.2e + data)\n",
            vmin, vmin2);
        for (n = 0; n < ntot; n++) data[n] = log (vmin2 + data[n]);
        vmin = log (vmin2 + vmin);
        vmax = log (vmin2 + vmax);
      } else {
        for (n = 0; n < ntot; n++) data[n] = log (data[n]);
        vmin = log (vmin);
        vmax = log (vmax);
      }
      if (verbose) printf ("  logs scaled between %f and %f\n", vmin, vmax);
      kstat += check_and_set_key_double (orec, "LOG_BASE", exp (1.0));
    }
  
    kstat += check_and_set_key_float (orec, "DataMin", vmin);
    kstat += check_and_set_key_float (orec, "DataMax", vmax);
    kstat += check_and_set_key_str   (orec, "Module", module_ident);
    kstat += check_and_set_key_str   (orec, "BLD_VERS", jsoc_version);
    kstat += check_and_set_key_str   (orec, "Source", source);
    kstat += check_and_set_key_str   (orec, "Input", inser);
    kstat += check_and_set_key_time  (orec, "Created", CURRENT_SYSTEM_TIME);
    if (kstat) {
      fprintf (stderr, "Error writing key value(s) to record in series %s\n",
          out_series);
      fprintf (stderr, "      series may not have appropriate structure\n");
      fprintf (stderr, "      record %d skipped\n", rgn);
      continue;
    }
    bscale = (vmax - vmin) / scale_range;
    bzero = 0.5 * (vmax + vmin);
    pspec = drms_array_create (DRMS_TYPE_FLOAT, 3, axes, data, &status);
    if (status) {
      fprintf (stderr, "Error: couldn't create output array\n");
      return cleanup (1, irecs, orec, orig, pspec, dispose);
    }
    pspec->bscale = bscale;
    pspec->bzero = bzero;
    if (bigcube) {
      fprintf (stderr,
          "Warning: call to drms_segment_write()  may fail for %lld values\n",
          ntot);
      fprintf (stderr, "         or results may be garbage\n");
    }
    status = drms_segment_write (oseg, pspec, 0);
    if (status) {
      drms_sprint_rec_query (recid, orec);
      fprintf (stderr, "Error writing data to record %s\n", recid);
      fprintf (stderr, "      series may not have appropriate structure\n");
      return cleanup (1, irecs, orec, orig, pspec, dispose);
    }
    drms_free_array (pspec);
    pspec = NULL;

    if (verbose) {
      drms_sprint_rec_query (recid, orec);
      if (no_save)
        printf ("output data record would be %s\n", recid);
      else {
        drms_segment_filename (oseg, pathname);
        printf ("output data record is %s\n", recid);
        printf ("output data cube at %s\n", pathname);
      }
    }
    drms_close_record (orec, dispose);
  }
                                                      /* Clean up FFTW stuff */
  fftwf_destroy_plan (fplan);
/*  this is kind of silly, why not just remove error argument from cleanup?  */
  return cleanup (0, irecs, NULL, orig, pspec, dispose);
}

/*
 *  Revision History (all mods by Rick Bogart unless otherwise noted)
 *
 *  08.03.24    created this file, based on SOI powrspec3_v20
 *  v 0.6 frozen 2009.07.14
 *  v 0.7:      fix to run with NetDRMS 2.0 (09.07.14);
 *              Modify and add to list of keywords copied from input
 *  v 0.7 frozen 2009.09.08.05
 *  v 0.8:      general code cleanup; added copying of all standard (and
 *      additional "standard") keys;
 *              removed defaults for I/O;
 *              removed CGI output and special processing for exports;
 *              fixed bug in setting of CDELT3 and associated keywords;
 *              modified key copying to use list and to be more careful;
 *              added setting of Cadence keywords from input CDELT3;
 *              put in (approximately) correct normalization;
 *              generalized scaling of logs (and direct values) for output
 *      precision
 *              added looping over multiple input and output records
 *              added no_save option for diagnostics
 *  v 0.8 frozen 2010.01.18
 *  v 0.9:      fixed use of params_get_str
 *              added setting of several keys: bld_vers, created
 *              added setting of keys APOD_MIN, APOD_MAX
 *              added Width & Height to and removed Cadence from default
 *      propagation list
 *              added optional segment argument and removed restriction on
 *      its name
 *              added warning about large input cubes
 *  v 0.9 frozen 2010.04.23
 *  v 1.0       Added code for reading and apodizing large cubes by slices;
 *      however, there is still a problem with the memory required for the
 *      output array at the same time the complex array is in place
 *              Removed an extraneous keyword propagation
 *              Added Ident to list of default propagated keywords; added option
 *      for providing alternate or additional list of keywords for propagation
 *  v 1.0 frozen 2010.08.19
 *  v 1.1       Removed some old commented out code
 *      Fixed bug in test for minimum values out of scaling range (zero)
 *  v 1.1 frozen 2012.12.24
 *
 */
Revision:	1.4
Committed:	Mon Apr 8 22:27:55 2013 UTC (10 years, 5 months ago) by rick
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Ver_9-1, Ver_9-3, Ver_9-2, Ver_8-8, Ver_8-2, Ver_8-3, Ver_8-0, Ver_8-1, Ver_8-6, Ver_8-7, Ver_8-4, Ver_8-5, Ver_8-10, Ver_8-11, Ver_8-12, Ver_9-0
Changes since 1.3:	+25 -40 lines
Log Message:	modified for JSOC release 8.0