00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include <jsoc_main.h>
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <math.h>
00046 #include <sys/time.h>
00047 #include <sys/resource.h>
00048
00049 #define DTOR (M_PI / 180.)
00050
00051 #define ARRLENGTH(ARR) (sizeof(ARR) / sizeof(ARR[0]))
00052 #define DIE(msg) {fflush(stdout); fprintf(stderr, "%s, status=%d\n", msg, status); return(status);}
00053 #define SHOW(msg) {printf("%s", msg); fflush(stdout);}
00054
00055 #include "pfss.h"
00056 #include "pfss_pkg.h"
00057 #include "timing.c"
00058
00059 char *module_name = "jgh2pfss";
00060
00061 ModuleArgs_t module_args[] =
00062 {
00063 {ARG_STRING, "in", NULL, "Input data series."},
00064 {ARG_STRING, "out", NULL, "Output data series."},
00065 {ARG_INT, "LMAX", "72", "Maximum l number desired."},
00066 {ARG_INT, "SINLAT", "0", "Set to zero for latitude grid; lat by default."},
00067 {ARG_INT, "NP", "360", "Number of Grid points in longitude."},
00068 {ARG_INT, "NT", "180", "Number of Grid points in latitude."},
00069 {ARG_INT, "NR", "16", "Number of Grid points in radial direction."},
00070 {ARG_INT, "VERB", "1", "Level of verbosity: 0=errors/warnings; 1=all messages"},
00071 {ARG_END}
00072 };
00073
00074
00075
00076
00077 int DoIt(void)
00078 {
00079 int status = DRMS_SUCCESS;
00080 char *inQuery, *outQuery;
00081 DRMS_RecordSet_t *inRS, *outRS;
00082 int irec, nrecs;
00083 DRMS_Record_t *inRec, *outRec, *linkedRec;
00084 DRMS_Segment_t *inSeg_g, *inSeg_h;
00085 DRMS_Segment_t *outSegBr, *outSegBt, *outSegBp;
00086 DRMS_Array_t *inArray_g, *inArray_h;
00087 DRMS_Array_t *outArrayBr, *outArrayBt, *outArrayBp;
00088 DRMS_Link_t *synopLink, *ghLink;
00089 float *inData_g, *inData_h;
00090 float *g, *h;
00091 float *Br, *Bt, *Bp;
00092
00093 char *inchecklist[] = {"CAR_ROT", "LON0", "LMAX"};
00094 char *outchecklist[] = {"CAR_ROT", "LON0", "LMAX", "CTYPE2"};
00095 DRMS_Record_t *inRec_tmp, *outRec_tmp;
00096 DRMS_Keyword_t *inkeytest, *outkeytest;
00097
00098 int verbflag;
00099 int outDims[3];
00100 int lmax, lmax1;
00101 int want_lmax, have_lmax;
00102 int car_rot, sinlat;
00103 int np, nt, nr;
00104
00105 int i, j, n;
00106 int l, m, itest;
00107 int n0, n1, n2;
00108
00109 struct Grid *grid;
00110
00111
00112 double wt0, wt1, wt;
00113 double ut0, ut1, ut;
00114 double st0, st1, st;
00115 double ct0, ct1, ct;
00116 wt0 = getwalltime();
00117 ct0 = getcputime(&ut0, &st0);
00118
00119
00120 inQuery = (char *)params_get_str(&cmdparams, "in");
00121 outQuery = (char *)params_get_str(&cmdparams, "out");
00122 want_lmax = params_get_int(&cmdparams, "LMAX");
00123 sinlat = params_get_int(&cmdparams, "SINLAT");
00124 np = params_get_int(&cmdparams, "NP");
00125 nt = params_get_int(&cmdparams, "NT");
00126 nr = params_get_int(&cmdparams, "NR");
00127 verbflag = params_get_int(&cmdparams, "VERB");
00128
00129 outDims[0] = np; outDims[1] = nt; outDims[2] = nr;
00130
00131
00132 inRS = drms_open_records(drms_env, inQuery, &status);
00133 if (status || inRS->n == 0) DIE("No input data found");
00134 nrecs = inRS->n;
00135
00136
00137 inRec_tmp = inRS->records[0];
00138 for (itest = 0; itest < ARRLENGTH(inchecklist); itest++) {
00139 inkeytest = hcon_lookup_lower(&inRec_tmp->keywords, inchecklist[itest]);
00140 if (!inkeytest) {
00141 fprintf(stderr, "ERROR: required input keyword %s is missing\n", inchecklist[itest]);
00142 return 0;
00143 }
00144 }
00145 have_lmax = drms_getkey_int(inRec_tmp, "LMAX", &status);
00146 if (want_lmax > have_lmax) DIE("Request lmax greater than input lmax");
00147
00148
00149 outRec_tmp = drms_create_record(drms_env, outQuery, DRMS_TRANSIENT, &status);
00150 if (status) DIE("Unable to open record in output series");
00151 for (itest = 0; itest < ARRLENGTH(outchecklist); itest++) {
00152 outkeytest = hcon_lookup_lower(&outRec_tmp->keywords, outchecklist[itest]);
00153 if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1) {
00154 fprintf(stderr, "Output keyword %s is missing/constant/a link.\n", outchecklist[itest]);
00155 return 0;
00156 }
00157 }
00158 drms_close_record(outRec_tmp, DRMS_FREE_RECORD);
00159
00160
00161 outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00162 if (status) DIE("Output recordset not created");
00163
00164
00165 grid = (struct Grid *)(malloc(sizeof(struct Grid)));
00166 grid->np = np; grid->nt = nt; grid->nr = nr;
00167 grid->ph = (float *)(calloc(np, sizeof(float)));
00168 grid->th = (float *)(calloc(nt, sizeof(float)));
00169 grid->rr = (float *)(calloc(nr, sizeof(float)));
00170 for (i = 0; i < np; i++)
00171 grid->ph[i] = (i + 0.5) / np * 2 * M_PI;
00172 for (j = 0; j < nt; j++)
00173 grid->th[j] = sinlat ? acos(1.0 - 2.0 * (nt - j - 0.5) / nt) : ((nt - j - 0.5) / nt * M_PI);
00174 if (nr != 1) {
00175 for (n = 0; n < nr; n++)
00176 grid->rr[n] = (RS - R0) / (nr - 1) * n + R0;
00177 } else {
00178 grid->rr[0] = RS;
00179 printf("Only one point specified in radial direction, compute source surface.\n");
00180 }
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 for (irec = 0; irec < nrecs; irec++) {
00193
00194 if (verbflag) {
00195 wt1 = getwalltime();
00196 ct1 = getcputime(&ut1, &st1);
00197 printf("processing record %d...\n", irec);
00198 }
00199
00200
00201 inRec = inRS->records[irec];
00202 inSeg_g = drms_segment_lookup(inRec, "g");
00203 inArray_g = drms_segment_read(inSeg_g, DRMS_TYPE_FLOAT, &status);
00204 if (status) DIE("No g data file found. \n");
00205 inSeg_h = drms_segment_lookup(inRec, "h");
00206 inArray_h = drms_segment_read(inSeg_h, DRMS_TYPE_FLOAT, &status);
00207 if (status) DIE("No h data file found. \n");
00208 inData_g = (float *)inArray_g->data; inData_h = (float *)inArray_h->data;
00209
00210
00211 car_rot = drms_getkey_int(inRec, "CAR_ROT", &status);
00212 n0 = inArray_g->axis[0];
00213 n1 = inArray_g->axis[1];
00214 if (inArray_g->axis[0] != inArray_g->axis[1] ||
00215 inArray_h->axis[0] != inArray_h->axis[1] ||
00216 (inArray_g->axis[1] - 1) != have_lmax ||
00217 (inArray_h->axis[1] - 1) != have_lmax)
00218 DIE("Wrong input data dimension");
00219 if (want_lmax > have_lmax)
00220 DIE("Request lmax greater than available lmax");
00221 lmax = want_lmax;
00222 lmax1 = lmax + 1;
00223 g = (float *)(calloc(lmax1 * lmax1, sizeof(float)));
00224 h = (float *)(calloc(lmax1 * lmax1, sizeof(float)));
00225 for (int l = 0; l <= lmax; l++) {
00226 for (int m = 0; m <= l; m++) {
00227 g[l * lmax1 + m] = inData_g[l * n0 + m];
00228 h[l * lmax1 + m] = inData_h[l * n0 + m];
00229 }
00230 }
00231 g[0] = 0.0; h[0] = 0.0;
00232 drms_free_array(inArray_g); drms_free_array(inArray_h);
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 outArrayBr = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00251 outArrayBt = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00252 outArrayBp = drms_array_create(DRMS_TYPE_FLOAT, 3, outDims, NULL, &status);
00253 Br = (float *)outArrayBr->data;
00254 Bt = (float *)outArrayBt->data;
00255 Bp = (float *)outArrayBp->data;
00256
00257
00258 if (verbflag) SHOW("start pfss computation... ");
00259 pfss(g, h, grid, Br, Bt, Bp, lmax, 0.0);
00260 if (verbflag) SHOW("done.\n");
00261 free(g); free(h);
00262
00263
00264 outRec = outRS->records[irec];
00265 outSegBr = drms_segment_lookupnum(outRec, 0);
00266 outSegBt = drms_segment_lookupnum(outRec, 1);
00267 outSegBp = drms_segment_lookupnum(outRec, 2);
00268 for (i = 0; i < 3; i++) {
00269 outSegBr->axis[i] = outArrayBr->axis[i];
00270 outSegBt->axis[i] = outArrayBt->axis[i];
00271 outSegBp->axis[i] = outArrayBp->axis[i];
00272 }
00273 outArrayBr->parent_segment = outSegBr;
00274 outArrayBt->parent_segment = outSegBt;
00275 outArrayBp->parent_segment = outSegBp;
00276
00277
00278 status = drms_segment_write(outSegBr, outArrayBr, 0);
00279 if (status) DIE("problem writing Br");
00280 drms_free_array(outArrayBr);
00281 status = drms_segment_write(outSegBt, outArrayBt, 0);
00282 if (status) DIE("problem writing Bt");
00283 drms_free_array(outArrayBt);
00284 status = drms_segment_write(outSegBp, outArrayBp, 0);
00285 if (status) DIE("problem writing Bp");
00286 drms_free_array(outArrayBp);
00287
00288
00289 drms_setkey_string(outRec, "BLD_VERS", jsoc_version);
00290 drms_setkey_time(outRec, "DATE", CURRENT_SYSTEM_TIME);
00291 drms_copykey(outRec, inRec, "CAR_ROT");
00292 drms_copykey(outRec, inRec, "LON0");
00293 drms_setkey_int(outRec, "LMAX", lmax);
00294 drms_setkey_float(outRec, "RSS", RS);
00295 drms_setkey_float(outRec, "CDELT1", - 360. / np);
00296 drms_setkey_float(outRec, "CRPIX1", 1.0);
00297 drms_setkey_float(outRec, "CRVAL1", car_rot * 360. - 180. / np);
00298 drms_setkey_float(outRec, "CRPIX2", nt / 2.0 + 0.5);
00299 drms_setkey_float(outRec, "CRVAL2", 0.0);
00300 if (sinlat) {
00301 drms_setkey_string(outRec, "CTYPE2", "CRLT-CEA");
00302 drms_setkey_float(outRec, "CDELT2", 2. / nt);
00303 drms_setkey_string(outRec, "CUNIT2", "Sine latitude");
00304 } else {
00305 drms_setkey_string(outRec, "CTYPE2", "CRLT-CAR");
00306 drms_setkey_float(outRec, "CDELT2", 180. / nt);
00307 drms_setkey_string(outRec, "CUNIT2", "degree");
00308 }
00309 drms_setkey_float(outRec, "CDELT3", (RS - R0) / (nr - 1));
00310 drms_setkey_float(outRec, "CRPIX3", 1.0);
00311 drms_setkey_float(outRec, "CRVAL3", 1.0);
00312
00313
00314 ghLink = hcon_lookup_lower(&outRec->links, "HARMONICS");
00315 if (ghLink) {
00316 drms_setlink_static(outRec, "HARMONICS", inRec->recnum);
00317 linkedRec = drms_link_follow(outRec, "HARMONICS", &status);
00318 }
00319 if (!status) {
00320 synopLink = hcon_lookup_lower(&linkedRec->links, "SYNOP");
00321 if (synopLink) {
00322 drms_setlink_static(outRec, "SYNOP", linkedRec->recnum);
00323 }
00324 }
00325
00326
00327 if (verbflag) {
00328 wt = getwalltime();
00329 ct = getcputime(&ut, &st);
00330 printf("record %d done, %.2f ms wall time, %.2f ms cpu time\n",
00331 irec, wt - wt1, ct - ct1);
00332 }
00333 }
00334
00335 free(grid->rr); free(grid->th); free(grid->ph); free(grid);
00336
00337 drms_close_records(inRS, DRMS_FREE_RECORD);
00338 drms_close_records(outRS, DRMS_INSERT_RECORD);
00339
00340 if (verbflag) {
00341 wt = getwalltime();
00342 ct = getcputime(&ut, &st);
00343 printf("total time spent: %.2f ms wall time, %.2f ms cpu time\n",
00344 wt - wt0, ct - ct0);
00345 }
00346
00347 return(DRMS_SUCCESS);
00348
00349 }
