1 tplarson 1.1 #include <stdio.h>
2 #include <math.h>
3 #include <stdlib.h>
4 #include <time.h>
5 #include <sys/time.h>
6 #include <sys/resource.h>
7 #include <fftw3.h>
8 #include "jsoc_main.h"
9 #include "fresize.h"
10
11 #define kNOTSPECIFIED "not specified"
12 #define ARRLENGTH(ARR) (sizeof(ARR)/sizeof(ARR[0]))
13 #define PI (M_PI)
14 #define RADSINDEG (PI/180)
|
15 tplarson 1.5 #define RTRUE (6.96000000e8) // meters
16 #define AU (149597870691) // meters/au
|
17 tplarson 1.2 #define TAU_A (499.004783806) // light travel time in seconds, = 1 AU/c
18 #define C (299792458) // speed of light, m/s
|
19 tplarson 1.1
20 char *module_name = "mkylms";
|
21 tplarson 1.5 char *cvsinfo_mkylms = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/mkylms.c,v 1.4 2013/04/28 09:14:04 tplarson Exp $";
|
22 tplarson 1.1
23 ModuleArgs_t module_args[] =
24 {
25 {ARG_STRING, "out", NULL, "output data series"},
26 {ARG_STRING, "histlink", "HISTORY", "name of link to ancillary dataseries for processing metadata"},
27 {ARG_INT, "PERM", "1", "set to 0 for transient records, nonzero for permanent records"},
28 {ARG_INT, "VERB", "1", "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL},
29 {ARG_STRING, "MODELIST", kNOTSPECIFIED, "file specifying modes to generate"},
30 {ARG_INT, "XPIXELS", "1024", "number of pixels in x direction"},
31 {ARG_INT, "YPIXELS", "1024", "number of pixels in y direction"},
32 {ARG_FLOAT, "SOLARP", "0.0", "p-angle in degrees"},
33 {ARG_FLOAT, "OBSB0", "0.0", "b-angle in degrees"},
34 {ARG_FLOAT, "OBSDIST", "1.0", "distance from the sun in au"},
35 {ARG_FLOAT, "XOFFSET", "0.0", "offset in pixels from center in x direction"},
36 {ARG_FLOAT, "YOFFSET", "0.0", "offset in pixels from center in y direction"},
37 {ARG_INT, "IINTEN", "0", "flag for making intensity images, make velocity images otherwise"},
|
38 tplarson 1.2 {ARG_FLOAT, "PHASE", "0.0", "phase in radians"},
39 {ARG_FLOAT, "FREQUENCY", "0.0", "frequency in Hz"},
40 {ARG_INT, "ILIMBDARK", "0", "flag to correct for limb darkening in intensity images"},
41 {ARG_DOUBLES,"coefs", "0.0", "limb darkening coeficients, 6 needed"},
42 {ARG_INT, "IHEIGHT", "0", "flag to correct for height of formation"},
|
43 tplarson 1.1
44 {ARG_FLOAT, "AIRY_CDOWN", "0.0", "ratio of the pixel nyquist frequency to the cutoff frequency of the Airy function"},
45 {ARG_INT, "AIRY_IAP", "0", "option for type of apodization"},
46 {ARG_INT, "AIRY_NAP", "1", "distance between sampled points in pixels"},
47 {ARG_INT, "AIRY_WID", "-1", "half width of kernel"},
48 {ARG_INT, "AIRY_NSUB", "1", "distance between sampled points in pixels"},
49 {ARG_INT, "AIRY_XOFF", "0", "offset in pixels to add to x index of input image"},
50 {ARG_INT, "AIRY_YOFF", "0", "offset in pixels to add to y index of input image"},
51
52 {ARG_INT, "NBIN", "-1", "factor for binning"},
53 {ARG_INT, "BIN_XOFF", "0", "offset in pixels in x direction to apply before binning"},
54 {ARG_INT, "BIN_YOFF", "0", "offset in pixels in y direction to apply before binning"},
55
56 {ARG_FLOAT, "GAUSS_SIG", "-1.0", "width of gaussian"},
57 {ARG_INT, "GAUSS_WID", "-1", "half width of kernel"},
58 {ARG_INT, "GAUSS_NSUB", "1", "distance between sampled points in pixels"},
59 {ARG_INT, "GAUSS_XOFF", "0", "offset in pixels to add to x index of input image"},
60 {ARG_INT, "GAUSS_YOFF", "0", "offset in pixels to add to y index of input image"},
61
62 /*
63 {ARG_INT, "IHORIZ", "0", "flag for outputting horizontal component only, otherwise output radial component only"},
64 tplarson 1.1 {ARG_INT, "IBOX", "0", "flag for convolving psf with box"},
65 {ARG_INT, "IGAUSS", "0", "flag for convolving psf with gaussian"},
66 {ARG_FLOAT, "WPSFX", "1.0", "width of gaussian psf in x direction"},
67 {ARG_FLOAT, "WPSFY", "1.0", "width of gaussian psf in y direction"},
68 {ARG_INT, "ILORENTZ", "0", "flag for convolving psf with lorentzian"},
69 {ARG_FLOAT, "WLORENTZ", "1.0", "width of lorenztian psf"},
70 {ARG_INT, "ILININT", "0", "flag for convolving psf with linear interpolation"},
71 {ARG_INT, "ICUBINT", "0", "flag for convolving psf with cubic interpolation"},
72 {ARG_INT, "IBIN", "0", "flag for binning"},
73 {ARG_INT, "NBIN", "1", "factor for binning"},
74 {ARG_INT, "ISKIP", "0", "flag for subsampling"},
75 {ARG_INT, "NOFF", "2", "offset for subsampling"},
76 */
77 {ARG_TIME, "TSTART", "1996.05.26_16:00:00_TAI", "first output time"},
78 {ARG_INT, "DTOFF", "0", "number of timesteps to offset first image"},
79 {ARG_INT, "NDT", "0", "number of time points to generate"},
80
81 {ARG_END}
82 };
83
84 #include "saveparm.c"
85 tplarson 1.1 #include "timing.c"
86 #include "set_history.c"
87 #include "setplm2.c"
88
89 static void Ccker(double *u, double s);
90 static double Ccintd(double *f, int nx, double x);
91
92 int DoIt(void)
93 {
94 int newstat=0;
95 int status = DRMS_SUCCESS;
96 DRMS_Record_t *outrec = NULL;
97 DRMS_Segment_t *segout = NULL;
98 DRMS_Array_t *outarr = NULL;
99 DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
100 DRMS_RecLifetime_t lifetime;
101 long long histrecnum=-1;
102 int length[2];
103 TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
104 int ipsf=0;
105 struct fresize_struct fress;
|
106 tplarson 1.2 static double defaultcoefs[] = {1.0, 0.459224, 0.132395, 0.019601, 0.000802, -4.31934E-05 };
|
107 tplarson 1.1
108 int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
109 int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);
110
111 /*
112 double wt0, wt1, wt2, wt3, wt;
113 double ut0, ut1, ut2, ut3, ut;
114 double st0, st1, st2, st3, st;
115 double ct0, ct1, ct2, ct3, ct;
116
117 wt0=getwalltime();
118 ct0=getcputime(&ut0, &st0);
119 */
120
121 char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
122 int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
123 int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
124 if (permflag)
125 lifetime = DRMS_PERMANENT;
126 else
127 lifetime = DRMS_TRANSIENT;
128 tplarson 1.1
129 char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);
130
131 char *modefile = (char *)cmdparams_save_str(&cmdparams, "MODELIST", &newstat);
132 int xpixels = cmdparams_save_int(&cmdparams, "XPIXELS", &newstat);
133 int ypixels = cmdparams_save_int(&cmdparams, "YPIXELS", &newstat);
134 double solarp = cmdparams_save_float(&cmdparams, "SOLARP", &newstat);
135 double obsb0 = cmdparams_save_float(&cmdparams, "OBSB0", &newstat);
136 double obsdist = cmdparams_save_float(&cmdparams, "OBSDIST", &newstat);
137 float xoffset = cmdparams_save_float(&cmdparams, "XOFFSET", &newstat);
138 float yoffset = cmdparams_save_float(&cmdparams, "YOFFSET", &newstat);
139 int iinten = cmdparams_save_int(&cmdparams, "IINTEN", &newstat);
|
140 tplarson 1.2 float phasein = cmdparams_save_float(&cmdparams, "PHASE", &newstat);
141 float freqin = cmdparams_save_float(&cmdparams, "FREQUENCY", &newstat);
142 int ilimbdark = cmdparams_save_int(&cmdparams, "ILIMBDARK", &newstat);
143 double coefs[6];
144 int n_user_coefs = cmdparams_get_int(&cmdparams, "coefs_nvals", &status);
145 if (ilimbdark)
146 {
147 if (n_user_coefs == 6)
148 {
149 double *cmdcoefs;
150 int i;
151 cmdparams_get_dblarr(&cmdparams, "coefs", &cmdcoefs, &status);
152 for (i=0; i<6; i++)
153 coefs[i] = cmdcoefs[i];
154 }
155 else
156 {
157 int i;
158 for (i=0; i<6; i++)
159 coefs[i] = defaultcoefs[i];
160 }
161 tplarson 1.2 }
162 int iheight = cmdparams_save_int(&cmdparams, "IHEIGHT", &newstat);
|
163 tplarson 1.1
164 int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
165 int bxoff = cmdparams_save_int(&cmdparams, "BIN_XOFF", &newstat);
166 int byoff = cmdparams_save_int(&cmdparams, "BIN_YOFF", &newstat);
167
168 float sigma = cmdparams_save_float(&cmdparams, "GAUSS_SIG", &newstat);
169 int hwidth = cmdparams_save_int(&cmdparams, "GAUSS_WID", &newstat);
170 int nsub = cmdparams_save_int(&cmdparams, "GAUSS_NSUB", &newstat);
171 int gxoff = cmdparams_save_int(&cmdparams, "GAUSS_XOFF", &newstat);
172 int gyoff = cmdparams_save_int(&cmdparams, "GAUSS_YOFF", &newstat);
173 // for now use this patch so i can specify integers for GAUSS_SIG
174 sigma/=sqrt(2);
175
176 /*
177 int ibox = cmdparams_save_int(&cmdparams, "IBOX", &newstat);
178 int igauss = cmdparams_save_int(&cmdparams, "IGAUSS", &newstat);
179 float wpsfx = cmdparams_save_float(&cmdparams, "WPSFX", &newstat);
180 float wpsfy = cmdparams_save_float(&cmdparams, "WPSFY", &newstat);
181 int ilorentz = cmdparams_save_int(&cmdparams, "ILORENTZ", &newstat);
182 float wlorentz = cmdparams_save_float(&cmdparams, "WLORENTZ", &newstat);
183 int ilinint = cmdparams_save_int(&cmdparams, "ILININT", &newstat);
184 tplarson 1.1 int icubint = cmdparams_save_int(&cmdparams, "ICUBINT", &newstat);
185 int ibin = cmdparams_save_int(&cmdparams, "IBIN", &newstat);
186 int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
187 int iskip = cmdparams_save_int(&cmdparams, "ISKIP", &newstat);
188 int noff = cmdparams_save_int(&cmdparams, "NOFF", &newstat);
189 */
190 double tstart = cmdparams_save_time(&cmdparams, "TSTART", &newstat);
191 int dtoff = cmdparams_save_int(&cmdparams, "DTOFF", &newstat);
192 int ndt = cmdparams_save_int(&cmdparams, "NDT", &newstat);
193
194 double pangle = RADSINDEG*solarp;
195 double b0 = RADSINDEG*obsb0;
196
197 // if (ibox+igauss+ilorentz+ilinint+icubint > 0) ipsf=1;
198
199 FILE *fptr = fopen(modefile,"r");
200 if (fptr == NULL)
201 {
202 fprintf(stderr, "ERROR: can't open mode file %s\n",modefile);
203 return 1;
204 }
205 tplarson 1.1
206 if (newstat)
207 {
208 fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
209 cpsave_decode_error(newstat);
210 return 1;
211 }
212 else if (savestrlen != strlen(savestr))
213 {
214 fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
215 return 1;
216 }
217
218 // set up ancillary dataseries for processing metadata
219 DRMS_Record_t *tempoutrec = drms_create_record(drms_env,
220 outseries,
221 DRMS_TRANSIENT,
222 &status);
223
224 if (status != DRMS_SUCCESS)
225 {
226 tplarson 1.1 fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
227 return 1;
228 }
229
230 DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
|
231 tplarson 1.5 // char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/mkylms.c,v 1.4 2013/04/28 09:14:04 tplarson Exp $");
|
232 tplarson 1.1 if (histlink != NULL)
233 {
|
234 tplarson 1.3 histrecnum=set_history(histlink);
|
235 tplarson 1.1 if (histrecnum < 0)
236 {
237 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
238 return 1;
239 }
240 }
241 else
242 {
243 fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
244 }
245
246
247 // these must be present in the output dataseries and variable, not links or constants
248 char *outchecklist[] = {"T_REC", "QUALITY", "T_OBS", "CRLT_OBS", "CRLN_OBS",
249 //"SAT_ROT", "INST_ROT", "IM_SCALE",
250 "CDELT1", "CDELT2"};
251
252 int itest;
253 for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
254 {
255 DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
256 tplarson 1.1 if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)
257 {
258 fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
259 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
260 return 1;
261 }
262 }
263 float cadence = drms_getkey_float(tempoutrec, "T_REC_step", &status);
264 drms_close_record(tempoutrec, DRMS_FREE_RECORD);
265
266 int i, hold;
267 long nmodes=0;
268 while ((hold = getc(fptr)) != EOF)
269 if (hold == '\n')
270 nmodes++;
271 fclose(fptr);
272
273 fptr = fopen(modefile,"r");
274 int *lim, *mim;
275 int lin, min;
276 int mmax=0;
277 tplarson 1.1 int lmax=0;
278 int *marr_lmax, *marr_count, **marr_index;
279 // fscanf(fptr,"%d",&nmodes);
280 lim=(int *)malloc(nmodes*sizeof(int));
281 mim=(int *)malloc(nmodes*sizeof(int));
282 for (i=0;i<nmodes;i++)
283 {
284 if (fscanf(fptr,"%d %d\n", &lin, &min) != 2)
285 {
286 fprintf(stderr,"ERROR: something went wrong on line i=%d\n", i);
287 return 1;
288 }
289 // fscanf(fptr,"%d",&lin);
290 // fscanf(fptr,"%d",&min);
291 if (min > lin || lin < 0)
292 {
293 fprintf(stderr,"ERROR: mode file has m > l or l < 0\n");
294 return 1;
295 }
296 /*
297 if (min > mmax)
298 tplarson 1.1 mmax=min;
299 if (lin > lmax)
300 lmax=lin;
301 */
302 lim[i]=lin;
303 mim[i]=min;
304 }
305 fclose(fptr);
306
307 if (verbflag)
308 printf("nmodes=%ld, dtoff=%d, ndt=%d\n", nmodes, dtoff, ndt);
309
310 int *lim2, *mim2;
311 if (ndt > 0)
312 {
313 if (dtoff + ndt > nmodes)
314 nmodes=nmodes-dtoff;
315 else
316 nmodes=ndt;
317
318 lim2=(int *)malloc(nmodes*sizeof(int));
319 tplarson 1.1 mim2=(int *)malloc(nmodes*sizeof(int));
320 for (i=0;i<nmodes;i++)
321 {
322 lim2[i]=lim[dtoff+i];
323 mim2[i]=mim[dtoff+i];
324 }
325 free(lim);
326 free(mim);
327 lim=lim2;
328 mim=mim2;
329 }
330
331 if (nmodes <= 0)
332 {
333 fprintf(stderr, "ERROR: no modes selected\n");
334 return 1;
335 }
336
337 for (i=0;i<nmodes;i++)
338 {
339 if (mim[i] > mmax)
340 tplarson 1.1 mmax=mim[i];
341 if (lim[i] > lmax)
342 lmax=lim[i];
343 }
344
345 marr_lmax=(int *)(malloc((mmax+1)*sizeof(int)));
346 marr_count=(int *)(malloc((mmax+1)*sizeof(int)));
347 marr_index=(int **)(malloc((mmax+1)*sizeof(int *)));
348 for (i=0;i<=mmax;i++)
349 {
350 marr_lmax[i]=-1;
351 marr_count[i]=0;
352 marr_index[i]=(int *)(malloc((lmax+1)*sizeof(int)));
353 }
354 for (i=0;i<nmodes;i++)
355 {
356 if (lim[i] > marr_lmax[mim[i]])
357 marr_lmax[mim[i]]=lim[i];
358 marr_index[mim[i]][marr_count[mim[i]]]=i;
359 marr_count[mim[i]]++;
360 }
361 tplarson 1.1
362 // double obsdist = 1.0; //in au
363 double distobs = obsdist*AU/RTRUE; //in solar radii
364 double dsunobs = obsdist*AU; //in meters
365 double x0=(double)xpixels/2-0.5 + xoffset;
366 double y0=(double)ypixels/2-0.5 + yoffset;
367 long npixels=xpixels*ypixels;
368 double imscale=1.97784*1024/xpixels;
369 double scale=imscale/(180*60*60/PI);
370 length[0]=xpixels;
371 length[1]=ypixels;
372 float obsl0=211.67;
373 int obscr=1784;
|
374 tplarson 1.3 // double rsunobs=atan(asin(RTRUE/AU/obsdist))*60*60/RADSINDEG;
375 // double rsun=rsunobs/imscale;
376 double rsunobs=asin(RTRUE/AU/obsdist)*60*60/RADSINDEG;
377 double rsun=tan(asin(RTRUE/AU/obsdist))/scale;
378 double cdelt=rsunobs/rsun; // Use mean image scale across solar image
|
379 tplarson 1.1 /*
380 ; Sun is sitting at the center of the main coordinate system
381 ; and has radius 1.
382 ; Observer is at robs=(xobs,yobs,zobs) moving with a velocity
383 ; vobs.
384 ; Start by setting robs from distobs and b0
385 */
386
387 double robs_x = distobs * cos(b0);
388 double robs_y = 0.0;
389 double robs_z = distobs * sin(b0);
390
391 // Start with CCD coordinates
392 double *x6 = (double *)(malloc(npixels*sizeof(double)));
393 double *y6 = (double *)(malloc(npixels*sizeof(double)));
394 int j;
395 for (i=0;i<ypixels;i++)
396 for (j=0;j<xpixels;j++)
397 x6[i*ypixels+j]=(double)j;
398 for (i=0;i<ypixels;i++)
399 for (j=0;j<xpixels;j++)
400 tplarson 1.1 y6[i*ypixels+j]=(double)i;
401
402 double *x4 = (double *)(malloc(npixels*sizeof(double)));
403 double *y4 = (double *)(malloc(npixels*sizeof(double)));
404 for (i=0;i<npixels;i++)
405 {
406 x4[i]=scale*(x6[i] - x0);
407 y4[i]=scale*(y6[i] - y0);
408 }
409
410
411 /*
412 ; Rotate by the P-angle. New coordinate system has the y-axis pointing
413 ; towards the solar north pole.
414 */
415
416 double *x2 = (double *)(malloc(npixels*sizeof(double)));
417 double *y2 = (double *)(malloc(npixels*sizeof(double)));
418 for (i=0;i<npixels;i++)
419 {
420 x2[i]= x4[i]*cos(pangle) + y4[i]*sin(pangle);
421 tplarson 1.1 y2[i]=-x4[i]*sin(pangle) + y4[i]*cos(pangle);
422 }
423
424
425 /*
426 ; Now transform to put the coordinates into the solar coordinate
427 ; system.
428 ; First find the directions (vecx and vecy) the x2 and y2 coordinate
429 ; axis correspond to. vecsun points towards the Sun. Note that the
430 ; (x2,y2,Sun) system is left handed. These vectors are unit vectors.
431 */
432 double vecx_x=0.0;
433 double vecx_y=1.0;
434 double vecx_z=0.0;
435 double vecy_x=-sin(b0);
436 double vecy_y=0.0;
437 double vecy_z=cos(b0);
438 double vecsun_x=-cos(b0);
439 double vecsun_y=0.0;
440 double vecsun_z=-sin(b0);
441 // Now the proper direction can be found. These are not unit vectors.
442 tplarson 1.1 double *x1 = (double *)(malloc(npixels*sizeof(double)));
443 double *y1 = (double *)(malloc(npixels*sizeof(double)));
444 double *z1 = (double *)(malloc(npixels*sizeof(double)));
445 double *q1 = (double *)(malloc(npixels*sizeof(double)));
446 for (i=0;i<npixels;i++)
447 {
448 x1[i]=vecx_x*x2[i] + vecy_x*y2[i] + vecsun_x;
449 y1[i]=vecx_y*x2[i] + vecy_y*y2[i] + vecsun_y;
450 z1[i]=vecx_z*x2[i] + vecy_z*y2[i] + vecsun_z;
451 q1[i]=1/sqrt(x1[i]*x1[i] + y1[i]*y1[i] + z1[i]*z1[i]);
452 }
453 // Make them unit vectors.
454 for (i=0;i<npixels;i++)
455 {
456 x1[i]*=q1[i];
457 y1[i]*=q1[i];
458 z1[i]*=q1[i];
459 }
460
461
462
463 tplarson 1.1 // Now find intersection with the Sun.
|
464 tplarson 1.2 // Solve quadratic equation |robs+q*[x1,y1,z1]|=1 for q
465 // a, b and c are terms in a*x^2+bx+c=0. a==1 since [x1,y1,z1] is unit vector.
466 // double *b = (double *)(malloc(npixels*sizeof(double)));
467 // double *d = (double *)(malloc(npixels*sizeof(double)));
468 double b,d;
469 double *q = (double *)(malloc(npixels*sizeof(double)));
|
470 tplarson 1.1 int *inflag = (int *)(malloc(npixels*sizeof(int)));
471 double c = robs_x*robs_x + robs_y*robs_y + robs_z*robs_z -1;
|
472 tplarson 1.2 double *xsun = (double *)(malloc(npixels*sizeof(double)));
473 double *ysun = (double *)(malloc(npixels*sizeof(double)));
474 double *zsun = (double *)(malloc(npixels*sizeof(double)));
475 double *cosrho = (double *)(malloc(npixels*sizeof(double)));
|
476 tplarson 1.1 for (i=0;i<npixels;i++)
477 {
|
478 tplarson 1.2 // b[i]=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
479 // d[i]=b[i]*b[i] - 4*c;
480 b=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
481 d=b*b - 4*c;
482 if (d >= 0)
483 {
|
484 tplarson 1.1 inflag[i]=1;
|
485 tplarson 1.2 q[i]=(-b - sqrt(d))/2;
486 xsun[i]=robs_x + x1[i]*q[i];
487 ysun[i]=robs_y + y1[i]*q[i];
488 zsun[i]=robs_z + z1[i]*q[i];
489 cosrho[i]=-(x1[i]*xsun[i]+y1[i]*ysun[i]+z1[i]*zsun[i]);
490 }
|
491 tplarson 1.1 else
|
492 tplarson 1.2 {
|
493 tplarson 1.1 inflag[i]=0;
|
494 tplarson 1.2 xsun[i]=0;
495 ysun[i]=0;
496 zsun[i]=0;
497 cosrho[i]=0;
498 }
499 }
500
501 if (iheight)
502 {
503 // Now calculate observing height in units of the solar radius
504 double pc3[]={0.24047045,-0.76504325,0.86252178,-0.33859142};
505 double height;
506 double c0 = robs_x*robs_x + robs_y*robs_y + robs_z*robs_z;
507 for (i=0;i<npixels;i++)
508 {
509 double x=1.0;
510 height=pc3[0];
511 for (j=1;j<4;j++)
512 {
513 x*=cosrho[i];
514 height+=pc3[j]*x;
515 tplarson 1.2 }
516 height*=1e6/RTRUE;
517 // Recalculate coordinates.
518 b=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
519 c=c0-(1+height)*(1+height);
520 d=b*b - 4*c;
521 if (d >= 0)
522 {
523 inflag[i]=1;
524 q[i]=(-b - sqrt(d))/2;
525 // Remember to divide by new target radius
526 xsun[i]=(robs_x + x1[i]*q[i])/(1+height);
527 ysun[i]=(robs_y + y1[i]*q[i])/(1+height);
528 zsun[i]=(robs_z + z1[i]*q[i])/(1+height);
529 }
530 else
531 {
532 inflag[i]=0;
533 xsun[i]=0;
534 ysun[i]=0;
535 zsun[i]=0;
536 tplarson 1.2 }
537 }
|
538 tplarson 1.1 }
539
540 double *phisun = (double *)(malloc(npixels*sizeof(double)));
541 double *cphisun = (double *)(malloc(npixels*sizeof(double)));
542 double *sphisun = (double *)(malloc(npixels*sizeof(double)));
543 double *clatsun = (double *)(malloc(npixels*sizeof(double)));
544 double *slatsun = (double *)(malloc(npixels*sizeof(double)));
545 double *prad = (double *)(malloc(npixels*sizeof(double)));
546 double *pphi = (double *)(malloc(npixels*sizeof(double)));
547 double *plat = (double *)(malloc(npixels*sizeof(double)));
548
549 for (i=0;i<npixels;i++)
550 {
551 phisun[i] = atan2(ysun[i],xsun[i]);
552 cphisun[i] = cos(phisun[i]);
553 sphisun[i] = sin(phisun[i]);
554 slatsun[i] = zsun[i];
555 clatsun[i] = sqrt(1 - zsun[i]*zsun[i]);
556 // Set velocity projection factors
557 prad[i]=xsun[i]*x1[i]+ysun[i]*y1[i]+zsun[i]*z1[i];
558 pphi[i]=-sphisun[i]*x1[i]+cphisun[i]*y1[i];
559 tplarson 1.1 plat[i]=-cphisun[i]*slatsun[i]*x1[i]-sphisun[i]*slatsun[i]*y1[i]+clatsun[i]*z1[i];
|
560 tplarson 1.2 }
561
562 double *ld;
563 if (ilimbdark)
564 {
565 ld=(double *)malloc(npixels*sizeof(double));
566 int iy,ix;
567 for (iy=0; iy<ypixels; iy++)
568 for (ix=0; ix<xpixels; ix++)
569 {
570 double costheta2;
571 double xi, mu, z, ldi;
572 double x, y, R2;
573 int ord;
574 int index=iy*xpixels + ix;
575
576 if (!inflag[index])
577 continue;
578
579 /* get coordinates of point */
580 x = ((double)ix - x0) / rsun;
581 tplarson 1.2 y = ((double)iy - y0) / rsun;
582
583
584 R2 = x*x + y*y;
585 if (R2 <= 1.0)
586 costheta2 = 1.0 - R2;
587 else
588 costheta2 = 0.0;
|
589 tplarson 1.1
|
590 tplarson 1.2 mu = sqrt(costheta2);
591 xi = log(mu);
592 z = 1.0;
593 ldi = 1.0;
594 for (ord=1; ord<6; ord++)
595 {
596 z *= xi;
597 ldi += coefs[ord] * z;
598 }
599 // not sure what to do with this...
600 if (ldi <= 0.0 || isnan(ldi))
601 {
602 //*Op = missval;
603 //ncropped++;
604 ld[index]=0.0;
605 }
606 else
607 ld[index] = ldi;
608 }
|
609 tplarson 1.1 }
610
611 free(x6);free(y6);free(x4);free(y4);free(x2);free(y2);
612 free(x1);free(y1);free(z1);free(q1);
|
613 tplarson 1.2 // free(b);free(d);
614 free(cosrho);
|
615 tplarson 1.1 free(xsun);free(ysun);free(zsun);
616
617 if (verbflag)
618 printf("coordinate arrays set up\n");
619
620 int lmax1, l, m, il;
621 long nplm=10001;
622 double *plm = (double *)(malloc(nplm*(lmax+1)*sizeof(double)));
623 double *dplm = (double *)(malloc(nplm*(lmax+1)*sizeof(double)));
624 double *plms = (double *)(malloc(nplm*nmodes*sizeof(double)));
625 double *dplms = (double *)(malloc(nplm*nmodes*sizeof(double)));
626 double *dxplm = (double *)(malloc(nplm*sizeof(double)));
627 for (i=0;i<nplm;i++)
628 dxplm[i]= -1 + (double)i * 2 / (nplm - 1);
629 int *indx = (int *)(malloc ((lmax+1)*sizeof(int)));
630 for (l=0; l<=lmax; l++)
631 indx[l]=l;
632
633 for (m=0;m<=mmax;m++)
634 {
635 if (marr_lmax[m] == -1)
636 tplarson 1.1 continue;
637 lmax1=marr_lmax[m];
638 setplm2(m, lmax1, m, nplm, indx, dxplm, nplm, plm, dplm);
639 for (i=0;i<marr_count[m];i++)
640 {
641 il=marr_index[m][i];
642 l=lim[il];
643 for (j=0;j<nplm;j++)
644 {
645 plms[il*nplm + j]=plm[l*nplm + j];
646 dplms[il*nplm + j]=dplm[l*nplm + j];
647 }
648 }
649 }
650
651 free(plm);
652 free(dplm);
653 free(dxplm);
654
655 if (verbflag)
656 printf("plms set up\n");
657 tplarson 1.1
658 // float cadence=60.0;
659 // double time0=612201600.0 + dtoff*cadence;
660 double time0=tstart + dtoff*cadence;
661 double *vradr=(double *)(malloc(npixels*sizeof(double)));
662 double *vradi=(double *)(malloc(npixels*sizeof(double)));
663
664 //reuse arrays malloc'ed above to save memory
665 double *vlatr, *vlati;
666 vlatr=vradr;
667 vlati=vradi;
668 // double *vlatr=(double *)(malloc(npixels*sizeof(double)));
669 // double *vlati=(double *)(malloc(npixels*sizeof(double)));
670 double *vphir=(double *)(malloc(npixels*sizeof(double)));
671 double *vphii=(double *)(malloc(npixels*sizeof(double)));
672
673 /*
674 double *velr=(double *)(malloc(npixels*sizeof(double)));
675 double *veli=(double *)(malloc(npixels*sizeof(double)));
676 double *velsum=(double *)(malloc(npixels*sizeof(double)));
677 */
678 tplarson 1.1 float *velr=(float *)(malloc(npixels*sizeof(float)));
679 float *veli=(float *)(malloc(npixels*sizeof(float)));
680 float *velsum=(float *)(malloc(npixels*sizeof(float)));
681 float *velrsave=velr;
682 float *velisave=veli;
683 float *velssave=velsum;
684
685 double *plminterp=(double *)(malloc(npixels*sizeof(double)));
686 double ll, dinterp;
687
688 float *binptr;
689 int xpix1=xpixels;
690 int ypix1=ypixels;
691 if (nbin > 0)
692 {
693 x0=(x0 - bxoff + 0.5)/nbin - 0.5;
694 y0=(y0 - byoff + 0.5)/nbin - 0.5;
|
695 tplarson 1.3 // imscale*=nbin;
696 cdelt*=nbin;
|
697 tplarson 1.1 xpix1=xpixels/nbin;
698 ypix1=ypixels/nbin;
699 length[0]=xpix1;
700 length[1]=ypix1;
701 binptr=(float *)malloc(xpix1*ypix1*sizeof(float));
702 }
703
704 float *gaussptr;
705 int xpix2, ypix2;
706 if (hwidth > 0)
707 {
708 x0=(x0 - gxoff)/nsub;
709 y0=(y0 - gyoff)/nsub;
|
710 tplarson 1.3 // imscale*=nsub;
711 cdelt*=nsub;
|
712 tplarson 1.1 xpix2=xpix1/nsub;
713 ypix2=ypix1/nsub;
714 length[0]=xpix2;
715 length[1]=ypix2;
716 gaussptr=(float *)malloc(xpix2*ypix2*sizeof(float));
717 init_fresize_gaussian(&fress, sigma, hwidth, nsub);
718 }
|
719 tplarson 1.2
720 double *phase=(double *)malloc(npixels*sizeof(double));
721 double rsecs = RTRUE/C;
722 for (i=0;i<npixels;i++)
723 {
724 if (inflag[i])
725 {
726 double deltat;
727 deltat=(q[i]-distobs+1)*rsecs;
728 phase[i]=phasein-2*PI*freqin*deltat;
729 }
730 }
731
|
732 tplarson 1.1 for (i=0;i<nmodes;i++)
733 {
734 trec=time0+i*cadence;
735 if (verbflag > 1)
736 printf("i=%d, trec=%f, l=%d, m=%d\n",i,trec,lim[i],mim[i]);
737
738 outrec = drms_create_record(drms_env, outseries, lifetime, &status);
739 if (status != DRMS_SUCCESS || outrec == NULL)
740 {
741 fprintf(stderr,"ERROR: unable to open record in output dataseries, i = %d, status = %d\n", i, status);
742 return 0;
743 }
744 if (histlink)
745 drms_setlink_static(outrec, histlinkname, histrecnum);
746
747 for (j=0;j<npixels;j++)
748 plminterp[j]=Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
749
750 for (j=0;j<npixels;j++)
751 {
752 if (inflag[j])
753 tplarson 1.1 {
|
754 tplarson 1.2 vradr[j]=1000 * (cos(mim[i]*phisun[j]+phase[j])) * plminterp[j]; //Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
755 vradi[j]=1000 * (sin(mim[i]*phisun[j]+phase[j])) * plminterp[j]; //Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
|
756 tplarson 1.1 // needed to make main leaks positive since DATASIGN=-1
757 if (iinten)
758 {
|
759 tplarson 1.2 if (ilimbdark)
760 {
761 velr[j]=-1*vradr[j]*ld[j];
762 veli[j]=-1*vradi[j]*ld[j];
763 velsum[j]=velr[j]+veli[j];
764 }
765 else
766 {
767 velr[j]=-1*vradr[j];
768 veli[j]=-1*vradi[j];
769 velsum[j]=velr[j]+veli[j];
770 }
|
771 tplarson 1.1 }
772 else
773 {
774 velr[j]=prad[j]*vradr[j];
775 veli[j]=prad[j]*vradi[j];
|
776 tplarson 1.2 velsum[j]=velr[j]+veli[j];
|
777 tplarson 1.1 }
778 /*
779 if (j == 524800)
780 velsum[j]=1.0;
781 else
782 velsum[j]=0.0;
783 */
784 }
785 else
786 {
787 velr[j]=0;
788 veli[j]=0;
789 velsum[j]=0;
790 }
791 }
792
793 if (nbin > 0)
794 {
795 fbin(velr, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
796 velr=binptr;
797 }
798 tplarson 1.1
799 if (hwidth > 0)
800 {
801 fresize(&fress, velr, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
802 velr=gaussptr;
803 }
804
805 segout = drms_segment_lookup(outrec, "vradr");
806 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velr, &status);
807 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
808 {
809 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
810 return 0;
811 }
812 drms_segment_write(segout, outarr, 0);
813 free(outarr);
814
815 if (nbin > 0)
816 {
817 fbin(veli, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
818 veli=binptr;
819 tplarson 1.1 }
820
821 if (hwidth > 0)
822 {
823 fresize(&fress, veli, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
824 veli=gaussptr;
825 }
826
827 segout = drms_segment_lookup(outrec, "vradi");
828 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, veli, &status);
829 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
830 {
831 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
832 return 0;
833 }
834 drms_segment_write(segout, outarr, 0);
835 free(outarr);
836
837 if (nbin > 0)
838 {
839 fbin(velsum, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
840 tplarson 1.1 velsum=binptr;
841 }
842
843 if (hwidth > 0)
844 {
845 fresize(&fress, velsum, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
846 velsum=gaussptr;
847 }
848
849 segout = drms_segment_lookup(outrec, "vradsum");
850 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velsum, &status);
851 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
852 {
853 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
854 return 0;
855 }
856 drms_segment_write(segout, outarr, 0);
857 free(outarr);
858
859 if (!iinten)
860 {
861 tplarson 1.1 velr=velrsave;
862 veli=velisave;
863 velsum=velssave;
864 for (j=0;j<npixels;j++)
865 {
866 if (inflag[j])
867 {
868 if (lim[i] != 0)
869 {
870 ll=sqrt(lim[i]*(lim[i]+1.0));
871 dinterp=Ccintd(&dplms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
|
872 tplarson 1.2 vlatr[j]=1000 * (cos(mim[i]*phisun[j]+phase[j]))*dinterp/ll*clatsun[j];
873 vlati[j]=1000 * (sin(mim[i]*phisun[j]+phase[j]))*dinterp/ll*clatsun[j];
874 vphir[j]=1000 * mim[i]*(-sin(mim[i]*phisun[j]+phase[j]))*plminterp[j]/ll/clatsun[j];
875 vphii[j]=1000 * mim[i]*( cos(mim[i]*phisun[j]+phase[j]))*plminterp[j]/ll/clatsun[j];
|
876 tplarson 1.1 }
877 else
878 {
879 vlatr[j]=vlati[j]=0;
880 vphir[j]=vphii[j]=0;
881 }
882 velr[j]=pphi[j]*vphir[j]+plat[j]*vlatr[j];
883 veli[j]=pphi[j]*vphii[j]+plat[j]*vlati[j];
|
884 tplarson 1.2 velsum[j]=velr[j]+veli[j];
|
885 tplarson 1.1 }
886 else
887 {
888 velr[j]=0;
889 veli[j]=0;
890 velsum[j]=0;
891 }
892 }
893
894 if (nbin > 0)
895 {
896 fbin(velr, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
897 velr=binptr;
898 }
899
900 if (hwidth > 0)
901 {
902 fresize(&fress, velr, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
903 velr=gaussptr;
904 }
905
906 tplarson 1.1 segout = drms_segment_lookup(outrec, "vhorr");
907 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velr, &status);
908 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
909 {
910 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
911 return 0;
912 }
913 drms_segment_write(segout, outarr, 0);
914 free(outarr);
915
916 if (nbin > 0)
917 {
918 fbin(veli, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
919 veli=binptr;
920 }
921
922 if (hwidth > 0)
923 {
924 fresize(&fress, veli, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
925 veli=gaussptr;
926 }
927 tplarson 1.1
928 segout = drms_segment_lookup(outrec, "vhori");
929 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, veli, &status);
930 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
931 {
932 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
933 return 0;
934 }
935 drms_segment_write(segout, outarr, 0);
936 free(outarr);
937
938 if (nbin > 0)
939 {
940 fbin(velsum, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
941 velsum=binptr;
942 }
943
944 if (hwidth > 0)
945 {
946 fresize(&fress, velsum, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
947 velsum=gaussptr;
948 tplarson 1.1 }
949
950 segout = drms_segment_lookup(outrec, "vhorsum");
951 outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velsum, &status);
952 if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
953 {
954 fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
955 return 0;
956 }
957 drms_segment_write(segout, outarr, 0);
958 free(outarr);
959
960 }
961
962 velr=velrsave;
963 veli=velisave;
964 velsum=velssave;
965
966 drms_setkey_time(outrec, "T_REC", trec);
967 drms_setkey_int(outrec, "L", lim[i]);
968 drms_setkey_int(outrec, "M", mim[i]);
969 tplarson 1.1 drms_setkey_float(outrec, "XOFF", xoffset);
970 drms_setkey_float(outrec, "YOFF", yoffset);
971 drms_setkey_float(outrec, "PANGLE", solarp);
972 drms_setkey_float(outrec, "BANGLE", obsb0);
973 //obviously these are redundant, but give series designer freedom to use either
974 drms_setkey_float(outrec, "DISTOBS", distobs);
975 drms_setkey_float(outrec, "OBSDIST", obsdist);
976
977 drms_setkey_string(outrec, "FILE", modefile);
978
979 drms_setkey_float(outrec, "CRPIX1", x0+1);
980 drms_setkey_float(outrec, "CRVAL1", 0.0);
|
981 tplarson 1.3 drms_setkey_float(outrec, "CDELT1", cdelt);
|
982 tplarson 1.1
983 drms_setkey_float(outrec, "CRPIX2", y0+1);
984 drms_setkey_float(outrec, "CRVAL2", 0.0);
985 drms_setkey_float(outrec, "CROTA2", -solarp);
|
986 tplarson 1.3 drms_setkey_float(outrec, "CDELT2", cdelt);
|
987 tplarson 1.1
988 drms_setkey_time(outrec, "T_OBS", trec);
989 drms_setkey_int(outrec, "QUALITY", 0);
990 drms_setkey_float(outrec, "CADENCE", cadence);
991 drms_setkey_float(outrec, "CRLT_OBS", obsb0);
992 drms_setkey_float(outrec, "CRLN_OBS", obsl0);
993 drms_setkey_int(outrec, "CAR_ROT", obscr);
994
|
995 tplarson 1.5 drms_setkey_double(outrec, "DSUN_OBS", dsunobs);
996 drms_setkey_double(outrec, "RSUN_OBS", rsunobs);
|
997 tplarson 1.1
|
998 tplarson 1.5 drms_setkey_double(outrec, "OBS_VR", 0.0);
999 drms_setkey_double(outrec, "OBS_VW", 0.0);
1000 drms_setkey_double(outrec, "OBS_VN", 0.0);
|
1001 tplarson 1.1
1002 drms_setkey_int(outrec, "DATASIGN", -1);
1003
1004 tnow = (double)time(NULL);
1005 tnow += UNIX_epoch;
1006 drms_setkey_time(outrec, "DATE", tnow);
1007 drms_close_record(outrec, DRMS_INSERT_RECORD);
1008
1009 }
1010
1011 if (hwidth > 0)
1012 free_fresize(&fress);
1013
1014 if (verbflag)
1015 printf("module %s successful completion\n", cmdparams.argv[0]);
1016
1017 return 0;
1018 }
1019
1020
1021
1022 tplarson 1.1 /* Calculate the interpolation kernel. */
1023 void Ccker(double *u, double s)
1024 {
1025 double s2, s3;
1026
1027 s2= s * s;
1028 s3= s2 * s;
1029 u[0] = s2 - 0.5 * (s3 + s);
1030 u[1] = 1.5*s3 - 2.5*s2 + 1.0;
1031 u[2] = -1.5*s3 + 2.0*s2 + 0.5*s;
1032 u[3] = 0.5 * (s3 - s2);
1033 }
1034
1035 /* Cubic convolution interpolation */
1036 double Ccintd(double *f, int nx, double x)
1037 {
1038 double ux[4];
1039 /* Sum changed to double to speed up things */
1040 double sum;
1041
1042 int ix, ix1, i;
1043 tplarson 1.1
1044 if (x < 1. || x >= (double)(nx-2))
1045 return 0.0;
1046
1047 ix = (int)x;
1048 Ccker(ux, x - (double)ix);
1049
1050 ix1 = ix - 1;
1051 sum = 0.;
1052 for (i = 0; i < 4; i++)
1053 sum = sum + f[ix1 + i] * ux[i];
1054 return sum;
1055 }
1056
1057
|
Return to
mkylms.c
CVS log
Up to [Development] / JSOC / proj / globalhs / apps