JSOC_Documentation: proj/mag/d4vm/apps/dave4vm4velocity.c Source File

00001 /*
00002  * Module name:     test_d4vm.c
00003  *
00004  * Description:
00005  *          DAVE4VM C version
00006  *
00007  * Calling:
00008  *                  
00009  *
00010  * Original by:     Peter Schuck, Jacob Hageman
00011  * Adapted by:      Xudong Sun (xudongs@stanford.edu)
00012  *
00013  * Version:
00014  *          v1.0        Mar 01 2010
00015  *          v1.1        Mar 15 2010
00016  *
00017 
00018  * Issues:
00019  *          v1.0
00020  *          1. Designed to deal with a continuous time series, computing result
00021  *          with each pair of neighbouring records.
00022  *          2. Added "cadence" keyword to stride over a minimum fixed period between
00023  *          two records, so the result can be computed at a lower cadence.
00024  *          3. Commented out #infdef USE_INTEL_LIB and DEBUG, seems not properly set up
00025  *          4. Will have to make sure the stack size is big enough (at least ~70MB)
00026  *          run "limit s u" to set unlimited
00027  *          v1.1
00028  *          Updated according to Jacob's new code
00029  *
00030  *
00031  * Example:
00032  *          test_d4vm_old 'in=su_xudong.d4vmin[2004.08.31_23:55:00_TAI/10m]'
00033  */
00034 
00035 
00036 #include <jsoc_main.h>
00037 #include <stdio.h>
00038 #include <stdlib.h>
00039 #include <math.h>
00040 #include <sys/time.h>
00041 #include <sys/resource.h>
00042 #include <fstats.h>
00043 
00044 #define PI  (M_PI)
00045 #define DTOR    (PI / 180.)
00046 
00047 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00048 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00049 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00050 
00051 #include "timing.c"
00052 
00053 
00054 
00055 /* ################## Prototypes ################## */
00056 
00057 extern void d4vm_preproc_(float *bx0, float *bx1, 
00058                           float *by0, float *by1, 
00059                           float *bz0, float *bz1, 
00060                           double *dx, double *dy, double *dt, 
00061                           float *bx, float *bxx, float *bxy, 
00062                           float *by, float *byx, float *byy, 
00063                           float *bz, float *bzx, float *bzy, float *bzt,
00064                           long *xsize, long *ysize);
00065 
00066 void d4vm_kernel_(double *dx, double *dy, 
00067                   long *ksize, float *k_th,
00068                   float *k_x, float *k_y, 
00069                   float *k_xx, float *k_yy, float *k_xy);
00070 
00071 void d4vm_matrix_(float *a, 
00072                   float *bx, float *bxx, float *bxy, 
00073                   float *by, float *byx, float *byy, 
00074                   float *bz, float *bzx, float *bzy, float *bzt, 
00075                   float *k_th, float *k_x, float *k_y,
00076                   float *k_xx, float *k_yy, float *k_xy,
00077                   long *xsize, long *ysize, long *ksize);
00078 
00079 void d4vm_solver_(float *a, 
00080                   float *u0, float *v0, float *w0, 
00081                   float *ux, float *vx, float *wx, 
00082                   float *uy, float *vy, float *wy,
00083                   long *xsize, long *ysize, double *threshold);
00084 
00085 
00086 char *module_name = "dave4vm4velocity"; /* Module name */
00087 
00088 ModuleArgs_t module_args[] =
00089 {
00090     {ARG_STRING, "in", "hmi.sharp_cea_720s", "Input data series."},
00091     {ARG_STRING, "out", "su_yang.hmi_dave4vm4velo",  "Output data series."},
00092     {ARG_INT, "ksz0", "19", "Width of window for convolution, must be odd."},
00093     {ARG_INT, "ksz1", "19", "Height of window for convolution, must be odd."},
00094     {ARG_DOUBLE, "threshold", "0.", "Threshold for resolving aperture problem."},
00095     {ARG_DOUBLE, "cadence", "240.0", "Minimum stride between two records."},
00096     {ARG_INT, "VERB", "2", "Level of verbosity: 0=errors/warnings; 1=status; 2=all"},
00097     {ARG_END}
00098 };
00099 
00100 
00101 /* ################## Main Module ################## */
00102 
00103 int DoIt(void)
00104 {
00105     int status = DRMS_SUCCESS;
00106     char *inQuery, *outQuery;
00107     DRMS_RecordSet_t *inRS;
00108     int irec, nrecs;
00109     DRMS_Record_t *inRec0, *inRec1, *outRec;
00110     DRMS_Segment_t *inSeg0, *inSeg1;
00111     DRMS_Segment_t *outSegV; 
00112     DRMS_Array_t *inArrayX0, *inArrayY0, *inArrayZ0, *inArrayX1, *inArrayY1, *inArrayZ1;
00113     DRMS_Array_t *outArrayV, *outArrayVy; 
00114     float *outDataVx, *outDataVy, *outDataVz; 
00115     int harpnum;
00116     
00117     float *Bx0, *Bx1, *By0, *By1, *Bz0, *Bz1;               // ancillary
00118     float *Bx, *Bxx, *Bxy, *By, *Byx, *Byy, *Bz, *Bzx, *Bzy, *Bzt;  // processed
00119     float *k_th, *k_x, *k_y, *k_xx, *k_yy, *k_xy;           // convolution kernels
00120     float *a;                               // matrix
00121     float *Vx, *Vy, *Vz;                        // outputs
00122     float *DVxDx, *DVyDx, *DVzDx;
00123     float *DVxDy, *DVyDy, *DVzDy;
00124     float onedegree4L = 6.955e5 * DTOR; // 1 degree corresponse to km
00125 
00126     char *inchecklist[] = {"T_REC","T_OBS","HARPNUM"};      // For now
00127     char *outchecklist[] = {"T_REC", "K0", "K1", "THRESH","HARPNUM"};
00128     DRMS_Record_t *inRec_tmp, *outRec_tmp;
00129     DRMS_Keyword_t *inkeytest, *outkeytest;
00130 
00131     int outDims[2] = {0, 0};
00132     long nx, ny, nxny, ksize[2];    // Seems compatible with long in Fortran?
00133     int verbflag;
00134     double threshold, cadence;
00135     TIME t0, t1;    // Need to check format
00136     double dx, dy, dt;
00137 
00138     int i, l, m, itest;
00139     int nsuccess = 0;
00140     
00141     // Time measuring
00142     double wt0, wt1, wt;
00143     double ut0, ut1, ut;
00144     double st0, st1, st;
00145     double ct0, ct1, ct;
00146     wt0 = getwalltime();
00147     ct0 = getcputime(&ut0, &st0);
00148 
00149     /* Get parameters */
00150     inQuery = (char *)params_get_str(&cmdparams, "in");
00151     outQuery = (char *)params_get_str(&cmdparams, "out");
00152     ksize[0] = params_get_int(&cmdparams, "ksz0");
00153     ksize[1] = params_get_int(&cmdparams, "ksz1");
00154     threshold = params_get_double(&cmdparams, "threshold");
00155     cadence = params_get_double(&cmdparams, "cadence");
00156     verbflag = params_get_int(&cmdparams, "VERB");
00157     for (i = 0; i < 2; i++) {
00158         if (ksize[i] % 2 != 1) {
00159             DIE("Kernel size must be odd");
00160         }
00161     }
00162 
00163     /* Open input */
00164     inRS = drms_open_records(drms_env, inQuery, &status);
00165     if (status || inRS->n == 0) DIE("No input data found");
00166     if (inRS->n < 2) DIE("At least 2 records are needed");
00167     nrecs = inRS->n;
00168 
00169 //printf("nrecs=%d\n", nrecs);
00170 
00171     /* Check for essential input keywords */
00172     inRec_tmp = inRS->records[0];
00173     for (itest = 0; itest < ARRLENGTH(inchecklist); itest++) {
00174         inkeytest = hcon_lookup_lower(&inRec_tmp->keywords, inchecklist[itest]);
00175         if (!inkeytest) {
00176             fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
00177             return 0;
00178         }
00179     }
00180 
00181     /* Check for essential output keywords */
00182     outRec_tmp = drms_create_record(drms_env, outQuery, DRMS_TRANSIENT, &status);
00183     if (status) DIE("Unable to open record in output series");
00184     for (itest = 0; itest < ARRLENGTH(outchecklist); itest++) {
00185         outkeytest = hcon_lookup_lower(&outRec_tmp->keywords, outchecklist[itest]);
00186         if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1) {
00187             fprintf(stderr, "Output keyword %s is missing/constant/a link.\n", outchecklist[itest]);
00188             return 0;
00189         }
00190     }
00191     drms_close_record(outRec_tmp, DRMS_FREE_RECORD);
00192 
00193     /* Read in the first record and save it */
00194     inRec0 = inRS->records[0];
00195     t0 = drms_getkey_time(inRec0, "T_REC", NULL);   // NEED TO CHECK!!!!!
00196     harpnum = drms_getkey_int(inRec0, "HARPNUM", NULL);
00197     
00198     inSeg0 = drms_segment_lookup(inRec0, "Bp");
00199     inArrayX0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00200     Bx0 = (float *)inArrayX0->data;
00201     nx = inArrayX0->axis[0]; ny = inArrayX0->axis[1]; nxny = nx * ny;
00202     outDims[0] = nx; outDims[1] = ny;
00203     inSeg0 = drms_segment_lookup(inRec0, "Bt");
00204     inArrayY0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00205     By0 = (float *)inArrayY0->data;
00206     for (int ii = 0; ii < nxny; ii++) By0[ii] *= -1.;
00207     inSeg0 = drms_segment_lookup(inRec0, "Br");
00208     inArrayZ0 = drms_segment_read(inSeg0, DRMS_TYPE_FLOAT, &status);
00209     Bz0 = (float *)inArrayZ0->data;
00210 
00211     /* Do this for record #1 to #nrecs-1 */
00212     for (irec = 1; irec < nrecs; irec++)
00213     {
00214         /* Input record and data */
00215         inRec1 = inRS->records[irec];
00216         t1 = drms_getkey_time(inRec1, "T_REC", NULL);   // NEED TO CHECK!!!!!
00217         dt = t1 - t0;
00218         if (dt < cadence) {
00219             printf("Cadence higher than limit, skipping current record. \n");
00220             drms_free_array(inArrayX1);
00221             continue;
00222         }
00223         if (harpnum != drms_getkey_time(inRec1, "HARPNUM", NULL)) {
00224             printf("HARP number incorrect, skipping current record. \n");
00225             drms_free_array(inArrayX1);
00226             continue;
00227         }
00228 
00229         inSeg1 = drms_segment_lookup(inRec1, "Bp"); // Single segment
00230         inArrayX1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00231 //printf("t0=%f, t1=%f\n", t0, t1);
00232 //printf("dt=%f\n",t1-t0);
00233         Bx1 = (float *)inArrayX1->data;
00234         if (inArrayX1->axis[0] != nx || inArrayX1->axis[1] != ny)
00235             DIE("Wrong data dimension for current record, stop. \n"); 
00236         inSeg1 = drms_segment_lookup(inRec1, "Bt");
00237         inArrayY1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00238         By1 = (float *)inArrayY1->data;
00239         for (int ii = 0; ii < nxny; ii++) By1[ii] *= -1.;
00240         inSeg1 = drms_segment_lookup(inRec1, "Br");
00241         inArrayZ1 = drms_segment_read(inSeg1, DRMS_TYPE_FLOAT, &status);
00242         Bz1 = (float *)inArrayZ1->data;
00243         
00244         /* Time measure */
00245         if (verbflag) {
00246             wt1 = getwalltime();
00247             ct1 = getcputime(&ut1, &st1);
00248             printf("processing record %d...\n", nsuccess);
00249         }
00250 
00251         /* ======================== */
00252         /* This is the working part.*/
00253         /* ======================== */
00254 
00255         // Preproc
00256         Bx = (float *)calloc(sizeof(float), nxny);
00257         Bxx = (float *)calloc(sizeof(float), nxny);
00258         Bxy = (float *)calloc(sizeof(float), nxny);
00259         By = (float *)calloc(sizeof(float), nxny);
00260         Byx = (float *)calloc(sizeof(float), nxny);
00261         Byy = (float *)calloc(sizeof(float), nxny);
00262         Bz = (float *)calloc(sizeof(float), nxny);
00263         Bzx = (float *)calloc(sizeof(float), nxny);
00264         Bzy = (float *)calloc(sizeof(float), nxny);
00265         Bzt = (float *)calloc(sizeof(float), nxny);
00266 
00267         DVxDx = (float *)calloc(sizeof(float), nxny);
00268         DVyDx = (float *)calloc(sizeof(float), nxny);
00269         DVzDx = (float *)calloc(sizeof(float), nxny);
00270         DVxDy = (float *)calloc(sizeof(float), nxny);
00271         DVyDy = (float *)calloc(sizeof(float), nxny);
00272         DVzDy = (float *)calloc(sizeof(float), nxny);
00273 
00274         // Kernel
00275         k_th = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00276         k_x = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00277         k_y = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00278         k_xx = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00279         k_yy = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00280         k_xy = (float *)calloc(sizeof(float), ksize[0] * ksize[1]);
00281 
00282         // Matrix
00283         a = (float *)calloc(sizeof(float), 10 * 10 * nxny);
00284 
00285         // Output
00286         Vx = (float *)calloc(sizeof(float), nxny);
00287         Vy = (float *)calloc(sizeof(float), nxny);
00288         Vz = (float *)calloc(sizeof(float), nxny);
00289 
00290         dx = onedegree4L * drms_getkey_float(inRec0, "CDELT1", &status);
00291         dy = onedegree4L * drms_getkey_float(inRec0, "CDELT2", &status);
00292 
00293 //printf("dx = %f, dy = %f\n", dx, dy);
00294 //SHOW("111\n");
00295         // Preprocess data
00296         d4vm_preproc_(Bx0, Bx1, By0, By1, Bz0, Bz1, 
00297                       &dx, &dy, &dt, 
00298                       Bx, Bxx, Bxy, By, Byx, Byy, Bz, Bzx, Bzy, Bzt,
00299                       &nx, &ny);
00300 //SHOW("112\n");
00301         // Generate kernel
00302         d4vm_kernel_(&dx, &dy, ksize,
00303                      k_th, k_x, k_y, k_xx, k_yy, k_xy);
00304 //SHOW("113\n");
00305 
00306         // Process data
00307         d4vm_matrix_(a, 
00308                      Bx, Bxx, Bxy, By, Byx, Byy, Bz, Bzx, Bzy, Bzt,
00309                      k_th, k_x, k_y, k_xx, k_yy, k_xy,
00310                      &nx, &ny, ksize);
00311 //SHOW("114\n");
00312         // Solve matrix
00313         d4vm_solver_(a, 
00314                      Vx, Vy, Vz, 
00315                      DVxDx, DVyDx, DVzDx, 
00316                      DVxDy, DVyDy, DVzDy, 
00317                      &nx, &ny, &threshold);
00318 //SHOW("115\n");
00319 
00320         /* Output record, create one at a time */
00321         outRec = drms_create_record(drms_env, outQuery, DRMS_PERMANENT, &status);
00322         if (status) DIE("Output record not created");
00323 
00324         outSegV = drms_segment_lookup(outRec, "Vp");
00325         outArrayV = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, Vx, &status);
00326         outArrayV->israw = 0;
00327         outArrayV->bzero = 0.0;
00328         outArrayV->bscale = 0.0001;
00329         /* Result writing */
00330         status = drms_segment_write(outSegV, outArrayV, 0);
00331         if (status) DIE("Problem writing file V");
00332         free(outArrayV);
00333 
00334         outSegV = drms_segment_lookup(outRec, "Vt");
00335         int ii;
00336         for (ii = 0; ii < nxny; ii++) Vy[ii] *= -1.0;
00337         outArrayV = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, Vy, &status);
00338         outArrayV->israw = 0;
00339         outArrayV->bzero = 0.0;
00340         outArrayV->bscale = 0.0001;
00341         /* Result writing */
00342         status = drms_segment_write(outSegV, outArrayV, 0);
00343         if (status) DIE("Problem writing file V");
00344         free(outArrayV);
00345 
00346         outSegV = drms_segment_lookup(outRec, "Vr");
00347         outArrayV = drms_array_create(DRMS_TYPE_FLOAT, 2, outDims, Vz, &status);
00348         outArrayV->israw = 0;
00349         outArrayV->bzero = 0.0;
00350         outArrayV->bscale = 0.0001;
00351         /* Result writing */
00352         status = drms_segment_write(outSegV, outArrayV, 0);
00353         if (status) DIE("Problem writing file V");
00354         free(outArrayV);
00355 
00356         /* Set keywords */
00357         drms_copykey(outRec, inRec0, "HARPNUM");
00358         drms_copykey(outRec, inRec0, "T_REC");
00359         drms_copykey(outRec, inRec0, "T_OBS");
00360 // add aditional keywords
00361         drms_copykey(outRec, inRec0, "DATE__OBS");
00362         drms_copykey(outRec, inRec0, "CRPIX1");
00363         drms_copykey(outRec, inRec0, "CRPIX2");
00364         drms_copykey(outRec, inRec0, "CRVAL1");
00365         drms_copykey(outRec, inRec0, "CRVAL2");
00366         drms_copykey(outRec, inRec0, "CDELT1");
00367         drms_copykey(outRec, inRec0, "CDELT2");
00368         drms_copykey(outRec, inRec0, "IMCRPIX1");
00369         drms_copykey(outRec, inRec0, "IMCRPIX2");
00370         drms_copykey(outRec, inRec0, "IMCRVAL1");
00371         drms_copykey(outRec, inRec0, "IMCRVAL2");
00372         drms_copykey(outRec, inRec0, "CROTA2");
00373         drms_copykey(outRec, inRec0, "CRDER1");
00374         drms_copykey(outRec, inRec0, "CRDER2");
00375         drms_copykey(outRec, inRec0, "CSYSER1");
00376         drms_copykey(outRec, inRec0, "CSYSER2");
00377         drms_copykey(outRec, inRec0, "DSUN_OBS");
00378         drms_copykey(outRec, inRec0, "CRLN_OBS");
00379         drms_copykey(outRec, inRec0, "CRLT_OBS");
00380         drms_copykey(outRec, inRec0, "CAR_ROT");
00381         drms_copykey(outRec, inRec0, "OBS_VR");
00382         drms_copykey(outRec, inRec0, "OBS_VW");
00383         drms_copykey(outRec, inRec0, "OBS_VN");
00384         drms_copykey(outRec, inRec0, "RSUN_OBS");
00385         drms_copykey(outRec, inRec0, "INSTRUME");
00386         drms_copykey(outRec, inRec0, "CAMERA");
00387         drms_copykey(outRec, inRec0, "LATDTMIN");
00388         drms_copykey(outRec, inRec0, "LONDTMIN");
00389         drms_copykey(outRec, inRec0, "LATDTMAX");
00390         drms_copykey(outRec, inRec0, "LONDTMAX");
00391         drms_copykey(outRec, inRec0, "OMEGA_DT");
00392         drms_copykey(outRec, inRec0, "LAT_MIN");
00393         drms_copykey(outRec, inRec0, "LON_MIN");
00394         drms_copykey(outRec, inRec0, "LAT_MAX");
00395         drms_copykey(outRec, inRec0, "LON_MAX");
00396         drms_copykey(outRec, inRec0, "NOAA_AR");
00397         drms_copykey(outRec, inRec0, "NOAA_NUM");
00398         drms_copykey(outRec, inRec0, "NOAA_ARS");
00399 
00400         drms_setkey_time(outRec, "T_REC_STAR", drms_getkey_time(inRec0, "T_REC", NULL));
00401         drms_setkey_time(outRec, "T_OBS_STAR", drms_getkey_time(inRec0, "T_OBS", NULL));
00402         drms_setkey_time(outRec, "T_REC_STOP", drms_getkey_time(inRec1, "T_REC", NULL));
00403         drms_setkey_time(outRec, "T_OBS_STOP", drms_getkey_time(inRec1, "T_OBS", NULL));
00404 
00405 // add statistical keywords
00406 
00407         double statMin, statMax, statMedn, statMean, statSig, statSkew, statKurt;
00408         int statNgood, ipixels;
00409          if (fstats(outDims[0] * outDims[1], Vx, &statMin, &statMax, &statMedn, &statMean, &statSig,
00410             &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00411 
00412         ipixels = outDims[0] * outDims[1];
00413         drms_setkey_int(outRec, "TOTVALS_1", ipixels);
00414         drms_setkey_int(outRec, "DATAVALS_1", statNgood);
00415         ipixels = outDims[0] * outDims[1]-statNgood;
00416         drms_setkey_int(outRec, "MISSVALS_1", ipixels);
00417         drms_setkey_double(outRec, "DATAMIN_1", statMin);
00418         drms_setkey_double(outRec, "DATAMAX_1", statMax);
00419         drms_setkey_double(outRec, "DATAMEDN_1", statMedn);
00420         drms_setkey_double(outRec, "DATAMEAN_1", statMean);
00421         drms_setkey_double(outRec, "DATARMS_1", statSig);
00422         drms_setkey_double(outRec, "DATASKEW_1", statSkew);
00423         drms_setkey_double(outRec, "DATAKURT_1", statKurt);
00424 
00425          if (fstats(outDims[0] * outDims[1], Vy, &statMin, &statMax, &statMedn, &statMean, &statSig,
00426             &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00427 
00428         ipixels = outDims[0] * outDims[1];
00429         drms_setkey_int(outRec, "TOTVALS_2", ipixels);
00430         drms_setkey_int(outRec, "DATAVALS_2", statNgood);
00431         ipixels = outDims[0] * outDims[1]-statNgood;
00432         drms_setkey_int(outRec, "MISSVALS_2", ipixels);
00433         drms_setkey_double(outRec, "DATAMIN_2", statMin);
00434         drms_setkey_double(outRec, "DATAMAX_2", statMax);
00435         drms_setkey_double(outRec, "DATAMEDN_2", statMedn);
00436         drms_setkey_double(outRec, "DATAMEAN_2", statMean);
00437         drms_setkey_double(outRec, "DATARMS_2", statSig);
00438         drms_setkey_double(outRec, "DATASKEW_2", statSkew);
00439         drms_setkey_double(outRec, "DATAKURT_2", statKurt);
00440 
00441          if (fstats(outDims[0] * outDims[1], Vz, &statMin, &statMax, &statMedn, &statMean, &statSig,
00442             &statSkew, &statKurt, &statNgood)) printf("\n Statistics computation failed\n");
00443 
00444         ipixels = outDims[0] * outDims[1];
00445         drms_setkey_int(outRec, "TOTVALS_3", ipixels);
00446         drms_setkey_int(outRec, "DATAVALS_3", statNgood);
00447         ipixels = outDims[0] * outDims[1]-statNgood;
00448         drms_setkey_int(outRec, "MISSVALS_3", ipixels);
00449         drms_setkey_double(outRec, "DATAMIN_3", statMin);
00450         drms_setkey_double(outRec, "DATAMAX_3", statMax);
00451         drms_setkey_double(outRec, "DATAMEDN_3", statMedn);
00452         drms_setkey_double(outRec, "DATAMEAN_3", statMean);
00453         drms_setkey_double(outRec, "DATARMS_3", statSig);
00454         drms_setkey_double(outRec, "DATASKEW_3", statSkew);
00455         drms_setkey_double(outRec, "DATAKURT_3", statKurt);
00456 
00457 // end of statistical keywords
00458 
00459 
00460         drms_setkey_int(outRec, "K0", ksize[0]);
00461         drms_setkey_int(outRec, "K1", ksize[1]);
00462         drms_setkey_double(outRec, "THRESH", threshold);
00463         drms_close_record(outRec, DRMS_INSERT_RECORD);
00464 
00465         /* Clean up, move to next record */
00466         drms_free_array(inArrayX0);
00467         drms_free_array(inArrayY0);
00468         drms_free_array(inArrayZ0);
00469         free(Bx); free(Bxx); free(Bxy);
00470         free(By); free(Byx); free(Byy);
00471         free(Bz); free(Bzx); free(Bzy); free(Bzt);
00472         free(DVxDx); free(DVyDx); free(DVzDx);
00473         free(DVxDy); free(DVyDy); free(DVzDy);
00474         free(a);
00475         free(k_th); free(k_x); free(k_y);
00476         free(k_xx); free(k_yy); free(k_xy);
00477         free(Vy); free(Vz); free(Vx);
00478 
00479 // set the second data to the first for next calculation
00480 
00481         inRec0 = inRec1;
00482         inSeg0 = inSeg1;
00483         inArrayX0 = inArrayX1;
00484         inArrayY0 = inArrayY1;
00485         inArrayZ0 = inArrayZ1;
00486         Bx0 = Bx1;
00487         By0 = By1;
00488         Bz0 = Bz1;
00489         t0 = t1;
00490         
00491         /* Time measure */
00492         if (verbflag) {
00493             wt = getwalltime();
00494             ct = getcputime(&ut, &st);
00495             printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n", 
00496                      nsuccess, wt - wt1, ct - ct1);
00497         }
00498 
00499         nsuccess++;
00500     }
00501 
00502     drms_free_array(inArrayX1);
00503     drms_free_array(inArrayY1);
00504     drms_free_array(inArrayZ1);
00505     drms_close_records(inRS, DRMS_FREE_RECORD);
00506     return(DRMS_SUCCESS);
00507 }