JSOC_Documentation: proj/rings/apps/rdfitc.c Source File

00001 /*
00002  *  rdfitc.c                        ~rick/src/ringfit/rdfitc
00003  *
00004  *  Responsible:  Rick Bogart                       RBogart@spd.aas.org
00005  *     
00006  *  "Ring" fitting based on the code developed and used by Sarbani Basu
00007  *    and H. M. Antia for fitting 3-d solar acoustic power spectra; based
00008  *    on SOI module ringfit_v0.86
00009  *
00010  *  Main program begins around line 800
00011  *
00012  *  Bugs:
00013  *    Timing is not re-zeroed for each I/O record pair processed
00014  *    default values of lmin and lmax are appropriate for HMI resolution
00015  *  and areas of size 15 deg
00016  *    non-default value of nmax seems to force premature termination; in
00017  *  particular, nmax <= 3 causes crash due to going beyond array bounds
00018  *  of freq, wid
00019  *    It is implicitly assumed (not checked) that the spectrum is centred,
00020  *  with CRPIX1 = cols / 2 + 0.5, CRPIX2 = rows / 2 + 0.5, CRPIX3 = 1,
00021  *  CRVALi = 0 (or the equivlent), and that the units are in Mm-1 and
00022  *  sec-1
00023  *
00024  *  Planned development:
00025  *    parametrize test criteria for bad fits
00026  *    Use drms_keyword_lookup()
00027  *    Remove unnecessary defined constants, predefined size declarations
00028  *    Add parameters for fitting coefficients
00029  *    Remove unused variables, provide more descriptive names
00030  *    Generalize bfgs() to use function call in parameter
00031  *
00032  *  Revision history is at end of file.
00033  */
00034 
00035 #include <jsoc_main.h>
00036                               /*  module identifier  */
00037 char *module_name = "ringfit_bba";
00038 char *module_desc = "ring fitting using method of Basu and Antia";
00039 char *version_id = "1.2";
00040 
00041 ModuleArgs_t module_args[] = {
00042   {ARG_DATASET, "in", "", "input dataset"}, 
00043   {ARG_STRING, "guessfile", "",
00044       "name of file containing initial frequency guesses"},
00045   {ARG_STRING, "out", "fit.out",
00046       "name of output data series (or file) containing fit coefficients"},
00047   {ARG_INT, "nmin", "0", "lowest radial order to fit"},
00048   {ARG_INT, "nmax", "10", "highest radial order to fit"},
00049   {ARG_INT, "lmin", "80", "lowest degree to fit"},
00050   {ARG_INT, "lmax", "3000", "highest degree to fit"},
00051   {ARG_INT, "fmin", "900", "lowest frequency to fit [uHz]"},
00052   {ARG_INT, "fmax", "5500", "highest frequency index to fit"},
00053   {ARG_INT, "bfgsct", "125", "number of iterations for bfgs"},
00054   {ARG_INT, "linminct", "15", "number of iterations for linmin"},
00055   {ARG_FLOAT, "ux", "0.0", "initial guess for zonal flow speed [m/s]"},
00056   {ARG_FLOAT, "uy", "0.0", "initial guess for merridional flow speed [m/s]"},
00057   {ARG_FLOAT, "A1", "fit", "fixed value for amplitude term A1 in place of fit"},
00058   {ARG_FLOAT, "S", "fit", "fixed value for asymmetry term in place of fit"},
00059   {ARG_FLOAT, "A1_guess", "0.0", "initial guess for amplitude term A1"},
00060   {ARG_FLOAT, "S_guess", "-200.0", "initial guess for asymmetry term S"},
00061   {ARG_FLOAT, "kxmin", "unspecified", "minimum cutoff for kx in spectrum"},
00062   {ARG_FLOAT, "kxmax", "unspecified", "maximum cutoff for kx in spectrum"},
00063   {ARG_FLOAT, "kymin", "unspecified", "minimum cutoff for ky in spectrum"},
00064   {ARG_FLOAT, "kymax", "unspecified", "maximum cutoff for ky in spectrum"},
00065   {ARG_STRING, "copy",  "+",
00066       "comma separated list of keys to propagate forward"},
00067   {ARG_FLAG, "n", "0", "no fitting (diagnostics only)"},      
00068   {ARG_FLAG, "v", "0", "verbose output"},      
00069   {ARG_FLAG, "x", "0", "extra reporting of number of iterations"},      
00070   {ARG_FLAG, "2", "0", "calculate second derivatives explicitly"},      
00071   {}
00072 };
00073        /*  list of keywords to propagate (if possible) from input to output  */
00074 char *propagate[] = {"CarrTime", "CarrRot", "CMLon", "LonHG", "LatHG", "LonCM",
00075     "MidTime", "Duration", "Cadence", "LonSpan",
00076     "T_START", "T_STOP", "Coverage", "Size", "Width", "Height",
00077     "ZonalTrk", "ZonalVel", "MeridTrk", "MeridVel",
00078     "MapScale", "MapProj", "Map_PA", "RSunRef", "PosAng", "MAI", "Ident",
00079     "Apode_f", "Apode_k_min", "Apode_k_max", "APODIZNG", "APOD_MIN", "APOD_MAX"};
00080 
00081 #include <unistd.h>
00082 #include <math.h>
00083 #include <time.h>
00084 #include "keystuff.c"
00085 #include "rdutil.c"
00086 
00087 #define VSCALE  (-0.4372)
00088 #define RSUN_MM (696.0)
00089 
00090 static char *time_elapsed (int cct, int *sum) {
00091 /*
00092  *  print a string containing utime reported to millisec plus wall-clock
00093  *    time since last zeroing call
00094  *  allows for one wrap-around in clock, which occurs every 36 minutes,
00095  *    but results will still be garbage if consecutive calls are spaced
00096  *    more than 72 minutes apart.
00097  */
00098   static clock_t t_last, t_curr;
00099   static time_t wt_last, wt_curr;
00100   static double cloadavg[3], loadavg[3];
00101   static double t_count, scale = 1.0 / CLOCKS_PER_SEC;
00102   static int loadct, first = 1;
00103   static char str[64];
00104   int n;
00105 
00106   if (first) {
00107     t_last = clock ();
00108     wt_last = time (NULL);
00109     first = 0;
00110     *sum = 0;
00111     for (n = 0; n < 3; n++) loadavg[n] = 0.0;
00112     loadct = 0;
00113   }
00114   t_curr = clock ();
00115   t_count += scale * (t_curr - t_last);
00116 /*
00117   if (t_curr < t_last) t_count += 4294.9673;
00118 */
00119   wt_curr = time (NULL);
00120   getloadavg (cloadavg, 3);
00121   for (n = 0; n < 3; n++) loadavg[n] += cloadavg[n];
00122   loadct++;
00123   
00124   sprintf (str, "%.3f (%d) sec (avg load: %.2f %.2f %.2f)", t_count,
00125       (int)(wt_curr - wt_last), loadavg[0]/loadct, loadavg[1]/loadct,
00126       loadavg[2]/loadct);
00127   t_last = t_curr;
00128   if (!cct) {
00129     *sum += wt_curr - wt_last;
00130     t_last = clock ();
00131     wt_last = time (NULL);
00132     t_count = 0.0;
00133     for (n = 0; n < 3; n++) loadavg[n] = 0.0;
00134     loadct = 0;
00135   }
00136   return str;
00137 }
00138 
00139 static int gauelm (int n, int num, double *a, double *x, double *det,
00140     int *indx, int lj, int flag) {
00141 
00142   double r1, t1;
00143   int i, j, k, km, l;
00144 
00145   if (n <= 0 || n > lj) return 101;
00146 
00147   if (flag < 2) {
00148                             /*  perform elimination  */
00149     *det = 1.0;
00150     for (k = 0; k < n-1; k++) {
00151                        /*  find maximum element in column k  */
00152       r1 = 0.0;
00153       km = k;
00154       for (l = k; l < n; l++)
00155     if (fabs (a[l*lj + k]) > r1) {
00156       r1 = fabs (a[l*lj + k]);
00157       km = l;
00158     }
00159       indx[k] = km;
00160       if (km != k) {
00161                       /*  interchange the rows if necessary  */
00162     for (l = k; l < n; l++) {
00163       t1 = a[k*lj + l];
00164       a[k*lj + l] = a[km*lj + l];
00165       a[km*lj + l] = t1;
00166     }
00167     *det *= -1;
00168       }
00169       *det *= a[k*lj + k];
00170       if (a[k*lj + k] == 0.0) return 121;
00171       for (l = k+1; l < n; l++) {
00172     a[l*lj + k] /= a[k*lj + k];
00173     for (i = k+1; i < n; i++)
00174       a[l*lj + i] -= a[l*lj + k] * a[k*lj + i];
00175       }
00176     }
00177     *det *= a[(n-1)*lj + n-1];
00178     indx[n-1] = n - 1;
00179     if (a[(n-1)*lj + n-1] == 0) return 121;
00180     if (flag == 1) return 0;
00181   }
00182             /*  solution for the num different right-hand sides  */
00183   for (j = 0; j < num; j++) {
00184                            /*  forward-substitution  */
00185     int ljj = lj * j;
00186     for (k = 0; k < n-1; k++) {
00187       if (k != indx[k]) {
00188     t1 = x[ljj + k];
00189     x[ljj + k] = x[ljj + indx[k]];
00190     x[ljj + indx[k]] = t1;
00191       }
00192       for (l = k+1; l < n; l++)
00193     x[ljj + l] -= a[l*lj + k] * x[ljj + k];
00194     }
00195                               /*  back-substitution  */
00196     x[ljj + n-1] /= a[(n-1)*lj + n-1];
00197     for (k = n-2; k >= 0; k--) {
00198       for (l = n-1; l >= k+1; l--)
00199     x[ljj + k] -= x[ljj + l] * a[k*lj + l];
00200       x[ljj + k] /= a[k*lj + k];
00201     }
00202   }
00203   return 0;
00204 }
00205 
00206 static void fcn1 (int n, double *x, double *f, double *g,
00207     double *kx, double *ky, double *k2, double *omeg, double *ppow,
00208     double pk, double *fm, double ki, int npix) {
00209 /*
00210  *  Given values for the set of fitting parameters x[13], form the sums:
00211  *    log (likelihood function) = f = sum of log M_i + O_i / M_i
00212  *  (M_i are fits for the n data values kx, ky, omeg;
00213  *   O_i are the n values of power, ppow[i]);
00214  *    merit function = fm = sum of (O_i / M_i - 1)^2;
00215  *    g[13] = sum of ? (partial derivatives of fm w.r.t parameters?)
00216  */
00217   static double vf = 1.0e-4;
00218   double pk1, k, kx2, ky2, bk, dk, om1, wd, o1, d1, p1;
00219   double logk, kpk1, bko1p1, bktrm;
00220   double p1x0, p1x1, p1x2, p1x3, p1x4, p1x5, p1x8, p1x9,
00221       p1x10, p1x11, p1x12, p1x13, pix11, pix12, pi, dpix;
00222   int i;
00223                             /*  initialize sums  */
00224   *f = 0.0;
00225   *fm = 0.0;
00226   for (i = 0; i < n; i++) g[i] = 0.0;
00227                      /*  limit values of fitting parameters  */
00228 /*
00229   if (x[0] > 100.0) x[0] = 100.0;
00230   if (x[0] < 1.0) x[0] = 1.0;
00231   if (x[5] > 100.0) x[5] = 100.0;
00232   if (x[5] < 1.0) x[5] = 1.0;
00233   if (x[6] > 100.0) x[6] = 100.0;
00234   if (x[6] < 1.0) x[6] = 1.0;
00235   if (x[6] < 20.0 && x[5] < 20.0) x[6] = 20.0;
00236 */
00237   if (x[7] > 1000.0) x[7] = 1000.0;
00238   if (x[7] < -1000.0) x[7] = -1000.0;
00239   if (x[8] > 8.0e4) x[8] = 8.0e4;
00240   if (x[8] < -8.0e4) x[8] = -8.0e4;
00241   if (x[9] > 5.0e4) x[9] = 5.0e4;
00242   if (x[9] < -5.0e4) x[9] = -5.0e4;
00243   if (x[10] > 4000.0) x[10] = 4000.0;
00244   if (x[10] < -2000.0) x[10] = -2000.0;
00245   if (x[11] > 2.0e4) x[11] = 2.0e4;
00246   if (x[11] < -2.0e4) x[11] = -2.0e4;
00247   if (x[12] > 5000.0) x[12] = 5000.0;
00248   if (x[12] < -5000.0) x[12] = -5000.0;
00249 
00250   pk1 = pk + x[10] * vf;
00251   bk = x[12] * vf;
00252 
00253   for (i = 0; i < npix; i++) {
00254     kx2 = kx[i] * kx[i];
00255 /*
00256     k = hypot (kx[i], ky[i]);
00257     k2 = k * k;
00258     k2 = kx2 + ky[i] * ky[i];
00259 */
00260     k = sqrt (k2[i]);
00261     kx2 /= k2[i];
00262     ky2 = kx[i] * ky[i] / k2[i];
00263     dk = k - ki;
00264 /*
00265     kpk1 = pow (k, pk1);
00266 */
00267     logk = log (k);
00268     kpk1 = exp (pk1 * logk);
00269 
00270     om1 = x[1] * kpk1 * logk * vf;
00271     wd = x[4] + x[11] * dk * vf;
00272     o1 = (omeg[i] - x[1] * kpk1 - x[2] * vf * kx[i]
00273         - x[3] * vf * ky[i]) / wd;
00274     bko1p1 = 1.0 + bk * o1;
00275     bktrm = bk * bk + bko1p1 * bko1p1;
00276     d1 = o1 * o1 + 1.0;
00277     p1 = exp (x[0] + vf * (x[7] * dk + x[8] * kx2
00278         + x[9] * ky2)) * bktrm / d1;
00279     p1x1 = p1;
00280     p1x0 = p1 * 2 * bk * bko1p1 / bktrm - 2 * p1 * o1 / d1;
00281     p1x0 /= wd;
00282     p1x2 = -p1x0 * kpk1;
00283     p1x3 = -p1x0 * kx[i] * vf;
00284     p1x4 = -p1x0 * ky[i] * vf;
00285     p1x5 = -p1x0 * o1;
00286     p1x8 = p1 * dk * vf;
00287     p1x9 = p1 * kx2 * vf;
00288     p1x10 = p1 * ky2 * vf;
00289     p1x11 = -p1x0 * om1;
00290     p1x12 = -p1x0 * o1 * dk * vf;
00291     p1x13 = p1 * 2 * (bk + bko1p1 * o1) * vf / bktrm;
00292     pix11 = exp (x[5]) / (k2[i] * k);
00293     pix12 = exp (x[6]) / (k2[i] * k2[i]);
00294 
00295     pi = p1 + pix11 + pix12;
00296     *f += log (pi) + ppow[i] / pi;
00297     dpix = 1.0 - ppow[i] / pi;
00298     *fm += dpix * dpix;
00299     dpix /= pi;
00300     g[0] += p1x1 * dpix;
00301     g[1] += p1x2 * dpix;
00302     g[2] += p1x3 * dpix;
00303     g[3] += p1x4 * dpix;
00304     g[5] += pix11 * dpix;      
00305     g[6] += pix12 * dpix;      
00306     g[4] += p1x5 * dpix;
00307     g[7] += p1x8 * dpix;
00308     g[8] += p1x9 * dpix;
00309     g[9] += p1x10 * dpix;
00310     g[10] += p1x11 * dpix;
00311     g[11] += p1x12 * dpix;
00312     g[12] += p1x13 * dpix;
00313   }
00314 }
00315 
00316 static void d2fcn (int n, double *x, double *f, double *g, double *dd,
00317     double *k,  double *kx, double *ky, double *k2, double *omeg,
00318     double *ppow, double pk, double ki, int npix) {
00319   static double vf = 1.0e-4;
00320   double dp1[20*20];
00321   double bk, bkf, d1, om1, o1, pi, pix11, pix12, pk1, wd;
00322   double p1, p1x0, p1x1, p1x2, p1x3, p1x4, p1x5, p1x8, p1x9,
00323       p1x10, p1x11, p1x12, p1x13;
00324   double pi2, spow, pow2m1, scale;
00325   double dk, dp1xx, dp1x6;
00326   double bko1p1, dptrm;
00327   double k3, k4, kx2, ky2;
00328   double kxk, kpk1, logk;
00329   int i, j, p, col, row;
00330 
00331   *f = 0.0;
00332   for (j = 0; j < n; j++) {
00333     row = n * j;
00334     g[j] = 0.0;
00335     for (i = 0; i < n; i++)
00336       dd[i + row] = 0.0;
00337   }
00338   pk1 = x[10] * vf + pk;
00339 
00340   for (i = 0; i < npix; i++) {
00341 /*
00342     k2 = kx[i] * kx[i] + ky[i] * ky[i];
00343     k = sqrt (k2[i]);
00344 */
00345     k3 = k2[i] * k[i];
00346     k4 = k2[i] * k2[i];
00347     kxk = kx[i] / k[i];
00348     kx2 = vf * kxk * kxk;
00349     ky2 = vf * kxk * ky[i] / k[i];
00350     dk = k[i] - ki;
00351     logk = log (k[i]);
00352     kpk1 = exp (pk1 * logk);
00353 
00354     bk = x[12] * vf;
00355     om1 = x[1] * kpk1 * logk * vf;
00356     wd = x[4] + x[11] * dk * vf;
00357     o1 = (omeg[i] - x[1] * kpk1 - x[2] * vf * kx[i] -
00358         x[3] * vf * ky[i]) / wd;
00359     d1 = o1 * o1 + 1.0;
00360     bko1p1 = 1.0 + bk * o1;
00361     bkf = bk * bk + bko1p1 * bko1p1;
00362     p1 = exp (x[0] + vf * x[7] * dk + x[8] * kx2 +
00363         x[9] * ky2) * bkf / d1;
00364     p1x0 = 2 * ((p1 * bk * bko1p1 / bkf) - p1 * o1 / d1);
00365     p1x0 /= wd;
00366     p1x1 = p1;
00367     p1x2 = -p1x0 * kpk1;
00368     p1x3 = -p1x0 * (kx[i] * vf);
00369     p1x4 = -p1x0 * (ky[i] * vf);
00370     p1x5 = -p1x0 * o1;
00371     p1x8 = p1 * dk * vf;
00372     p1x9 = p1 * kx2;
00373     p1x10 = p1 * ky2;
00374     p1x11 = -p1x0 * om1;
00375     p1x12 = -p1x0 * o1 * dk * vf;
00376     p1x13 = p1 * 2 * (bk + bko1p1 * o1) * vf / bkf;
00377     dp1xx = (2 * p1) * (bk * bk / bkf -
00378         4 * o1 * bk * bko1p1 / (d1 * bkf) +
00379     4 * o1 * o1 / d1 / d1 - 1 / d1);
00380     dp1x6 = (4*o1*bk + 2 - 2*o1*(2*bk + 2*o1*bko1p1) / d1) * vf * p1 / bkf;
00381     dp1[0 + 20*0] = p1x1;
00382     dp1[0 + 20*1] = p1x2;
00383     dp1[0 + 20*2] = p1x3;
00384     dp1[0 + 20*3] = p1x4;
00385     dp1[0 + 20*4] = p1x5;
00386     dp1[0 + 20*7] = p1x8;
00387     dp1[0 + 20*8] = p1x9;
00388     dp1[0 + 20*9] = p1x10;
00389     dp1[0 + 20*10] = p1x11;
00390     dp1[0 + 20*11] = p1x12;
00391     dp1[0 + 20*12] = p1x13;
00392     dp1[1 + 20*(0)] = p1x2;
00393     dp1[2 + 20*(0)] = p1x3;
00394     dp1[3 + 20*(0)] = p1x4;
00395     dp1[4 + 20*(0)] = p1x5;
00396     dp1[7 + 20*(0)] = p1x8;
00397     dp1[8 + 20*(0)] = p1x9;
00398     dp1[9 + 20*(0)] = p1x10;
00399     dp1[10 + 20*(0)] = p1x11;
00400     dp1[11 + 20*(0)] = p1x12;
00401     dp1[12 + 20*(0)] = p1x13;
00402     dptrm = dp1xx * o1 / wd + p1x0;
00403     dp1[1 + 20*(1)] = dp1xx * (kpk1 / wd) * (kpk1 / wd);
00404     dp1[1 + 20*(2)] = dp1xx * kpk1 * kx[i] * vf / wd / wd;
00405     dp1[1 + 20*(3)] = dp1xx * kpk1 * ky[i] * vf / wd / wd;
00406     dp1[1 + 20*(4)] = dptrm * kpk1 / wd;
00407     dp1[1 + 20*(7)] = p1x2 * dk * vf;
00408     dp1[1 + 20*(8)] = p1x2 * kx2;
00409     dp1[1 + 20*(9)] = p1x2 * ky2;
00410     dp1[1 + 20*(10)] = dp1xx * om1 * kpk1 / wd / wd +
00411         p1x0 * om1 / x[1];
00412     dp1[1 + 20*(11)] = dptrm * dk * vf * kpk1 / wd;
00413     dp1[1 + 20*(12)] = -dp1x6 * kpk1 / wd;
00414     dp1[2 + 20*(1)] = dp1[1 + 20*(2)];
00415     dp1[3 + 20*(1)] = dp1[1 + 20*(3)];
00416     dp1[4 + 20*(1)] = dp1[1 + 20*(4)];
00417     dp1[7 + 20*(1)] = dp1[1 + 20*(7)];
00418     dp1[8 + 20*(1)] = dp1[1 + 20*(8)];
00419     dp1[9 + 20*(1)] = dp1[1 + 20*(9)];
00420     dp1[10 + 20*(1)] = dp1[1 + 20*(10)];
00421     dp1[11 + 20*(1)] = dp1[1 + 20*(11)];
00422     dp1[12 + 20*(1)] = dp1[1 + 20*(12)];
00423     dp1[2 + 20*(2)] = dp1xx * (kx[i] * vf / wd) * (kx[i] * vf / wd);
00424     dp1[2 + 20*(3)] = dp1xx * kx[i] * ky[i] * (vf / wd) * (vf / wd);
00425                /*  changed from presumed erroneous code in original  */
00426     dp1[2 + 20*(4)] = dptrm * kx[i] * vf / wd;
00427     dp1[2 + 20*(7)] = p1x3 * dk * vf;
00428     dp1[2 + 20*(8)] = p1x3 * kx2;
00429     dp1[2 + 20*(9)] = p1x3 * ky2;
00430     dp1[2 + 20*(10)] = dp1xx * kx[i] * vf * om1 / wd / wd;
00431     dp1[2 + 20*(11)] = dptrm * dk * kx[i] * vf *vf / wd;
00432     dp1[2 + 20*(12)] = -dp1x6 * kx[i] * vf / wd;
00433     dp1[3 + 20*(2)] = dp1[2 + 20*(3)];
00434     dp1[4 + 20*(2)] = dp1[2 + 20*(4)];
00435     dp1[7 + 20*(2)] = dp1[2 + 20*(7)];
00436     dp1[8 + 20*(2)] = dp1[2 + 20*(8)];
00437     dp1[9 + 20*(2)] = dp1[2 + 20*(9)];
00438     dp1[10 + 20*(2)] = dp1[2 + 20*(10)];
00439     dp1[11 + 20*(2)] = dp1[2 + 20*(11)];
00440     dp1[12 + 20*(2)] = dp1[2 + 20*(12)];
00441     dp1[3 + 20*(3)] = dp1xx * (ky[i] * vf / wd) * (ky[i] * vf / wd);
00442     dp1[3 + 20*(4)] = dptrm * ky[i] * vf / wd;
00443     dp1[3 + 20*(7)] = p1x4 * dk * vf;
00444     dp1[3 + 20*(8)] = p1x4 * kx2;
00445     dp1[3 + 20*(9)] = p1x4 * ky2;
00446     dp1[3 + 20*(10)] = dp1xx * ky[i] * vf * om1 / wd / wd;
00447     dp1[3 + 20*(11)] = dptrm * dk * ky[i] * vf * vf / wd;
00448     dp1[3 + 20*(12)] = -dp1x6 * ky[i] * vf / wd;
00449     dp1[4 + 20*(3)] = dp1[3 + 20*(4)];
00450     dp1[7 + 20*(3)] = dp1[3 + 20*(7)];
00451     dp1[8 + 20*(3)] = dp1[3 + 20*(8)];
00452     dp1[9 + 20*(3)] = dp1[3 + 20*(9)];
00453     dp1[10 + 20*(3)] = dp1[3 + 20*(10)];
00454     dp1[11 + 20*(3)] = dp1[3 + 20*(11)];
00455     dp1[12 + 20*(3)] = dp1[3 + 20*(12)];
00456     dp1[4 + 20*(4)] = (dptrm + p1x0) * o1 / wd;
00457     dp1[4 + 20*(7)] = p1x5 * dk * vf;
00458     dp1[4 + 20*(8)] = p1x5 * kx2;
00459     dp1[4 + 20*(9)] = p1x5 * ky2;
00460     dp1[4 + 20*(10)] = dptrm * om1 / wd;
00461     dp1[4 + 20*(11)] = (dptrm + p1x0) * dk * o1 * vf / wd;
00462     dp1[4 + 20*(12)] = -dp1x6 * o1 / wd;
00463     dp1[7 + 20*(4)] = dp1[4 + 20*(7)];
00464     dp1[8 + 20*(4)] = dp1[4 + 20*(8)];
00465     dp1[9 + 20*(4)] = dp1[4 + 20*(9)];
00466     dp1[10 + 20*(4)] = dp1[4 + 20*(10)];
00467     dp1[11 + 20*(4)] = dp1[4 + 20*(11)];
00468     dp1[12 + 20*(4)] = dp1[4 + 20*(12)];
00469     dp1[7 + 20*(7)] = p1 * dk * dk * vf * vf;
00470     dp1[7 + 20*(8)] = p1 * dk * kx2;
00471     dp1[7 + 20*(9)] = p1 * dk * ky2;
00472     dp1[7 + 20*(10)] = p1x11 * dk * vf;
00473     dp1[7 + 20*(11)] = p1x12 * dk * vf;
00474     dp1[7 + 20*(12)] = p1x13 * dk * vf;
00475     dp1[8 + 20*(7)] = dp1[7 + 20*(8)];
00476     dp1[9 + 20*(7)] = dp1[7 + 20*(9)];
00477     dp1[10 + 20*(7)] = dp1[7 + 20*(10)];
00478     dp1[11 + 20*(7)] = dp1[7 + 20*(11)];
00479     dp1[12 + 20*(7)] = dp1[7 + 20*(12)];
00480     dp1[8 + 20*(8)] = p1 * kx2 * kx2;
00481     dp1[8 + 20*(9)] = p1 * kx2 * ky2;
00482     dp1[8 + 20*(10)] = p1x11 * kx2;
00483     dp1[8 + 20*(11)] = p1x12 * kx2;
00484     dp1[8 + 20*(12)] = p1x13 * kx2;
00485     dp1[9 + 20*(8)] = dp1[8 + 20*(9)];
00486     dp1[10 + 20*(8)] = dp1[8 + 20*(10)];
00487     dp1[11 + 20*(8)] = dp1[8 + 20*(11)];
00488     dp1[12 + 20*(8)] = dp1[8 + 20*(12)];
00489     dp1[9 + 20*(9)] = p1 * ky2 * ky2;
00490     dp1[9 + 20*(10)] = p1x11 * ky2;
00491     dp1[9 + 20*(11)] = p1x12 * ky2;
00492     dp1[9 + 20*(12)] = p1x13 * ky2;
00493     dp1[10 + 20*(9)] = dp1[9 + 20*(10)];
00494     dp1[11 + 20*(9)] = dp1[9 + 20*(11)];
00495     dp1[12 + 20*(9)] = dp1[9 + 20*(12)];
00496     dp1[10 + 20*(10)] = dp1xx * (om1 / wd) * (om1 / wd) -
00497         p1x0 * om1 * logk * vf;
00498     dp1[10 + 20*(11)] = dptrm * om1 * dk * vf / wd;
00499     dp1[10 + 20*(12)] = -dp1x6 * om1 / wd;
00500     dp1[11 + 20*(10)] = dp1[10 + 20*(11)];
00501     dp1[12 + 20*(10)] = dp1[10 + 20*(12)];
00502     dp1[11 + 20*(11)] = (dptrm + p1x0) * o1 *
00503         (dk * vf / wd) * dk * vf;
00504     dp1[11 + 20*(12)] = -dp1x6 * o1 * dk * vf / wd;
00505     dp1[12 + 20*(11)] = dp1[11 + 20*(12)];
00506     dp1[12 + 20*(12)] = 2 * p1 * vf * vf * d1 / bkf;
00507 
00508     pix11 = exp(x[5]) / k3;
00509     pix12 = exp(x[6]) / k4;
00510     pi = 1.0 / (p1 + pix11 + pix12);
00511     pi2 = pi * pi;
00512     spow = pi * ppow[i];
00513     *f += -log (pi) + spow;
00514 
00515     p1x1 *= pi;
00516     p1x2 *= pi;
00517     p1x3 *= pi;
00518     p1x4 *= pi;
00519     p1x5 *= pi;
00520     p1x8 *= pi;
00521     p1x9 *= pi;
00522     p1x10 *= pi;
00523     p1x11 *= pi;
00524     p1x12 *= pi;
00525     pix11 *= pi;
00526     pix12 *= pi;
00527     pow2m1 = spow + spow - 1.0;
00528     g[0] += p1x1 - spow * p1x1;
00529     g[1] += p1x2 - spow * p1x2;
00530     g[2] += p1x3 - spow * p1x3;
00531     g[3] += p1x4 - spow * p1x4;
00532     g[4] += p1x5 - spow * p1x5;
00533     g[5] += pix11 - spow * pix11;
00534     g[6] += pix12 - spow * pix12;
00535     g[7] += p1x8 - spow * p1x8;
00536     g[8] += p1x9 - spow * p1x9;
00537     g[9] += p1x10 - spow * p1x10;
00538     g[10] += p1x11 - spow * p1x11;
00539     g[11] += p1x12 - spow * p1x12;
00540     g[12] += p1x13 - spow * p1x13;
00541 
00542     for (j = 0; j < n; j++) {
00543       int rows = 20 * j;
00544       for (col = 0; col < n; col++) {
00545         int cols = 20 * col;
00546     if ((j != 6) && (j != 5) && (col != 6) && (col != 5))
00547       dd[j + n*col] += dp1[j + cols] * pi - dp1[rows] * dp1[cols] * pi2 -
00548           spow * dp1[j + cols] * pi +
00549           2 * spow * dp1[rows] * dp1[cols] * pi2;
00550       }
00551     }
00552 
00553     dd[5 + n*5] += pix11 * (1 + spow) + pix11 * pix11 * (2 * spow - 1);
00554     dd[6 + n*(6)] += pix12 * (1 + spow) + pix12 * pix12 * (2 * spow - 1);
00555     scale = pix11 * pow2m1;
00556 /*
00557     for (p = 0; p < 13; p++) {
00558       if (p == 5) continue;
00559       dd[5 + n*p] += scale * dpdx[p];
00560     }
00561 */
00562     dd[5 + n*(0)] += scale * p1x1;
00563     dd[5 + n*(1)] += scale * p1x2;
00564     dd[5 + n*(2)] += scale * p1x3;
00565     dd[5 + n*(3)] += scale * p1x4;
00566     dd[5 + n*(4)] += scale * p1x5;
00567     dd[5 + n*(6)] += scale * pix12;
00568     dd[5 + n*(7)] += scale * p1x8;
00569     dd[5 + n*(8)] += scale * p1x9;
00570     dd[5 + n*(9)] += scale * p1x10;
00571     dd[5 + n*(10)] += scale * p1x11;
00572     dd[5 + n*(11)] += scale * p1x12;
00573     dd[5 + n*(12)] += scale * p1x13;
00574     row = 5 * n;
00575     for (p = 0; p < 13; p++) {
00576       if (p == 5) continue;
00577       dd[p + row] = dd[5 + p*n];
00578     }
00579     scale = pix12 * pow2m1;
00580 /*
00581     for (p = 0; p < 13; p++) {
00582       if (p == 5) continue;
00583       dd[6 + n*p] += scale * dpdx[p];
00584     }
00585 */
00586     dd[6 + n*(0)] += scale * p1x1;
00587     dd[6 + n*(1)] += scale * p1x2;
00588     dd[6 + n*(2)] += scale * p1x3;
00589     dd[6 + n*(3)] += scale * p1x4;
00590     dd[6 + n*(4)] += scale * p1x5;
00591     dd[6 + n*(7)] += scale * p1x8;
00592     dd[6 + n*(8)] += scale * p1x9;
00593     dd[6 + n*(9)] += scale * p1x10;
00594     dd[6 + n*(10)] += scale * p1x11;
00595     dd[6 + n*(11)] += scale * p1x12;
00596     dd[6 + n*(12)] += scale * p1x13;
00597     row = 6 * n;
00598     for (p = 0; p < 13; p++) {
00599       if (p == 5) continue;
00600       dd[p + row] = dd[6 +p*n];
00601     }
00602   }
00603 }
00604 
00605 /*
00606 static double flnm (void (*fcn)(int, double *, double *, double *), 
00607     double x, double *df, double *v, double *x0, int n) {
00608 */
00609 static double flnm (double x, double *df, double *v, double *x0, int n,
00610     double *rkx, double *rky, double *k2, double *omeg, double *ppow,
00611     double pk, double *fm, double rki, int npix, int *flnm_ct) {
00612 
00613   int i, n2;
00614   double flnm_ret;
00615 
00616   (*flnm_ct)++;
00617   n2 = n + n;
00618   for (i = 0; i < n; i++) v[n+i] = x0[i] + x * v[i];
00619 /*
00620   (*fcn) (n, &v[n], &flnm_ret, &v[n2]);
00621 */
00622   fcn1 (n, &v[n], &flnm_ret, &v[n2], rkx, rky, k2, omeg, ppow, pk, fm, rki, npix);
00623   *df = 0.0;
00624   for (i = 0; i < n; i++) *df += v[i] * v[n2+i];
00625   return flnm_ret;
00626 }
00627 
00628 static int linmin (double *x1, double *x2, double *f1, double *df1,
00629     double reps, double aeps,
00630 /*
00631     void (*f)(int, double *, double *, double *),
00632 */
00633     double *v, double *xi, int n, double *rkx, double *rky,
00634     double *k2, double *omeg, double *ppow, double pk, double *fm, double rki,
00635     int npix, int nit, int *flnm_ct, int *linmin_iter) {
00636   double dfa, df0, df12, df2, dx, dx1, fa, fc, f0, f2, r, r1, xa, x0;
00637   int i, qb, status;
00638 
00639   double rho = 0.01, sigma = 0.1;
00640   double t1 = 9.0, t2 = 0.9, t3 = 0.5;
00641 /*
00642   int nit = 15;
00643 */
00644 
00645   status = 0;
00646   qb = 0;
00647 
00648 /*
00649   f2 = flnm ((*f), *x2, &df2, v, xi, n);
00650 */
00651   f2 = flnm (*x2, &df2, v, xi, n, rkx, rky, k2, omeg, ppow, pk, fm, rki, npix,
00652       flnm_ct);
00653 
00654   dx1 = *x2 - *x1;
00655   f0 = *f1;
00656   df0 = *df1;
00657   x0 = *x1;
00658 
00659   for (i = 0; i < nit; i++) {
00660     *linmin_iter = i + 1;
00661     fc = f0 + df0 * rho * (*x2 - x0);
00662     if ((fabs (df2) <= -sigma * df0) && (f2 <= fc)) {
00663       *x1 = *x2;
00664       *f1 = f2;
00665       *df1 = df2;
00666       return status;
00667     }
00668     if (!qb)
00669       if ((f2 > fc) || (f2 > *f1) || (df2 >= 0.0)) qb = 1;
00670     df12 = (f2 - *f1) / dx1;
00671     r = 2 * df2 + *df1 - 3 * df12;
00672     r1 = 3 * df12 - df2 - *df1;
00673     r1 *= r1;
00674     r1 -= *df1 * df2;
00675     dx = 0.0;
00676     if (r1 > 0.0) {
00677       r1 = sqrt (r1);
00678       if (r < 0.0) r1 *= -1;
00679       r += r1;
00680       dx = -df2 * dx1 / r;
00681     } else {
00682       double r2 = 2 * (df12 - df2);
00683       if (r2 != 0.0) dx = dx1 * df2 / r2;
00684     }
00685 
00686     if (qb) {
00687       if (dx < (-t2 * dx1)) dx = -t2 * dx1;
00688       if (dx > (-t3 * dx1)) dx = -t3 * dx1;
00689       xa = *x2 + dx;
00690 /*
00691       fa = flnm ((*f), xa, &dfa, v, xi, n);
00692 */
00693       fa = flnm (xa, &dfa, v, xi, n, rkx, rky, k2, omeg, ppow, pk, fm, rki,
00694           npix, flnm_ct);
00695       fc = f0 + df0 * rho * (xa - x0);
00696       if ((fabs (dfa) <= (-sigma * df0)) && (fa <= fc)) {
00697         *x1 = xa;
00698     *f1 = fa;
00699     *df1 = dfa;
00700     return status;
00701       } else if ((fa > fc) || (fa >= *f1) || (dfa > 0.0)) {
00702         *x2 = xa;
00703     f2 = fa;
00704     df2 = dfa;
00705       } else {
00706         *x1 = xa;
00707     *f1 = fa;
00708     *df1 = dfa;
00709       }
00710       dx1 = *x2 - *x1;
00711       if ((fabs (dx1) < (reps * fabs (*x2))) || (fabs (dx1) < aeps)) {
00712     status = 135;
00713     if (f2 <= f0) {
00714       *x1 = *x2;
00715       *f1 = f2;
00716       *df1 = df2;
00717       status = 33;
00718     }
00719     return status;
00720       }
00721     } else {
00722       if ((dx < (*x2 - *x1)) && (dx > 0.0)) dx = *x2 - *x1;
00723       if ((dx > t1 * (*x2 - *x1)) || (dx <= 0.0)) dx = t1 * (*x2 - *x1);
00724       *x1 = *x2;
00725       *x2 += dx;
00726       dx1 = dx;
00727       *f1 = f2;
00728       *df1 = df2;
00729 /*
00730       f2 = flnm (f, *x2, &df2, v, xi, n);
00731 */
00732       f2 = flnm (*x2, &df2, v, xi, n, rkx, rky, k2, omeg, ppow, pk, fm, rki,
00733           npix, flnm_ct);
00734     }
00735   }
00736   if (f2 <= f0) {
00737     *f1 = f2;
00738     *x1 = *x2;
00739     *df1 = df2;
00740     return 22;
00741   } else if ((*f1 <= f0) && (*x1 != x0)) {
00742     return 22;
00743   }
00744   return 133;
00745 }
00746 
00747 static int bfgs (int n, double *x, double *f, double *g, double *h,
00748     double reps, double aeps, double *wk,
00749 /*
00750     void (*fcn)(int, double *, double *, double *)
00751 */
00752     double *rkx, double *rky, double *k2, double *omeg, double *ppow, double pk,
00753     double *fm, double rki, int npix, int nit, int linmin_nit, int *flnm_ct,
00754     int *bfgs_iter, int *linmin_iter) {
00755 /*
00756  *  Implements (?) the BFGS (Broyden, Fletcher, Goldfarb, Shanno)
00757  *    algorithm for updating the Hessian matrix
00758  *  The function fcn() is invoked to form the log of the likelihood function
00759  *    f and fill the array g[i] for the parameter set x[n] and the data.
00760  *  The work array wk is internal and could be mallocd to 3n.
00761  */
00762 
00763   double df, dg, f1, ghg, gi, h1, h2, r1, x1, x2;
00764   int i, ier, it, j, n2, qc, status = 0;
00765 /*
00766   int nit = 200;
00767 */
00768 
00769   if (n < 1)
00770     return 111;
00771 /*
00772   (*fcn) (n, x, f, g);
00773 */
00774   fcn1 (n, x, f, g, rkx, rky, k2, omeg, ppow, pk, fm, rki, npix);
00775 
00776   df = fabs (*f);
00777   if (df == 0.0) df = 1.0;
00778   n2 = n + n;
00779   h2 = 1.0;
00780   for (it = 0; it < nit; it++) {
00781     double df1 = 0.0;
00782     *bfgs_iter = it + 1;
00783     h1 = h2;
00784     h2 = 0.0;
00785     for (i = 0; i < n; i++) {
00786       wk[i] = 0.0;
00787       for (j = 0; j < n; j++) {
00788     h2 += fabs (h[i + n*j]);
00789     wk[i] -= h[i + n*j] * g[j];
00790       }
00791       df1 += wk[i] * g[i];
00792     }
00793     if (df1 == 0.0) {
00794       if (fabs (h2 / h1) > 1.3) status = 44;
00795       return status;
00796     }
00797 
00798     x1 = 0.0;
00799     x2 = -2 * df / df1;
00800     if (x2 > 1.0) x2 = 1.0;
00801     f1 = *f;
00802     if (x2 <= 0.0) x2 = 1.0;
00803 /*
00804     linmin (&x1, &x2, &f1, &df1, reps, aeps, &ier, (*fcn),
00805         wk, x, n, linmin_nit);
00806 */
00807     ier = linmin (&x1, &x2, &f1, &df1, reps, aeps, wk, x, n, rkx, rky,
00808         k2, omeg, ppow, pk, fm, rki, npix, linmin_nit, flnm_ct, linmin_iter);
00809     if (ier > 0) status = ier;
00810     if (ier > 100) return status;
00811 
00812     qc = 1;
00813     for (i = 0; i < n; i++) {
00814       x[i] = wk[n+i];
00815       wk[n+i] = x1 * wk[i];
00816       if ((fabs (wk[n+i]) > aeps) &&
00817           (fabs (wk[n+i]) > reps * fabs (x[i]))) qc = 0;
00818       gi = wk[n2+i];
00819       wk[n2+i] -= g[i];
00820       g[i] = gi;
00821     }
00822     df = *f - f1;
00823     *f = f1;
00824     if (qc) {
00825       if (fabs (h2 / h1) > 1.3) status = 44;
00826       return status;
00827     }
00828 
00829     for (i = 0; i < n; i++) {
00830       wk[i] = 0.0;
00831       for (j = 0; j < n; j++) {
00832     wk[i] += h[i + n*j] * wk[n2+j];
00833       }
00834     }
00835     ghg = dg = 0.0;
00836     for (i = 0; i < n; i++) {
00837       dg += wk[n2+i] * wk[n+i];
00838       ghg += wk[i] * wk[n2+i];
00839     }
00840     r1 = (1.0 + (ghg / dg)) / dg;
00841     for (j = 0; j < n; j++)
00842       for (i = 0; i < n; i++)
00843         h[i + n*j] += r1 * wk[n+i] * wk[n+j] -
00844         (wk[n+i] * wk[j] + wk[i] * wk[n+j]) / dg;
00845   }
00846   return 122;
00847 }
00848 
00849 #define RSUNMM      (696.00)
00850 
00851 static int get_scaling_values (DRMS_Record_t *rec, double *dnu, double *dkx,
00852     double *dky) {
00853   double rsun;
00854   double key_dbl;
00855   int status;
00856   char *key_str;
00857   int need_dkx = 0, need_dky = 0, need_dnu = 0;
00858                         /*  do a little sanity checking  */
00859   key_str = drms_getkey_string (rec, "WCSNAME", &status);
00860   if (key_str) {
00861     if (strcasecmp (key_str, "k_x/k_y/omega")) {
00862       fprintf (stderr, "Warning: unknown WCS Type %s\n", key_str);
00863     } else {
00864       key_str = drms_getkey_string (rec, "CTYPE1", &status);
00865       if (strcasecmp (key_str, "WAVE-NUM"))
00866     fprintf (stderr, "Warning: inconsistent value of CTYPE1: %s\n", key_str);
00867       key_str = drms_getkey_string (rec, "CTYPE2", &status);
00868       if (strcasecmp (key_str, "WAVE-NUM"))
00869     fprintf (stderr, "Warning: inconsistent value of CTYPE2: %s\n", key_str);
00870       key_str = drms_getkey_string (rec, "CTYPE3", &status);
00871       if (strcasecmp (key_str, "FREQ-ANG"))
00872     fprintf (stderr, "Warning: inconsistent value of CTYPE3: %s\n", key_str);
00873     }
00874   } else {
00875     fprintf (stderr, "Warning: no WCSTYPE; k_x/k_y/omega assumed\n");
00876   }
00877  /*  should really be able to parse CUNIT and convert accordingly; for now..  */
00878   key_str = drms_getkey_string (rec, "CUNIT1", &status);
00879   if (key_str) {
00880     if (strcasecmp (key_str, "Mm-1")) {
00881       fprintf (stderr, "Warning: unparsed value of CUNIT1: %s\n", key_str);
00882       fprintf (stderr, "  using DELTA_K instead\n");
00883       need_dkx = 1;
00884     } else {
00885       *dkx = drms_getkey_double (rec, "CDELT1", &status);
00886       if (status || isnan (*dkx)) {
00887     fprintf (stderr, "Warning: no valid value for CDELT1\n");
00888     fprintf (stderr, "  using DELTA_K instead\n");
00889     need_dkx = 1;
00890       }
00891     }
00892   }
00893   key_str = drms_getkey_string (rec, "CUNIT2", &status);
00894   if (key_str) {
00895     if (strcasecmp (key_str, "Mm-1")) {
00896       fprintf (stderr, "Warning: unparsed value of CUNIT2: %s\n", key_str);
00897       fprintf (stderr, "  using DELTA_K instead\n");
00898       need_dky = 1;
00899     } else {
00900       *dky = drms_getkey_double (rec, "CDELT2", &status);
00901       if (status || isnan (*dky)) {
00902     fprintf (stderr, "Warning: no valid value for CDELT2\n");
00903     fprintf (stderr, "  using DELTA_K instead\n");
00904     need_dky = 1;
00905       }
00906     }
00907   }
00908   key_str = drms_getkey_string (rec, "CUNIT3", &status);
00909   if (key_str) {
00910     if (strcasecmp (key_str, "s-1")) {
00911       fprintf (stderr, "Warning: unparsed value of CUNIT3: %s\n", key_str);
00912       fprintf (stderr, "  using D_OMEGA or DELTA_NU instead\n");
00913       need_dnu = 1;
00914     } else {
00915       *dnu = drms_getkey_double (rec, "CDELT3", &status);
00916       if (status || isnan (*dnu)) {
00917     fprintf (stderr, "Warning: no valid value for CDELT3\n");
00918     fprintf (stderr, "  using D_OMEGA or DELTA_NU instead\n");
00919     need_dnu = 1;
00920       }
00921     }
00922   }
00923   
00924   if (need_dkx || need_dky) {
00925     key_dbl = drms_getkey_double (rec, "DELTA_K", &status);
00926     if (status || isnan (key_dbl)) {
00927       fprintf (stderr, "Error: no valid value for delta-k\n");
00928       return 1;
00929     }
00930     if (need_dkx) *dkx = key_dbl;
00931     if (need_dky) *dky = key_dbl;
00932   }
00933   rsun = drms_getkey_double (rec, "RSunRef", &status);
00934   if (status || isnan (rsun)) {
00935     rsun = RSUNMM;
00936     fprintf (stderr, "Warning: no valid value for RSunRef; ");
00937     fprintf (stderr, "using %f\n", rsun);
00938   }
00939   *dkx *= rsun;
00940   *dky *= rsun;
00941 
00942   if (need_dnu) {
00943     *dnu = drms_getkey_double (rec, "DELTA_NU", &status);
00944     if (status || isnan (*dnu)) {
00945       *dnu = drms_getkey_double (rec, "D_OMEGA", &status);
00946       if (status || isnan (*dnu)) {
00947     fprintf (stderr, "Error: no valid value for delta-nu\n");
00948     return 1;
00949       } else *dnu *= 0.5e6 / M_PI;
00950     }
00951   } else *dnu *= 0.5e6 / M_PI;
00952 
00953   return 0;
00954 }
00955 
00956 int DoIt (void) {
00957   CmdParams_t *params = &cmdparams;
00958   DRMS_RecordSet_t *ids;
00959   DRMS_Record_t *irec, *orec;
00960   DRMS_Segment_t *iseg, *oseg;
00961   DRMS_Array_t *pspec;
00962   FILE *unit11, *unit22, *runlog;
00963   double *rkl, *rk0, *rku;
00964   double *freq, *wid;
00965   double *frl, *rk;
00966   double *rkk, *rkx, *rky, *k2, *omeg, *ppow;
00967   double *d2, *h, *hs;
00968   double par[20], wk[900], par0[20];
00969   double g[20];
00970   double aeps, amp, back, bk1, bk2, det, dkx, dky, dnu;
00971   double f, fr, f0, f00;
00972   double pk1, reps, kmax, rkw, rl0, sn, wd;
00973   double pk, fm, rki, xf;
00974   double dval, kval;
00975   double *l_fit, *freq_fit, *ux_fit, *uy_fit, **par_fit, *rki_fit, *f_fit;
00976   double *fm_fit, **h_fit, **hs_fit;
00977   float *spec;
00978   int intg[20];
00979   int c;
00980   int i, j, jj, k, kx, k0, l;
00981   int col, row, pln;
00982   int bfgs_iter, bfgs_status, flnm_ct, linmin_iter;
00983   int iers, fp0, fplo, fpup, fp;
00984   int nkx, nky, nnu, nfr, nstp;
00985   int n, nk, nt, nx2, ny2;
00986   int fit_ct, cvg_ct, cvg_minf, cvg_maxf;
00987   int npx, npxmax, status, total_clock;
00988   int rgn, rgnct, segct, isegnum, osegnum, logsegnum, drms_output, dispose;
00989   int propct;
00990   int pspec_restrict;
00991   int *n_fit, *mask, *bfgs_fit;
00992   char **copykeylist;
00993   char logfile[DRMS_MAXPATHLEN], outfile[DRMS_MAXPATHLEN];
00994   char source[DRMS_MAXQUERYLEN], recid[DRMS_MAXQUERYLEN];
00995   char line[1024], module_ident[64], key[64];
00996 
00997   double scale[] = {1.0, 1.0, VSCALE, VSCALE, 1.0,
00998       1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0e-4};
00999   double rscale  = 1.0 / RSUN_MM;
01000   float dofit[13];
01001   int npar = 13;
01002   char *hdrstr[] = {
01003     "          A0", "           c", "          Ux", "          Uy",
01004     "          w0", "          B1", "          B2", "          A1",
01005     "          A2", "          A3", "           p", "          w1",
01006     "           S"};
01007   char *hdrstr2[] = {
01008     "        D_A0", "         D_c", "        D_Ux", "        D_Uy",
01009     "        D_w0", "        D_B1", "        D_B2", "        D_A1",
01010     "        D_A2", "        D_A3", "         D_p", "        D_w1",
01011     "         D_S"};
01012 
01013   int nmin = params_get_int (params, "nmin");
01014   int nmax = params_get_int (params, "nmax");
01015                       /*  should actually be nmax - nmin + 1  */
01016   int nct = nmax + 2;
01017   int lmin = params_get_int (params, "lmin");
01018   int lmax = params_get_int (params, "lmax");
01019                       /*  should actually be lmax - lmin + 1  */
01020   int lct = lmax + 1;
01021   int fmin = params_get_int (params, "fmin");
01022   int fmax = params_get_int (params, "fmax");
01023   int bfgs_nit = params_get_int (params, "bfgsct");
01024   int linmin_nit = params_get_int (params, "linminct");
01025 
01026   double A1_val = params_get_double (params, "A1");
01027   double B1_val = params_get_double (params, "B1");
01028   double B2_val = params_get_double (params, "B2");
01029   double asym_val = params_get_double (params, "S");
01030 
01031   double ux_guess = params_get_double (params, "ux") / VSCALE;
01032   double uy_guess = params_get_double (params, "uy") / VSCALE;
01033   double A1_guess = params_get_double (params, "A1_guess");
01034   double asym_guess = params_get_double (params, "S_guess");
01035 
01036   double kxmin = params_get_double (params, "kxmin");
01037   double kxmax = params_get_double (params, "kxmax");
01038   double kymin = params_get_double (params, "kymin");
01039   double kymax = params_get_double (params, "kymax");
01040 
01041   char *inds = strdup (params_get_str (params, "in"));
01042   char *guessfile =  strdup (params_get_str (params, "guessfile"));
01043   char *outser =  strdup (params_get_str (params, "out"));
01044   char *propagate_req = strdup (params_get_str (params, "copy"));
01045 
01046   int verbose = params_isflagset (params, "v");
01047   int extra_reporting = params_isflagset (params, "x");
01048   int calc_seconds = params_isflagset (params, "2");
01049   int no_fits = params_isflagset (params, "n");
01050  
01051   double p = 0.0;                  /*  uninitialized in original  */
01052 
01053 /*
01054   time_elapsed (0, &total_clock);
01055 */
01056   snprintf (module_ident, 64, "%s v %s", module_name, version_id);
01057   if (verbose) printf ("%s:\n", module_ident);
01058   propct = construct_stringlist (propagate_req, ',', &copykeylist);
01059   if ((unit11 = fopen (guessfile, "r")) == NULL) {
01060     fprintf (stderr, "Error: unable to read frequency guess table \"%s\"\n",
01061     guessfile);
01062     return 1;
01063   }
01064                 /*  optionally eliminate parameters from fit  */
01065   for (n = 0; n < npar; n++) dofit[n] = 1.0;
01066 /*
01067   if (isfinite (B1_val)) dofit[5] = 0.0;
01068   if (isfinite (B2_val)) dofit[6] = 0.0;
01069 */
01070   if (isfinite (A1_val)) {
01071     dofit[7] = 0.0;
01072     A1_guess = A1_val;
01073   }
01074 /*
01075   if (isfinite (A2_val)) dofit[8] = 0.0;
01076   if (isfinite (A3_val)) dofit[9] = 0.0;
01077 */
01078   if (isfinite (asym_val)) {
01079     dofit[12] = 0.0;
01080     asym_guess = asym_val;
01081   }
01082   pspec_restrict = (isfinite (kxmin) || isfinite (kxmax) || isfinite (kymin) ||
01083       isfinite (kymax));
01084 
01085   d2 = (double *)malloc (npar * npar * sizeof (double));
01086   h = (double *)malloc (npar * npar * sizeof (double));
01087   hs = (double *)malloc (npar * npar * sizeof (double));
01088   freq = (double *)calloc (nct * lct, sizeof (double));
01089   wid = (double *)calloc (nct * lct, sizeof (double));
01090           /*  Read in the initial guess for frequency and line-width  */
01091   i = j = 0;
01092       /*  Read frequencies and widths from table into [nmax*lmax] arrays  */
01093   while ((c = fgetc (unit11)) != EOF) {
01094     if (c == '\n') {
01095       line[j] = '\0';
01096       if (sscanf (line, "%d %d %lf %lf", &n, &l, &fr, &wd) == 4) {
01097         if ((n <= nmax + 1) && (n >= nmin) && (l <= lmax) && (l >= lmin)) {
01098       if (wd > 600.0) wd = 600.0;
01099 /* this should probably be removed  */
01100       if (fr <= freq[n + nct*(l-1)])
01101         fr = freq[n + nct*(l-1)] + 0.5;
01102 
01103       freq[n + l*nct] = fr;
01104       wid[n + l*nct] = 0.5 * wd;
01105     }
01106       }
01107       j = 0;
01108       i++;
01109     } else line[j++] = (char)c;
01110   }
01111   fclose (unit11);
01112   if (verbose) printf ("%d entries read from guess table\n", i);
01113 
01114   ids = drms_open_records (drms_env, inds, &status);
01115   if (status) {
01116     fprintf (stderr, "Error: drms_open_records returned %d for dataset:\n", status);
01117     fprintf (stderr, "  %s\n", inds);
01118     return 0;
01119   }
01120              /*  code originally written to process single spectrum 
01121               processing of multiple regions should be parallelized  */
01122   rgnct = ids->n;
01123   if (rgnct < 1) {
01124     fprintf (stderr, "Error: no records found in requested dataset:\n");
01125     fprintf (stderr, "  %s\n", inds);
01126     return 0;
01127   }
01128   irec = ids->records[0];
01129   segct = drms_record_numsegments (irec);
01130   if (segct != 1) {
01131     int found = 0;
01132     for (n = 0; n < segct; n++) {
01133       iseg = drms_segment_lookupnum (irec, n);
01134       if (!iseg) continue;
01135       if (iseg->info->naxis == 3) {
01136         if (!found) isegnum = n;
01137     found++;
01138       }
01139     }
01140     if (!found) {
01141       fprintf (stderr, "Error: found no segment of rank 3 in input dataset\n");
01142       drms_close_records (ids, DRMS_FREE_RECORD);
01143       return 1;
01144     }
01145     if (found > 1) {
01146       fprintf (stderr,
01147       "Warning: found multiple segments of rank 3 in input dataset\n");
01148       iseg = drms_segment_lookupnum (irec, isegnum);
01149       fprintf (stderr, "         using segment %s\n", iseg->info->name);
01150     }
01151   } else {
01152     segct = DRMS_MAXSEGMENTS;
01153     for (n = 0; n < segct; n++) {
01154       iseg = drms_segment_lookupnum (irec, n);
01155       if (!iseg) continue;
01156       if (iseg->info->naxis != 3) {
01157         fprintf (stderr, "Error: found no segment of rank 3 in input dataset\n");
01158         drms_close_records (ids, DRMS_FREE_RECORD);
01159         return 1;
01160       }
01161       isegnum = n;
01162       break;
01163     }
01164   }
01165                     /*  check output data series struct  */
01166   orec = drms_create_record (drms_env, outser, DRMS_TRANSIENT, &status);
01167   if (orec) {
01168     segct = drms_record_numsegments (orec);
01169     if (segct < 1) {
01170       fprintf (stderr,
01171       "Error: no data segment in output series %s\n", outser);
01172       drms_close_records (ids, DRMS_FREE_RECORD);
01173       drms_close_record (orec, DRMS_FREE_RECORD);
01174       return 1;
01175     }
01176     oseg = drms_segment_lookup (orec, "fit.out");
01177     if (oseg && oseg->info->protocol == DRMS_GENERIC)
01178       osegnum = oseg->info->segnum;
01179     else {
01180     /*  either segment fit.out does not exist, or it has the wrong protocol  */
01181       for (n = 0; n < segct; n++) {
01182     oseg = drms_segment_lookupnum (orec, n);
01183     if (oseg->info->protocol != DRMS_GENERIC) continue;
01184     osegnum = n;
01185     break;
01186       }
01187       if (n == segct) {
01188     fprintf (stderr,
01189         "Error: no generic data segment in output series %s\n", outser);
01190     drms_close_records (ids, DRMS_FREE_RECORD);
01191     return 1;
01192       }
01193       fprintf (stderr, "         writing to segment %s\n", oseg->info->name);
01194     }
01195     logsegnum = -1;
01196     oseg = drms_segment_lookup (orec, "Log");
01197     if (oseg && oseg->info->protocol == DRMS_GENERIC)
01198       logsegnum = oseg->info->segnum;
01199     drms_close_record (orec, DRMS_FREE_RECORD);
01200     drms_output = 1;
01201   } else {
01202        /*  Can't create in named series, assume it's a filename instead  */
01203     if (rgnct > 1) {
01204       fprintf (stderr,
01205       "Query produced %d matching records with output to single file;",
01206       rgnct);
01207       fprintf (stderr, " must limit to 1\n");
01208       drms_close_records (ids, DRMS_FREE_RECORD);
01209       return 1;
01210     }
01211     strcpy (outfile, outser);
01212     drms_output = 0;
01213   }
01214   dispose = (no_fits) ? DRMS_FREE_RECORD : DRMS_INSERT_RECORD;
01215     /*  loop over all input records, creating an output record for each  */
01216   for (rgn = 0; rgn < rgnct; rgn++) {
01217     time_elapsed (0, &total_clock);
01218     irec = ids->records[rgn];
01219     drms_sprint_rec_query (source, irec);
01220     iseg = drms_segment_lookupnum (irec, isegnum);
01221     pspec = drms_segment_read (iseg, DRMS_TYPE_FLOAT, &status);
01222     if (status) {
01223       fprintf (stderr, "Warning: could not read segment %d\n", isegnum);
01224       fprintf (stderr, "       in %s; skipped\n", source);
01225       continue;
01226     }
01227     if (drms_output) {
01228       char *suffix;
01229       orec = drms_create_record (drms_env, outser, DRMS_PERMANENT, &status);
01230       oseg = drms_segment_lookupnum (orec, osegnum);
01231       drms_segment_filename (oseg, outfile);
01232                   /*  strip off the .generic tag if present  */
01233       suffix = strstr (outfile, ".generic");
01234       if (suffix) *suffix = '\0';
01235       if (logsegnum >= 0) {
01236         oseg = drms_segment_lookupnum (orec, logsegnum);
01237         drms_segment_filename (oseg, logfile);
01238                   /*  strip off the .generic tag if present  */
01239         suffix = strstr (logfile, ".generic");
01240         if (suffix) *suffix = '\0';
01241     runlog = fopen (logfile, "w");
01242       }
01243     } else runlog = stderr;
01244     unit22 = fopen (outfile, "w");
01245     if (!unit22) {
01246       fprintf (stderr, "Error: could not open output file %s\n", outfile);
01247       return 1;
01248     }
01249 
01250     nkx = pspec->axis[0];
01251     nky = pspec->axis[1];
01252     nnu = pspec->axis[2];
01253           /*  Warning: this "22" must be at least 2 larger than nct  */
01254     rkl = (double *)malloc (22 * nnu * sizeof (double));
01255     rk0 = (double *)malloc (22 * nnu * sizeof (double));
01256     rku = (double *)malloc (22 * nnu * sizeof (double));
01257 
01258     if (get_scaling_values (irec, &dnu, &dkx, &dky)) {
01259       fprintf (stderr, "Error: region %d lacks appropriate scaling keywords\n",
01260       rgn);
01261       continue;
01262     }
01263     nstp = (int)(rint (8.0 / dnu + 0.5) + 0.1);
01264     if (nstp < 1) nstp = 1;
01265     nfr = (int)(rint (100.0 / dnu + 0.5) + 0.1);
01266     if (nfr < 10) nfr = 10;
01267 /*
01268     kscale = drms_getkey_double (irec, "MAPSCALE", &status);
01269     if (status || isnan (kscale)) {
01270       fprintf (stderr, "Warning: required keyword MAPSCALE missing: assumed 0.125\n");
01271       kscale = 0.125;
01272     }
01273     dlon = kscale * M_PI / 180.0;
01274     if (verbose) printf ("original image scale = %.8f rad/pixel\n", dlon); 
01275 */
01276                            /*  Assumes azimuthal mapping  */
01277 /*
01278     keyval = drms_getkey_string (irec, "MapProj", &status);
01279     if (keyval) {
01280       char *projection = (char *)keyval;
01281       if (strncmp (projection, "Postel", 6)) {
01282     fprintf (stderr,
01283         "Warning: unknown projection %s; azimuthal symmetry assumed\n",
01284         projection);
01285       }
01286     dlat = dlon;
01287     } else {
01288       fprintf (stderr,
01289       "Warning: unknown projection; azimuthal symmetry assumed\n");
01290       dlat = dlon;
01291     }
01292 */
01293     spec = (float *)pspec->data;
01294     dval = drms_getkey_double (irec, "LOG_BASE", &status);
01295     if (!status && isfinite (dval)) {
01296       double scaled_log_base = log (dval);
01297       int ntot = nkx * nky * nnu;
01298       for (n = 0; n < ntot; n++) spec[n] = (float)exp (scaled_log_base * spec[n]);
01299     }
01300             /*  if kx or ky restricted, zero out relevant values  */
01301     if (pspec_restrict) {
01302       if (isfinite (kymin)) {
01303     for (n = 0, pln = 0; pln < nnu; pln++) {
01304       for (row = 0; row < nky; row++) {
01305         kval = dky * (row - 0.5 * nky);
01306         if (kval > kymin) {
01307           row = nky;
01308           continue;
01309         }
01310         for (col = 0; col < nkx; col++, n++) spec[n] = 0.0;
01311       }
01312     }
01313       }
01314       if (isfinite (kymax)) {
01315     for (n = 0, pln = 0; pln < nnu; pln++) {
01316       for (row = 0; row < nky; row++) {
01317         kval = dky * (row - 0.5 * nky);
01318         if (kval <= kymax) {
01319           n += nkx;
01320           continue;
01321         }
01322         for (col = 0; col < nkx; col++, n++) spec[n] = 0.0;
01323       }
01324     }
01325       }
01326       if (isfinite (kxmin)) {
01327     for (n = 0, pln = 0; pln < nnu; pln++) {
01328       for (row = 0; row < nky; row++) {
01329         for (col = 0; col < nkx; col++, n++) {
01330           kval = dkx * (col - 0.5 * nkx);
01331           if (kval < kxmin) spec[n] = 0.0;
01332         }
01333       }
01334     }
01335       }
01336       if (isfinite (kxmax)) {
01337     for (n = 0, pln = 0; pln < nnu; pln++) {
01338       for (row = 0; row < nky; row++) {
01339         for (col = 0; col < nkx; col++, n++) {
01340           kval = dkx * (col - 0.5 * nkx);
01341           if (kval > kxmax) spec[n] = 0.0;
01342         }
01343       }
01344     }
01345       }
01346    }
01347 
01348     nt = nnu;
01349 /*
01350     dkx = 2.0 * M_PI / (dlon * nkx);
01351     dky = 2.0 * M_PI / (dlat * nky);
01352 */
01353             /*  dk is /rad (360 deg / map width), effectively dl  */
01354     kmax = dkx * (nkx/2 - 3);
01355     if (verbose) printf ("resolution: dkx: %.2f, dky: %.2f, dnu: %f, kmax = %.2f\n",
01356     dkx, dky, dnu, kmax);
01357     fprintf (unit22, "# fit of %s\n", inds);
01358     fprintf (unit22, "# %s version %s\n", module_name, version_id);
01359     fprintf (unit22, "# nmin=%d nmax=%d lmin=%d lmax=%d fmin=%d fmax=%d\n",
01360     nmin, nmax, lmin, lmax, fmin, fmax);
01361     fprintf (unit22, "# bfgsct=%d linminct=%d\n", bfgs_nit, linmin_nit);
01362     fprintf (unit22, "#n        l        k        nu    D_w0");
01363     for (n = 2; n < 4; n++) {
01364       fprintf (unit22, "%s%s", hdrstr[n], hdrstr2[n]);
01365       if (calc_seconds) fprintf (unit22, "            ");
01366     }
01367     fprintf (unit22, " qual");
01368     fprintf (unit22, "  l_guess bfgs           f          fm");
01369     for (n = 0; n < 2; n++) {
01370       fprintf (unit22, "%s%s", hdrstr[n], hdrstr2[n]);
01371       if (calc_seconds) fprintf (unit22, "            ");
01372     }
01373     for (n = 4; n < npar; n++) {
01374       fprintf (unit22, "%s%s", hdrstr[n], hdrstr2[n]);
01375       if (calc_seconds) fprintf (unit22, "             ");
01376     }
01377 /*
01378     fprintf (unit22, "#n        nu  l_guess        l bfgs");
01379     fprintf (unit22, "           f          fm");
01380     for (n = 0; n < npar; n++) fprintf (unit22, "%s", hdrstr[n]);
01381     if (calc_seconds)
01382       for (n = 0; n < npar; n++) fprintf (unit22, "%s            ", hdrstr2[n]);
01383     else
01384       for (n = 0; n < npar; n++) fprintf (unit22, "%s", hdrstr2[n]);
01385 */
01386 /*
01387     fprintf (unit22, "         pk          sn");
01388     if (extra_reporting) fprintf (unit22, " flnm bfgs linm");
01389 */
01390     fprintf (unit22, "\n");
01391 
01392     fplo = (int)(fmin / dnu);
01393     if (fplo <= nfr) fplo = nfr + 1;
01394     fpup = (int)(fmax / dnu);
01395     if (fpup > (nt - nfr)) fpup = nt - nfr;
01396     if (verbose) printf ("fplo = %d fpup = %d\n", fplo, fpup);
01397 
01398     reps = 1.0e-6;
01399     frl = (double *)malloc (lct * sizeof (double));
01400     rk = (double *)malloc (lct * sizeof (double));
01401 /*
01402  *  Interpolate l-values from guess table to target frequencies
01403  *    for each n, set up [lct] arays of l and f values
01404  *    interpolated values stored in [22 * nnu] array rk0
01405  *    based on ??, in arrays rkl and rku
01406  *
01407  *  Frequencies are assumed monotone increasing with l for given n
01408  */
01409     for (fp = fplo - nfr; fp <= fpup + nfr; fp++) {
01410       f0 = fp * dnu;
01411       for (n = 0; n <= nmax; n++) {
01412     nk = 0;
01413     for (l = 0; l < lct; l++) {
01414       if (freq[n + nct*l] > 0.0) {
01415         frl[nk] = freq[n + nct*l];
01416         rk[nk] = l;
01417         nk++;
01418       }
01419     }
01420     for (i = 0; i < nk; i++) if (frl[i] > f0) break;
01421     if (i >= nk) rk0[fp + nnu*n] = lct + n - 1;
01422     else {
01423       xf = (frl[i] - f0) / (frl[i] - frl[i-1]);
01424       rk0[fp + nnu*n] = rk[i] - xf * (rk[i] - rk[i-1]);
01425     }
01426     if (f0 > frl[nk-1]) {
01427     /* for frequencies beyond tabulated range in guess table, set to max  */
01428       rk0[fp + nnu*n] = lct + n - 1;
01429     }
01430       }
01431 /*
01432  *  Set target bounds in [22 * nnu] arrays rkl and rku
01433  *  For n > 0 ridges, boundaries are 1/3 of distance to adjoining ridges
01434  *  For n = 0, lower bound is as above, upper bound is equidistant from target
01435  */
01436       rkl[fp] = (2 * rk0[fp] + rk0[fp + nnu]) / 3.0;
01437       rku[fp] = rk0[fp] + (rk0[fp] - rkl[fp]);
01438       if (rku[fp] > kmax) rku[fp] = kmax;
01439       for (n = 1; n < 18; n++) {
01440     rku[fp + nnu*n] = (2 * rk0[fp + nnu*n] + rk0[fp + nnu*(n-1)]) / 3.0;
01441     if (rku[fp + nnu*n] > kmax) rku[fp + nnu*n] = kmax;
01442     rkl[fp + nnu*n] = (2 * rk0[fp + nnu*n] + rk0[fp + nnu*(n+1)]) / 3.0;
01443       }
01444       if (verbose) {
01445     printf ("%8.3f", f0);
01446     for (j = nmin; j <= nmax; j++) printf ("%11.4e", rk0[fp + nnu*j]);
01447     printf ("\n");
01448       }
01449     }
01450     free (frl);
01451     free (rk);
01452 
01453     nx2 = nkx / 2 + 1;
01454     ny2 = nky / 2 + 1;
01455     npxmax = nkx * nky * (2 * nfr + 1);
01456     rkk = (double *)malloc (npxmax * sizeof (double));
01457     rkx = (double *)malloc (npxmax * sizeof (double));
01458     rky = (double *)malloc (npxmax * sizeof (double));
01459     k2 = (double *)malloc (npxmax * sizeof (double));
01460     omeg = (double *)malloc (npxmax * sizeof (double));
01461     ppow = (double *)malloc (npxmax * sizeof (double));
01462 
01463                     /* added 03/24/2011 */
01464     n_fit = (int*) malloc(lct*sizeof(int));
01465     l_fit = (double*) malloc(lct*sizeof(double));
01466     freq_fit = (double*) malloc(lct*sizeof(double));
01467     ux_fit = (double*) malloc(lct*sizeof(double));
01468     uy_fit = (double*) malloc(lct*sizeof(double));
01469     par_fit = (double**) malloc(lct*sizeof(double*));
01470     rki_fit = (double*) malloc(lct*sizeof(double));
01471     bfgs_fit = (int*) malloc(lct*sizeof(int));
01472     f_fit = (double*) malloc(lct*sizeof(double));
01473     fm_fit = (double*) malloc(lct*sizeof(double));
01474     mask = (int*) malloc(lct*sizeof(int));
01475     h_fit = (double**) malloc(lct*sizeof(double*));
01476     hs_fit = (double**) malloc(lct*sizeof(double*));
01477     for(i=0; i<lct; i++)    {
01478       par_fit[i] = (double*) malloc(npar*sizeof(double));
01479       h_fit[i] = (double*) malloc(npar*npar*sizeof(double));
01480       hs_fit[i] = (double*) malloc(npar*npar*sizeof(double));
01481     }
01482 
01483     for (n = nmin; n <= nmax; n++) {
01484       int row = nnu * n;
01485       int reject_ux = 0, reject_bfgs =0, reject_k = 0, reject_fm = 0,
01486           reject_sn = 0, reject_l = 0, reject_tot = 0;
01487 
01488       pk = (n) ? 0.75 : 0.5;
01489       iers = 200;
01490       f00 = 6000.0;
01491       cvg_ct = 0;
01492       cvg_minf = fpup;
01493       cvg_maxf = fplo;
01494 
01495       if (no_fits) {
01496     int expfct = fpup - fplo + 1;
01497     printf
01498         ("n = %2d, row = %d, flo = %.2f fhi = %.2f fct = %d, k0 = %.2f, ku = %.2f\n",
01499         n, row, fplo * dnu, fpup * dnu, expfct,  rk0[fp0 + row],
01500         rku[fp0 + row]);
01501       }
01502       for (fp0 = fplo, fit_ct = 0; fp0 <= fpup; fp0 += nstp, fit_ct++) {
01503     double powmax = 0.0;
01504 
01505     if (no_fits) continue;
01506     if (iers == 300) continue;
01507                                   /*  Changed 2004.03.30  */
01508     f0 = fp0 * dnu;
01509     k0 = (int)rk0[fp0 + row];
01510     rki = rk0[fp0 + row];
01511     kx = (int)(rk0[fp0 + row] / dkx) + nx2;
01512     rkw = rk0[fp0 + row] - rkl[fp0 + row];
01513     if (rkw > (rku[fp0 + row] - rk0[fp0 + row]))
01514           rkw = rku[fp0 + row] - rk0[fp0 + row];
01515     if (verbose) {
01516           printf ("region %.3f %.3f %.3f %.3f\n", rkl[fp0 + row],
01517           rku[fp0 + row], f0, rk0[fp0 + row]);
01518           printf ("%d %d\n", kx, fp0);
01519           printf ("%.4f %.4f %.4f %.4f\n",
01520               rkl[fp0 + row] / dkx + nx2, rku[fp0 + row] / dkx + nx2,
01521           nx2 - rkl[fp0 + row] / dkx, nx2 - rku[fp0 + row] / dkx);
01522     }
01523     if (kx > (nkx - 3)) {
01524       if (verbose) printf ("kx = %d > %d; skipped\n", kx, nkx - 3);
01525       continue;
01526     }
01527     if ((k0 < lmin) || (k0 > lmax)) {
01528       if (verbose) printf ("k0 = %d outside [%d, %d]; skipped\n", k0,
01529           lmin, lmax);
01530       continue;
01531     }
01532     if (freq[n + nct*(k0-1)] > 0.0)
01533           p = (freq[n + nct*k0] - freq[n + nct*(k0-1)]) * k0 /
01534           freq[n + nct*(k0-1)];
01535     if (verbose) printf ("p = %f\n", p);
01536     if ((p > 0.3) && (p < 0.5) && (k0 < 495)) pk = p;
01537     npx = 0;
01538 /*
01539       oml = f0 - 3.6 * dnu;
01540       omu = f0 + 3.6 * dnu;
01541 */
01542     for (i = fp0 - nfr; i <= fp0 + nfr; i++) {
01543                                   /*  Changed 2007.08.10  */
01544       double om = i * dnu;
01545       double klow = rkl[i + row];
01546       double khigh = rku[i + row];
01547       for (j = 0; j < nkx; j++) {
01548         double rlx = (j + 1 - nx2) * dkx;
01549         double rlx2 = rlx * rlx;
01550                                   /*  Changed 2007.08.10  */
01551         int powrow = j + nkx*nky*i;
01552         for (k = 0; k < nky; k++) {
01553           double rly = (k + 1 - ny2) * dky;
01554           double kk = sqrt (rlx2 + rly * rly);
01555           if ((kk > klow) && (kk < khigh)) {
01556         rkk[npx] = kk;
01557         rkx[npx] = rlx;
01558         rky[npx] = rly;
01559         k2[npx] = rlx2 + rly * rly;
01560         omeg[npx] = om;
01561         ppow[npx] = spec[powrow + nky*k];
01562         if (ppow[npx] > powmax) powmax = ppow[npx];
01563         npx++;
01564           }
01565         }
01566       }
01567     }
01568 
01569     if ((npx < 20) || (powmax <= 0.0)) {
01570       if (verbose)
01571         printf ("npx = %d < 20 || powmax = %.1f <= 0.0; skipped\n", npx,
01572         powmax);
01573       continue;
01574     }
01575                             /*  The initial guesses  */
01576 /*  Parameter correspondences (notation of Basu et al. ApJ 512, 458 and
01577  *    Basu & Antia, Solar Phys 192, 469):
01578  *    par[0]    A_0     7   A_1
01579  *  1   c       8   A_2
01580  *  2   U_x     9   A_3
01581  *  3   U_y     10  p
01582  *  4   w_0     11  w_1
01583  *  5   B_1     12  S
01584  *  6   B_2
01585  */
01586     if ((fp0 == fplo) || (iers > 100)) {
01587       par[1] = f0 / exp (pk * log (rk0[fp0 + row]));
01588       par[2] = ux_guess;
01589           par[3] = uy_guess;
01590       par[7] = A1_guess;
01591 /*
01592  * change to NaN?
01593  */
01594       par[0] = log (powmax);
01595       par[4] = wid[n + nct*k0];
01596       if (par[4] < dnu) par[4] = dnu;
01597       if (par[4] > 100.0) par[4] = 100.0;
01598 /*
01599     par[6] = par[1] * pow (rk0 / rkp, pk);
01600 */
01601       par[10] = par[11] = 0.0;
01602       par[8] = 10.0 * (rki - 200.0);
01603       if (par[8] < 0.0) par[8] = 0.0;
01604       par[9] = 6.0 * (rki - 300.0);
01605       if (par[9] < 0.0) par[9] = 0.0;
01606                       /*  par[12] = 0 implies symmetric fit  */
01607       par[12] = asym_guess;
01608       par[5] = log (powmax) + 3 * log (rk0[fp0 + row]) - 5.0;
01609       par[6] = par[5] + log (rk0[fp0 + row]);
01610     } else {
01611       for (jj = 0; jj < npar; jj++) par[jj] = par0[jj];
01612       back = par[6] - log (rki);
01613       if (back < par[5]) back = par[5];
01614       par[5] = back;
01615       par[6] = back + log (rki);
01616 /*
01617     par[12] = 500.0;
01618 */
01619     }
01620     reps = 1.e-5;
01621     aeps = 1.e-7;
01622 
01623     for (i = 0; i < npar; i++) {
01624       for (j = 0; j < npar; j++) h[j + npar*i] = 0.0;
01625       h[i + npar*i] = dofit[i];
01626     }
01627     if (verbose)
01628           printf ("%d %f %f %f %f\n", npx, par[1], par[4], par[5], par[6]);
01629     flnm_ct = 0;
01630     bfgs_status = bfgs (npar, par, &f, g, h, reps, aeps, wk,
01631 /*
01632           (*fcn1),
01633 */
01634             rkx, rky, k2, omeg, ppow, pk, &fm, rki, npx, bfgs_nit, linmin_nit,
01635         &flnm_ct, &bfgs_iter, &linmin_iter);
01636     if (verbose) {
01637           printf ("%5d %5d %12.4e\n", bfgs_status, flnm_ct, f);
01638           for (i = 0; i < npar; i++) {
01639             printf ("%12.4e", par[i]);
01640         if (i%5 == 4) printf ("\n");
01641           }
01642           printf ("\n");
01643           for (i = 0; i < npar; i++) printf
01644           ("%12.4e", sqrt (scale[i] * scale[i] * fabs (h[i + npar*i])));
01645           printf ("\n");
01646     }
01647     fcn1 (npar, par, &f, g, rkx, rky, k2, omeg, ppow, pk, &fm, rki, npx);
01648     pk1 = pk + 1.0e-4 * par[10];
01649     fm /= npx - npar;
01650     rl0 = f0 / par[1];
01651     if (rl0 > 0.0) rl0 = exp (log (rl0) / pk1);
01652     amp = par[0];
01653     bk1 = par[5] - 3 * log (fabs (rl0));
01654     bk2 = par[6] - 4 * log (fabs (rl0));
01655     sn = amp - ((bk1 > bk2) ? bk1 : bk2);
01656     if (verbose) printf ("merit %12.4e %.4f\n", fm, sn);
01657 
01658     if (calc_seconds) {
01659           d2fcn (npar, par, &f, g, d2, rkk, rkx, rky, k2, omeg, ppow, pk, rki,
01660           npx);
01661           for (i = 0; i < npar; i++) {
01662         for (j = 0; j < npar; j++) hs[j + npar*i] = 0.0;
01663         hs[i + npar*i] = 1.0;
01664           }
01665       if (gauelm (npar, npar, d2, hs, &det, intg, npar, 0) > 0) return 111;
01666           if (verbose) {
01667             for (i = 0; i < npar; i++) printf ("%12.5e",
01668             sqrt (scale[i] * scale[i] * fabs (hs[i + npar*i])));
01669             printf ("\n");
01670           }
01671     }
01672 
01673     if ((f0 - f00) > 120) iers = 200;
01674     if (((f0 - f00) > 220) && (f0 > 3500) && (n <= 1)) iers = 300;
01675     if (((f0 - f00) > 300) && (f0 > 4000) && (n >= 2)) iers = 300;
01676 /*
01677     if ((fabs (par[2]) < 500.0) && (bfgs_status < 100) &&
01678 */
01679         if ((bfgs_status < 100) &&
01680             (fabs (rl0 - rki) < 50.0) && (fm < 26.5) && (sn > -2.0) &&
01681         (sn < 45.0) && (rl0 > 40.0)) {
01682       n_fit[cvg_ct] = n;
01683       l_fit[cvg_ct] = rl0;
01684       freq_fit[cvg_ct] = f0;
01685       ux_fit[cvg_ct] = par[2];
01686       uy_fit[cvg_ct] = par[3];
01687       rki_fit[cvg_ct] = rki;
01688       bfgs_fit[cvg_ct] = bfgs_status;
01689       f_fit[cvg_ct] = f;
01690       fm_fit[cvg_ct] = fm;
01691       for(i=0; i<npar; i++) {
01692         par_fit[cvg_ct][i] = par[i];
01693         for(j=0; j<npar; j++) {
01694           h_fit[cvg_ct][i+j*npar] = h[i+j*npar];
01695           hs_fit[cvg_ct][i+j*npar] = hs[i+j*npar];
01696         }
01697       }
01698       if (!cvg_ct) cvg_minf = fp0;
01699       cvg_maxf = fp0;
01700       cvg_ct++;
01701 /*
01702     fprintf (unit22, "#n         l      k       nu      w0");
01703     fprintf (unit22, "           f          fm");
01704     for (n = 2; n < 4; n++)
01705       fprintf (unit22, "%s  D_%s", hdrstr[n], hdrstr[n]);
01706     fprintf (unit22, " l_guess bfgs");
01707     for (n = 0; n < 2; n++)
01708       fprintf (unit22, "%s  D_%s", hdrstr[n], hdrstr[n]);
01709     for (n = 4; n < npar; n++) {
01710       fprintf (unit22, "%s  D_%s", hdrstr[n], hdrstr[n]);
01711       if (calc_seconds) fprintf (unit22, "            ");
01712     }
01713 */
01714       /* removed 3/24/2011 */
01715 //    fprintf (unit22, "%2d %8.3f %8.5f %9.3f %7.3f", n, rl0, rl0 * rscale,
01716 /*
01717           f0, par[4]*scale[4]);
01718 */
01719 /*        f0, sqrt (scale[4] * scale[4] * fabs (h[4 + npar*4])));
01720       for (i=2; i<4; i++) {
01721         fprintf (unit22, "%12.4e%12.4e", par[i] * scale[i],
01722         sqrt (scale[i] * scale[i] * fabs (h[i + npar*i])));
01723         if (calc_seconds) fprintf (unit22, "%12.4e",
01724             sqrt (scale[i] * scale[i] * fabs (hs[i + npar*i])));
01725       }*/
01726                            /*  set fit quality flag  */
01727       /*
01728       fprintf (unit22, "%5i", -1);
01729 
01730       fprintf (unit22, " %8.3f %4d%12.4e%12.4e", rki, bfgs_status, f, fm);
01731       for (i=0; i<2; i++) {
01732         fprintf (unit22, "%12.4e%12.4e", par[i] * scale[i],
01733         sqrt (scale[i] * scale[i] * fabs (h[i + npar*i])));
01734         if (calc_seconds) fprintf (unit22, "%12.4e",
01735             sqrt (scale[i] * scale[i] * fabs (hs[i + npar*i])));
01736       }
01737       dval = par[10];
01738       par[10] = pk1;
01739       h[10 + 10 *npar] *= 1.0e-8;
01740       for (i=4; i<npar; i++) {
01741         fprintf (unit22, "%12.4e%12.4e", par[i] * scale[i],
01742         sqrt (scale[i] * scale[i] * fabs (h[i + npar*i])));
01743         if (calc_seconds) fprintf (unit22, "%12.4e",
01744             sqrt (scale[i] * scale[i] * fabs (hs[i + npar*i])));
01745       }
01746       h[10 + 10 *npar] *= 1.0e8;
01747       par[10] = dval;*/
01748 /*
01749       fprintf (unit22, "%2d%10.3f%9.3f%9.3f%5d", n, f0, rki, rl0,
01750         bfgs_status);
01751       fprintf (unit22, "%12.4e%12.4e", f, fm);
01752 */
01753 /*
01754  *  scaling of asymmetry parameter, calibration of u_x, u_y parameters to m/s
01755  *    added 01.11.05
01756  */
01757 /*
01758       for (i=0; i<npar; i++)
01759         fprintf (unit22, "%12.4e", par[i] * scale[i]);
01760       for (i=0; i<npar; i++) {
01761         fprintf (unit22, "%12.4e",
01762         sqrt (scale[i] * scale[i] * fabs (h[i + npar*i])));
01763         if (calc_seconds) fprintf (unit22, "%12.4e",
01764           sqrt (scale[i] * scale[i] * fabs (hs[i + npar*i])));
01765       }
01766 */
01767 /*
01768       fprintf (unit22, "%12.4e%12.4e\n", pk, sn);
01769       if (extra_reporting)
01770         fprintf (unit22, "%5d%5d%5d", flnm_ct, bfgs_iter, linmin_iter);
01771 */
01772 //    fprintf (unit22, "\n");
01773       iers = 0;
01774       for (j = 0; j < npar; j++) par0[j] = par[j];
01775       f00 = f0;
01776     } else {
01777       if (fabs (par[2] - ux_guess) >= 500.0) reject_ux++;
01778       if (bfgs_status >= 100) reject_bfgs++;
01779       if (fabs (rl0 - rki) >= 50.0) reject_k++;
01780       if (fm >= 26.5) reject_fm++;
01781       if (sn <= -2.0 || sn >= 45.0) reject_sn++;
01782       if (rl0 <= 40.0) reject_l++;
01783       reject_tot++;
01784     }
01785     if ((fit_ct % 10) == 9) time_elapsed (fit_ct, &total_clock); 
01786       }
01787                 /* Added write statements (moved from inside n loop) 3/24/2011 */
01788       autoweed_vel(n_fit, l_fit, ux_fit, uy_fit, mask, cvg_ct);
01789       for(j=0; j<cvg_ct; j++)   {
01790     fprintf (unit22, "%2d %8.3f %8.5f %9.3f %7.3f", n, l_fit[j], l_fit[j] * rscale,
01791         freq_fit[j], sqrt(fabs(h_fit[j][4+npar*4]) * scale[4] * scale[4]));
01792     for (i=2; i<4; i++) {
01793       fprintf (unit22, "%12.4e%12.4e", par_fit[j][i] * scale[i],
01794           sqrt (scale[i] * scale[i] * fabs (h_fit[j][i + npar*i])));
01795       if (calc_seconds) fprintf (unit22, "%12.4e",
01796           sqrt (scale[i] * scale[i] * fabs (hs_fit[j][i + npar*i])));
01797     }
01798     if(mask[j]) fprintf (unit22, "%5i", 0);
01799     else fprintf (unit22, "%5i", 1);
01800     fprintf (unit22, " %8.3f %4d%12.4e%12.4e", rki_fit[j], bfgs_fit[j], f_fit[j], fm_fit[j]);
01801     for (i=0; i<2; i++) {
01802       fprintf (unit22, "%12.4e%12.4e", par_fit[j][i] * scale[i],
01803           sqrt (scale[i] * scale[i] * fabs (h_fit[j][i + npar*i])));
01804       if (calc_seconds) fprintf (unit22, "%12.4e",
01805           sqrt (scale[i] * scale[i] * fabs (hs_fit[j][i + npar*i])));
01806     }
01807     for (i=4; i<npar; i++) {
01808       fprintf (unit22, "%12.4e%12.4e", par_fit[j][i] * scale[i],
01809           sqrt (scale[i] * scale[i] * fabs (h_fit[j][i + npar*i])));
01810       if (calc_seconds) fprintf (unit22, "%12.4e",
01811           sqrt (scale[i] * scale[i] * fabs (hs_fit[j][i + npar*i])));
01812     }
01813     fprintf (unit22, "\n");
01814       }
01815 
01816 
01817       fprintf (runlog, "End iteration n = %d; time = %s\n", n,
01818           time_elapsed (0, &total_clock));
01819       fprintf (runlog, "  %d converged", cvg_ct);
01820       if (drms_output) {
01821     sprintf (key, "n_%d_fits", n);
01822     drms_setkey_int (orec, key, cvg_ct);
01823       }
01824       if (cvg_ct) {
01825     if (drms_output) {
01826       sprintf (key, "n_%d_fmin", n);
01827       drms_setkey_float (orec, key, cvg_minf * dnu);
01828       sprintf (key, "n_%d_fmax", n);
01829       drms_setkey_float (orec, key, cvg_maxf * dnu);
01830     }
01831     if (cvg_ct > 1)
01832       fprintf (runlog, ", first @ %.0f, last @ %.0f", cvg_minf * dnu, cvg_maxf * dnu);
01833     else fprintf (runlog, ", @ %.0f", cvg_minf * dnu);
01834       }
01835       fprintf (runlog, "\n");
01836       fprintf (runlog,
01837           "  rejected: %d (ux); %d (bfgs); %d (k); %d (fm); %d (sn); %d (l); %d (total)\n",
01838           reject_ux, reject_bfgs, reject_k, reject_fm, reject_sn, reject_l,
01839       reject_tot);
01840     }
01841     fclose (unit22);
01842 
01843     if (drms_output) {
01844       int kstat = 0;
01845       int keyct = sizeof (propagate) / sizeof (char *);
01846       char *key_str;
01847       double key_dbl;
01848       int key_n, key_int, crcl_known = 1;
01849       TIME key_tim;
01850                     /*  copy designated keywords from input  */
01851       for (key_n = 0; key_n < propct; key_n++) {
01852     if (strcmp (copykeylist[key_n], "+"))
01853       kstat += check_and_copy_key (orec, irec, copykeylist[key_n]);
01854     else kstat += propagate_keys (orec, irec, propagate, keyct);
01855       }
01856          /*  if necessary, construct CarrTime from CR:CL, or vice-versa  */
01857       key_str = drms_getkey_string (irec, "CarrTime", &status);
01858       if (status) {
01859         key_int = drms_getkey_int (irec, "CarrRot", &status);
01860     if (status) crcl_known = 0;
01861     else {
01862       key_dbl = drms_getkey_double (irec, "CMLon", &status);
01863       if (status) crcl_known = 0;
01864     }
01865     if (crcl_known) {
01866       char CarrTime[9];
01867       snprintf (CarrTime, 9, "%d:%03.0f", key_int, key_dbl);
01868       kstat += check_and_set_key_str (orec, "CarrTime", CarrTime);
01869     }
01870       } else {
01871     key_int = drms_getkey_int (irec, "CarrRot", &status);
01872     if (status) {
01873       sscanf (key_str, "%d:%lf", &key_int, &key_dbl);
01874       kstat += check_and_set_key_int (orec, "CarrRot", key_int);
01875       kstat += check_and_set_key_double (orec, "CMLon", key_dbl);
01876     }
01877       }
01878 
01879       kstat += check_and_set_key_str (orec, "Module", module_ident);
01880       kstat += check_and_set_key_str (orec, "Source", source);
01881       kstat += check_and_set_key_time (orec, "Created", CURRENT_SYSTEM_TIME);
01882       kstat += check_and_set_key_str (orec, "BLD_VERS", jsoc_version);
01883       kstat += check_and_set_key_float (orec, "ux_guess", ux_guess);
01884       kstat += check_and_set_key_float (orec, "uy_guess", uy_guess);
01885       kstat += check_and_set_key_float (orec, "asym_guess", asym_guess);
01886       if (kstat) {
01887     drms_sprint_rec_query (recid, orec);
01888     fprintf (stderr, "Error writing key values(s) to record %s\n", recid);
01889     fprintf (stderr, "      series may not have appropriate structure\n");
01890         drms_close_record (orec, DRMS_FREE_RECORD);
01891         continue;
01892       }
01893       gethostname (line, 1024);
01894       fprintf (runlog, "Total wall clock time = %d sec on %s\n", total_clock, line);
01895       if (runlog != stderr) fclose (runlog);
01896       drms_close_record (orec, dispose);
01897     }
01898   }
01899   drms_close_records (ids, DRMS_FREE_RECORD);
01900   return 0;
01901 }
01902 
01903 /*
01904  *  Revision History (all mods by Rick Bogart unless otherwise noted)
01905  *
01906  *  v 0.7  09.07.14     ported from SOI ringfit v 0.86
01907  *    0.8  09.08.25     fixed scaling of ux_guess and uy_guess and test
01908  *          on par[2]; trimmed output format slightly; fixed bug in
01909  *          setting of Dnu from Domega
01910  *    0.9  09.09.04     output to DRMS as possible, with processing of
01911  *          multiple records; copying and setting of a few keys;
01912  *          modified output format to conform to standard
01913  *  09.11.04        fixed key setting; modified to use and prefer
01914  *          WCS keywords for frequency scaling, and to no longer
01915  *          use propagated mapping keys from tracking
01916  *  09.11.16        corrected rscale
01917  *    1.0  09.12.09     added numerous keywords to propagation list
01918  *      10.01.07        removed constraints on values of A0, B1, and
01919  *              B2, so that fits are insensitive to power spectrum
01920  *              normalization
01921  *  10.02.01        fixed treatment of params_get_str()    
01922  *              fixed setting of Source key and diagnostics and
01923  *          error recovery for multiple region processing
01924  *  10.02.21        minor stderr format mods; added Created key
01925  *  10.03.05        added BLD_VERS key
01926  *  10.04.22        corrected rscale (again! maybe it's right this
01927  *          time); changed defaults for nmax, lmin, and lmax to be
01928  *          more appropriate for HMI data; added quality flag column
01929  *          to output; removed default values for input dataset and
01930  *          guessfile; added three more apodization keywords to
01931  *          propagation list
01932  *  v 1.0 frozen 2010.04.23
01933  *  10.06.12        Added setting of n_x_fits, n_x_fmin, n_x_fmax;
01934  *          keywords; added Ident to list of default propagated
01935  *          keywords; added option for copying values of different
01936  *          or additional keywords from input
01937  *  10.09.14        fixed bug in log base calculation affecting
01938  *          unnatural logs
01939  *              Added parameter for asymmetry term initial
01940  *          guess (including turning off fit); added setting of
01941  *          keywords for parameters
01942  *  10.11.01        Added parameter for asymmetry term (to turn
01943  *          off fit; just specifying initial guess does not do it)
01944  *          Also added parameters for turning off fit and initial
01945  *          guess for A1
01946  *  10.11.16        Fixed reporting value of parameter p; changed
01947  *          default fmax from 7000 to 5500
01948  *  10.11.30        Additional fix to take care of maintaining p
01949  *          parameter value
01950  *  v 1.1 frozen 2010.12.02
01951  *  11.03.23        Changed output column 5 from w0 to D_w0
01952  *  11.03.25        Include autoweeding of velocities to set values
01953  *          in qual column
01954  *  11.04.18        Fixed bug in determination of default input
01955  *          segment when multiple segments possible
01956  *  11.06.06        ditto
01957  *              Added optional rectangular filtering on spectrum
01958  *  v 1.2 frozen 2011.06.09
01959  */