1 arta 1.1 /*
2 I.Scholl
3
|
4 scholl 1.6 limbfit(): adapted from Solar Astrometry Program (solarSDO.c)
5 Marcelo Emilio - Jan 2010
|
6 arta 1.1
7
8 #define CODE_NAME "limbfit"
|
9 scholl 1.17 #define CODE_VERSION "V4.0r9"
10 #define CODE_DATE "Mon Aug 31 17:16:24 PDT 2015"
|
11 arta 1.1 */
12
13 #include "limbfit.h"
14
|
15 scholl 1.6 #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
16
|
17 scholl 1.11 void sav_b0(float *pf_sb0, float *pl_sb0, float *pf_b0)
18 {
19 float *p_sb0=pf_sb0;
20 float *p_b0=pf_b0;
21 while(p_sb0<=pl_sb0) *(p_sb0++)=*(p_b0++);
22 }
23
24 void sum_b0(float *beta, float *pf_b0, float *pl_b0)
25 {
26 float *p_b0=pf_b0;
27 float *p_beta=beta;
28 while(p_b0 <= pl_b0) *(p_b0++)=*(p_b0)+*(p_beta++);
29 }
30
|
31 scholl 1.13 int limbfit(LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios)
|
32 arta 1.1 {
|
33 scholl 1.4 static char *log_msg_code="limbfit";
|
34 scholl 1.10 char log_msg[200];
|
35 scholl 1.4
|
36 scholl 1.13 if (results->debug) lf_logmsg("DEBUG", "APP", 0, 0, "", log_msg_code, results->opf);
|
37 arta 1.1
38 /************************************************************************
39 INIT VARS coming from build_lev1.c
40 ************************************************************************/
|
41 scholl 1.6 float *data=input->data;
|
42 arta 1.1 // float MIN_VALUE=BAD_PIXEL_VALUE;
43 // float MAX_VALUE=(float)fabs(BAD_PIXEL_VALUE);
|
44 scholl 1.2
|
45 arta 1.1
46 int ret_code = 0;
47
|
48 scholl 1.2 int w = ANNULUS_WIDTH;
49 long S = MAX_SIZE_ANN_VARS;
|
50 scholl 1.11 int nang = NUM_LDF; //180
51 int nprf = NUM_RADIAL_BINS; //64
|
52 scholl 1.2 int nreg = NUM_AB_BINS;
53 float rsi = LO_LIMIT;
54 float rso = UP_LIMIT;
55 int r_size = NUM_FITPNTS;
56 int grange = GUESS_RANGE;
57 float dx = INC_X;
58 float dy = INC_Y;
|
59 scholl 1.6 int naxis_row= input->img_sz0;
60 int naxis_col= input->img_sz1;
|
61 scholl 1.7 float lahi = AHI;
|
62 scholl 1.11 int fldf;
|
63 scholl 1.10 //int skip = SKIPGC;
|
64 scholl 1.2 //float flag = BAD_PIXEL_VALUE;
|
65 scholl 1.9 int iter = NB_ITER;
66 if (input->iter!=NB_ITER) iter=input->iter;
|
67 scholl 1.11 fldf=input->fldf;
|
68 scholl 1.7
69 if (input->spe==1)
70 {
71 iter=NB_ITER2;
72 lahi=AHI2;
73 rsi=LO_LIMIT2;
74 rso=UP_LIMIT2;
75 }
|
76 arta 1.1
77 /************************************************************************/
78 /* set parameters */
79 /************************************************************************/
|
80 scholl 1.11 long npixels=naxis_row*naxis_col;
|
81 scholl 1.6 long ii, jj, jk, i,j;
|
82 scholl 1.2 float cmx, cmy,r;//, guess_cx, guess_cy, guess_r;
|
83 scholl 1.11 int nitr=0, ncut=0;
|
84 arta 1.1
85 r = (float)naxis_row/2;
|
86 scholl 1.4 float sav_max=0.;
|
87 arta 1.1
|
88 scholl 1.10 /* Initial guess estimate of center position from Richard's code*/
89 cmx=(float)input->ix;
90 cmy=(float)input->iy;
91 r=(float)input->ir;
|
92 arta 1.1
93
94 /************************************************************************/
95 /* Make annulus data */
96 /************************************************************************/
|
97 scholl 1.13 int nbc=3;
98 float *p_anls=ios->pf_anls;
|
99 arta 1.1
|
100 scholl 1.13 if (ios->is_firstobs == 0)
|
101 scholl 1.6 {
|
102 scholl 1.13 ios->is_firstobs=1;
103 float d;
|
104 scholl 1.16 float w2p=(float)r+w/2;
105 float w2m=(float)r-w/2;
|
106 scholl 1.13 double iimcmy2,jjmcmx;
107 /* Select Points */
108 jk=-1;
109
110 for(ii = 0; ii < naxis_row; ii++)
111 {
112 iimcmy2=(ii-cmy)*(ii-cmy);
113 for(jj = 0; jj < naxis_col; jj++)
114 {
115 jjmcmx=jj-cmx;
116 d=(float)sqrt(iimcmy2+(jjmcmx)*(jjmcmx));
117 if (d<=w2p && d>=w2m)
118 {
119 jk++;
120 *(p_anls++)=(float)jj;
121 *(p_anls++)=(float)ii;
122 *(p_anls++)=data[ii*naxis_col+jj];
123 }
|
124 arta 1.1 }
125 }
|
126 scholl 1.13 ios->anls_nbpix=jk;
127 ios->pl_anls=p_anls-1;
128
129 if (results->debug)
130 {
131 sprintf(log_msg," jk = %ld", jk);
132 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
133 }
134 if ((jk*3) >= S)
135 {
136 lf_logmsg("ERROR", "APP", ERR_SIZE_ANN_TOO_BIG, 0,"nbc>S", log_msg_code, results->opf);
137 return ERR_SIZE_ANN_TOO_BIG;
138 }
|
139 arta 1.1 }
|
140 scholl 1.13 else
|
141 arta 1.1 {
|
142 scholl 1.13 jk=ios->anls_nbpix;
143 ii=0;
144 p_anls=p_anls+2;
145 while(p_anls<=ios->pl_anls)
146 {
147 *(p_anls)=data[(int)ios->anls[nbc*ii+1]*naxis_col+(int)ios->anls[nbc*ii]];
148 ii++;
149 p_anls=p_anls+3;
150 }
|
151 scholl 1.11 }
|
152 arta 1.1
153 /************************************************************************/
154 /* Call Fortran Code Subrotine limb.f */
155 /************************************************************************/
|
156 scholl 1.13 float *rprf, *lprf, *alph, *beta, *b0, *sb0; //, *beta1, *beta2, *beta3;
157 long ab_nrow=nreg, ab_ncol;
158 rprf = (float *) malloc(sizeof(float)*(nprf));
159 if(!rprf)
160 {
161 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (rprf)", log_msg_code, results->opf);
162 return ERR_MALLOC_FAILED;
163 }
164 float *p_rprf, *pl_rprf;
165
166 lprf = (float *) malloc(sizeof(float)*((nang+1)*nprf));
167 if(!lprf)
168 {
169 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (lprf)", log_msg_code, results->opf);
170 return ERR_MALLOC_FAILED;
171 }
172 float *p_lprf, *pl_lprf;
173
174 alph = (float *) malloc(sizeof(float)*(nreg));
175 if(!alph)
176 {
177 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (alph)", log_msg_code, results->opf);
178 return ERR_MALLOC_FAILED;
179 }
180 float *pf_alph=&alph[0];
181 float *p_alph;
182 float *pl_alph=&alph[ab_nrow-1];
183
184 beta = (float *) malloc(sizeof(float)*(nreg));
185 if(!beta)
186 {
187 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (xbeta)", log_msg_code, results->opf);
188 return ERR_MALLOC_FAILED;
189 }
190 float *pf_beta=&beta[0];
191 float *p_beta;
192
193 b0 = (float *) malloc(sizeof(float)*(nreg));
194 if(!b0)
195 {
196 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (b0)", log_msg_code, results->opf);
197 return ERR_MALLOC_FAILED;
198 scholl 1.13 }
199 float *pf_b0=&b0[0];
200 float *p_b0;
201 float *pl_b0=&b0[ab_nrow-1];
202
203 sb0 = (float *) malloc(sizeof(float)*(nreg));
204 if(!sb0)
205 {
206 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (sb0)", log_msg_code, results->opf);
207 return ERR_MALLOC_FAILED;
208 }
209 float *pf_sb0=&sb0[0];
210 float *p_sb0;
211 float *pl_sb0=&sb0[ab_nrow-1];
212
213 p_b0=pf_b0;
214 while(p_b0<=pl_b0) *(p_b0++)=0.;
|
215 arta 1.1
216 // fortran call:
|
217 scholl 1.13 if (results->debug) lf_logmsg("DEBUG", "APP", 0, 0, "entering limb", log_msg_code, results->opf);
|
218 scholl 1.4
|
219 scholl 1.9 int centyp=0; //=0 do center calculation, =1 skip center calculation
|
220 arta 1.1 //#1
|
221 scholl 1.9 if (input->cc == 0) centyp=1;
|
222 scholl 1.11
223 //in limb_ call: beta contains the output beta, b0 is initialized with 0.0
|
224 scholl 1.13 int it=1, ifail=0;
|
225 scholl 1.11
226 do
|
227 scholl 1.4 {
|
228 scholl 1.13 // init everything that is passed to fortran
229 //
230 p_alph=pf_alph;
231 p_beta=pf_beta;
232 while(p_alph<=pl_alph)
233 {
234 *(p_beta++)=0.;
235 *(p_alph++)=0.;
236 }
237 p_rprf=&rprf[0];
238 pl_rprf=&rprf[nprf-1];
239 while(p_rprf<=pl_rprf) *(p_rprf++)=0.;
240 p_lprf=&lprf[0];
241 pl_lprf=&lprf[((nang+1)*nprf)-1];
242 while(p_lprf<=pl_lprf) *(p_lprf++)=0.;
243
244 if (results->debug)
|
245 scholl 1.4 {
|
246 scholl 1.11 sprintf(log_msg,"entering limb# %d", it);
|
247 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
248 scholl 1.4 }
|
249 scholl 1.11 ifail=0;
|
250 scholl 1.13 limb_(&ios->anls[0],&jk, &cmx, &cmy, &r, &nitr, &ncut, &rprf[0], &lprf[0], &rsi, &rso,
|
251 scholl 1.11 &dx, &dy, &alph[0], &beta[0], &ifail, &b0[0], ¢yp, &lahi);
|
252 scholl 1.13 if (results->debug)
|
253 scholl 1.11 {
|
254 scholl 1.13 sprintf(log_msg,"exiting limb ifail= %d", ifail);
255 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
256 arta 1.1 }
|
257 scholl 1.11 if(ifail==0)
258 {
259 centyp=1;
260 if(it==iter && it>1) sav_b0(pf_sb0,pl_sb0,pf_b0);
261 sum_b0(&beta[0],pf_b0,pl_b0);
262 }
263 it++;
|
264 arta 1.1 }
|
265 scholl 1.11 while(it<=iter && ifail==0);
266
267
|
268 arta 1.1
|
269 scholl 1.13 if (results->debug)
|
270 arta 1.1 {
|
271 scholl 1.4 sprintf(log_msg," nitr = %6d", nitr);
|
272 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
273 scholl 1.10 sprintf(log_msg," cmx = %6.2f, cmy = %6.2f", cmx, cmy);
|
274 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
275 arta 1.1 }
276
277 if(ifail == 0) // || ifail == 7)
278 {
279 /************************************************************************/
280 /* Compute the mean of the center of the image */
281 /************************************************************************/
282 double cmean,ctot=0.0;
283 long nbp=0;
284 float limx_m=cmx-250;
285 float limx_p=cmx+250;
286 float limy_m=cmy-250;
287 float limy_p=cmy+250;
288 for (i=limx_m;i<limx_p;i++)
289 {
290 for (j=limy_m;j<limy_p;j++)
291 {
292 ctot=ctot+data[i*naxis_col+j];
293 nbp++;
294 }
295 }
296 arta 1.1 cmean=ctot/nbp;
|
297 scholl 1.13 if (results->debug)
|
298 arta 1.1 {
|
299 scholl 1.4 sprintf(log_msg," cmean = %6.4f (ctot= %6.4f , nbp=%ld)", cmean,ctot,nbp);
|
300 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
301 arta 1.1 }
302
303 /************************************************************************/
304 /* Compute the Inflection Point */
305 /************************************************************************/
306 float *LDF, *D;
307 LDF = (float *) malloc(sizeof(float)*(nprf));
|
308 scholl 1.6 if(!LDF)
309 {
|
310 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (LDF)", log_msg_code, results->opf);
|
311 scholl 1.6 return ERR_MALLOC_FAILED;
312 }
313
|
314 arta 1.1 D = (float *) malloc(sizeof(float)*(nprf));
|
315 scholl 1.6 if(!D)
316 {
|
317 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (D)", log_msg_code, results->opf);
|
318 scholl 1.6 return ERR_MALLOC_FAILED;
319 }
|
320 arta 1.1
321 double ip, ip1, maxim;
|
322 scholl 1.11 double radius = 0.;
|
323 scholl 1.4 int cont, c, ret_gsl;
|
324 arta 1.1 float h;
325
326 // for saving them in FITS file
|
327 scholl 1.6 float *save_ldf, *save_alpha_beta;
328 double *save_params;
|
329 scholl 1.3 long nb_p_as=6, nb_p_es=6, nb_p_radius=1, nb_p_ip=1;
330 long params_nrow=nang+1, params_ncol=nb_p_as+nb_p_es+nb_p_radius+nb_p_ip;
|
331 arta 1.1 long ldf_nrow=nang+2, ldf_ncol=nprf; //ii -nbr of ldfs, jj -nbr of points for each ldf
|
332 scholl 1.11 if (iter > 1) ab_ncol=3; else ab_ncol=2;
|
333 scholl 1.15 int degf=6;
|
334 scholl 1.13 //if (iter > 2) degf=4; else degf=6; // problem with degf = 4...
335 double A[degf];
336 double erro[degf];
|
337 scholl 1.2
|
338 arta 1.1 save_ldf = (float *) malloc((ldf_nrow*ldf_ncol)*sizeof(float));
|
339 scholl 1.6 if(!save_ldf)
340 {
|
341 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_ldf)", log_msg_code, results->opf);
|
342 scholl 1.6 return ERR_MALLOC_FAILED;
343 }
344 save_alpha_beta = (float *) malloc((ab_nrow*ab_ncol)*sizeof(float));
345 if(!save_alpha_beta)
346 {
|
347 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_alpha_beta)", log_msg_code, results->opf);
|
348 scholl 1.6 return ERR_MALLOC_FAILED;
349 }
350 save_params = (double *) malloc((params_nrow*params_ncol)*sizeof(double));
351 if(!save_params)
352 {
|
353 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_params)", log_msg_code, results->opf);
|
354 scholl 1.6 return ERR_MALLOC_FAILED;
355 }
|
356 scholl 1.3
357 long zero_as=0;
358 long zero_es=params_nrow*nb_p_as;
359 long zero_r =params_nrow*(nb_p_as+nb_p_es);
360 long zero_ip=params_nrow*(nb_p_as+nb_p_es+nb_p_radius);
|
361 scholl 1.11
|
362 scholl 1.13 /* copy the last pixel of the mean ldf in the prev one to correct a fortran issue */
|
363 scholl 1.11 lprf[(nprf*nang)+63]=lprf[(nprf*nang)+62];
|
364 arta 1.1
|
365 scholl 1.17 for (cont=0; cont<=nang; cont++)
366 { // if only the mean ldf replace above line by below line
367 // cont=nang;
|
368 scholl 1.4 ret_gsl=0;
|
369 arta 1.1 ip1=0.;
370 /* dx */
371 h=(float)rprf[1]-rprf[0];
372
373 /* Take the last (mean) LDF from lprf vector */
|
374 scholl 1.13 for(ii = 0; ii < nprf ; ii++) LDF[ii]=lprf[(nprf*cont)+ii];
|
375 scholl 1.2
|
376 arta 1.1 /* Calculate the Derivative */
377 D[0]=(-3*LDF[4]+16*LDF[3]-36*LDF[2]+48*LDF[1]-25*LDF[0])/(12*h);
378 D[1]=(LDF[4]-6*LDF[3]+18*LDF[2]-10*LDF[1]-3*LDF[0])/(12*h);
|
379 scholl 1.13 for (i=2; i< nprf-2; i++) D[i]=(LDF[i-2]-8*LDF[i-1]+8*LDF[i+1]-LDF[i+2])/(12*h);
|
380 arta 1.1 D[nprf-2]=(-LDF[nprf-5]+6*LDF[nprf-4]-18*LDF[nprf-3]+10*LDF[nprf-2]+3*LDF[nprf-1])/(12*h);
381 D[nprf-1]=(3*LDF[nprf-5]-16*LDF[nprf-4]+36*LDF[nprf-3]-48*LDF[nprf-2]+25*LDF[nprf-1])/(12*h);
|
382 scholl 1.13
|
383 arta 1.1 /* square the result */
|
384 scholl 1.11 // for(ii = 0; ii < nprf ; ii++)
385 // D[ii]=D[ii]*D[ii]; //pointers here!
386 float* p_D=&D[0];
387 float* pl_D=&D[nprf-1];
388 while(p_D<=pl_D) *(p_D++)=(*(p_D) * *(p_D));
389
|
390 arta 1.1 /* find the maximum */
391 jj=0;
392 maxim=-1;
393 for(ii = 1; ii < nprf-1 ; ii++)
394 {
395 if (D[ii] > maxim)
396 {
397 maxim=(double)D[ii];
398 jj=ii;
399 }
400 }
401
402 /* improve the maximum estimation looking at the 2 neibors points */
403 ip=rprf[jj]-h/2.*(D[jj-1]-D[jj+1])/(2*D[jj]-D[jj-1]-D[jj+1]);
404
|
405 scholl 1.13 if (results->debug)
|
406 arta 1.1 {
407 if (cont==nang)
|
408 scholl 1.4 {
409 sprintf(log_msg," Inflection Point 1: %ld %d %8.5f %8.5f", jj, nprf, rprf[jj], ip);
|
410 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
411 scholl 1.4 }
|
412 arta 1.1 }
413 ip1=ip;
414 /* FIT A GAUSSIAN PLUS QUADRATIC FUNCTION */
|
415 scholl 1.13
|
416 arta 1.1 /* Select Region */
417 int m_in, m_ex, N;
418 /* Gaussian + quadratic function */
419 /* After Launch verify if this number is OK */
420 /* Should contain the pixels around the 3*FWHM level */
421 m_in=jj-r_size;
422 m_ex=jj+r_size;
423
424 if (m_in < 0) m_in = 0;
425 if (m_ex > nprf-1) m_ex =nprf-1;
426 N=m_ex-m_in+1;
|
427 scholl 1.13
|
428 arta 1.1 /* parameters:- initial guess */
429 A[0]= maxim;
430 A[1]= ip;
431 A[2]= 0.5;
432 A[3]= 12*maxim;
|
433 scholl 1.13 for (c=0;c<degf;c++) erro[c]=0.;
|
434 arta 1.1
|
435 scholl 1.13 if (N >= degf) // degree of freedom 6 parameters, 6 values minimum
|
436 arta 1.1 {
437 double px[N], der[N] , sigma[N];
438
439 jj=-1;
440 for(ii = m_in; ii <= m_ex ; ii++)
441 {
442 jj++;
443 px[jj]=rprf[ii];
444 der[jj]=D[ii];
445 sigma[jj] = 1.;
446 }
|
447 scholl 1.13 //A[4]= der[N-1]-der[0];
|
448 scholl 1.4 //if (debug==2) fprintf(opf,"%s DEBUG_INFO in limbfit: #: %d\n",LOGMSG1,cont);
|
449 scholl 1.13
450 ret_gsl=gaussfit(der, px, sigma, A, erro, N, degf,results->debug,results->opf);
451
|
452 scholl 1.4 if (ret_gsl < 0)
|
453 arta 1.1 {
|
454 scholl 1.4 if (cont==nang)
455 {
456 sprintf(log_msg," gaussfit failed for the averaged LDF %d err:%d", cont,ret_gsl);
|
457 scholl 1.13 lf_logmsg("ERROR", "APP", 0, ret_gsl, log_msg, log_msg_code, results->opf);
|
458 scholl 1.4 }
|
459 scholl 1.13 for (c=0;c<degf;c++)
|
460 arta 1.1 {
461 A[c]=0.;
462 erro[c]=0.;
463 }
464 radius=0.;
465 }
466 else
467 {
468 /* FIND THE MAXIMUM OF THE GAUSSIAN PLUS QUADRATIC FUNCTION */
469 radius = A[1]; /* Initial Guess */
|
470 scholl 1.16 radius = fin_min(A, radius, grange, results->debug,results->opf);
|
471 arta 1.1 if (radius < 0)
472 {
|
473 scholl 1.4 ret_gsl=(int)radius;
474 if (cont==nang)
475 {
476 sprintf(log_msg," fin_min failed for the averaged LDF %d err:%d", cont, ret_gsl);
|
477 scholl 1.13 lf_logmsg("ERROR", "APP", 0, ret_gsl, log_msg, log_msg_code, results->opf);
|
478 scholl 1.4 }
|
479 scholl 1.13 for (c=0;c<degf;c++)
|
480 arta 1.1 {
481 A[c]=0.;
482 erro[c]=0.;
483 }
484 radius=0.;
485 }
486 }
487 }
488 else
489 {
|
490 scholl 1.4 sprintf(log_msg," Inflection point too close to annulus data border %d", cont);
|
491 scholl 1.13 lf_logmsg("WARNING", "APP", 0, 0, log_msg, log_msg_code, results->opf);
492 for (c=0;c<degf;c++) erro[c]=0.;
|
493 arta 1.1 radius=ip;
|
494 scholl 1.6 save_params[zero_ip+cont]=0.;
|
495 arta 1.1 }
496 // save them
|
497 scholl 1.13 for (c=0;c<degf;c++)
|
498 arta 1.1 {
|
499 scholl 1.6 save_params[zero_as+params_nrow*c+cont]=A[c];
500 save_params[zero_es+params_nrow*c+cont]=erro[c];
501 }
502 save_params[zero_r+cont]=radius;
503 save_params[zero_ip+cont]=ip1;
|
504 arta 1.1 } // endfor-cont
|
505 scholl 1.13 if (results->debug)
|
506 arta 1.1 {
|
507 scholl 1.4 sprintf(log_msg," Inflection Point 2: %8.5f %8.5f %8.5f", A[1], erro[1], radius);
|
508 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
509 if (degf==4)
510 {
511 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f ", A[0],A[1],A[2],A[3]);
512 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
513 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f ", erro[0],erro[1],erro[2],erro[3]);
514 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
515 }
516 else
517 {
518 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f", A[0],A[1],A[2],A[3],A[4],A[5]);
519 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
520 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f", erro[0],erro[1],erro[2],erro[3],erro[4],erro[5]);
521 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
522 }
|
523 arta 1.1 }
524
525 //**********************************************************************
526 // Save LDFs & AB data in a FITS file
527 //**********************************************************************
|
528 scholl 1.11
529 int fldfr=0; //proc result: 0=ok; 1=failed; 2=cannot be processed; 3=malloc pb; 4=not processed;
530 int fulldf_nrows, fulldf_ncols=2;
531 int bins1=0, bins2=0;
532 int retcode=0;
533 float *save_full_ldf;
534 if (fldf==1)
535 {
536 // full LDF
537
538 //if (ret_gsl<0 || cmx < 0.01 || cmy < 0.01 || radius < 0.01) // that I don't remember what for
539 if (cmx < 0.01 || cmy < 0.01 || radius < 0.01)
540 {
541 fulldf_nrows=1;
542 bins1=0;
543 save_full_ldf = (float *) malloc((2)*sizeof(float));
544 if(!save_full_ldf)
545 {
|
546 scholl 1.13 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_full_ldf)", log_msg_code, results->opf);
|
547 scholl 1.11 return ERR_MALLOC_FAILED;
548 }
549 save_full_ldf[0]=0;
550 save_full_ldf[1]=0;
551 fldfr=3;
552 }
553 else
554 {
|
555 scholl 1.13 retcode=mk_fldfs(cmx, cmy, radius, naxis_row, naxis_col, npixels, data, &save_full_ldf, &bins1, &bins2, results->opf, results->debug);
|
556 scholl 1.11 if (retcode == ERR_MALLOC_FAILED)
557 {
558 fldfr=2;
559 return ERR_MALLOC_FAILED;
560 }
561 else
562 if (retcode == ERR_NR_STACK_TOO_SMALL)
563 {
564 ret_code=ERR_LIMBFIT_FLDF_FAILED;
565 fldfr=1;
566 //shouldnt exist?
567 }
568 fulldf_nrows=bins2;
569 }
570 }
|
571 scholl 1.6
572 // LDF // lprf & rprf come from limbfit.f
|
573 arta 1.1 float *p_sldf=&save_ldf[0];
574 float *pl_sldf=&save_ldf[(ldf_nrow*ldf_ncol)-1];
|
575 scholl 1.13 p_rprf=&rprf[0];
576 pl_rprf=&rprf[nprf-1];
|
577 arta 1.1 while (p_rprf <= pl_rprf)
578 *(p_sldf++)=*(p_rprf++);
579
|
580 scholl 1.13 p_lprf=&lprf[0];
|
581 arta 1.1 p_sldf=&save_ldf[ldf_ncol];
582 while (p_sldf <= pl_sldf)
583 *(p_sldf++)=*(p_lprf++);
584
|
585 scholl 1.6 // AB
|
586 scholl 1.13 p_alph=pf_alph;
|
587 scholl 1.11 p_b0=pf_b0;
|
588 scholl 1.6 float *p_save_alpha_beta=&save_alpha_beta[0];
|
589 scholl 1.3 while (p_alph <= pl_alph) *(p_save_alpha_beta++)=*(p_alph++);
|
590 scholl 1.11 while (p_b0 <= pl_b0) *(p_save_alpha_beta++)=*(p_b0++);
591 if (iter>1)
592 {
593 p_sb0=pf_sb0;
594 while (p_sb0 <= pl_sb0) *(p_save_alpha_beta++)=*(p_sb0++);
595 }
|
596 scholl 1.3
|
597 arta 1.1 // Update Returned Structure when process succeeded
598 results->cenx=cmx;
599 results->ceny=cmy;
600 results->radius=radius;
601 results->cmean=cmean;
|
602 scholl 1.4 results->max_limb=sav_max;
|
603 scholl 1.3 results->numext=3;
|
604 scholl 1.2 results->fits_ldfs_naxis1=ldf_ncol;
605 results->fits_ldfs_naxis2=ldf_nrow;
|
606 scholl 1.11 results->fits_fldfs_tfields=fulldf_ncols;
607 results->fits_fldfs_nrows=fulldf_nrows;
|
608 scholl 1.3 results->fits_ab_tfields=ab_ncol;
|
609 scholl 1.2 results->fits_ab_nrows=ab_nrow;
|
610 scholl 1.3 results->fits_params_tfields=params_ncol;
611 results->fits_params_nrows=params_nrow;
|
612 scholl 1.11 results->nb_fbins=bins1;
|
613 scholl 1.2 results->ann_wd=w;
614 results->mxszannv=S;
615 results->nb_ldf=nang;
616 results->nb_rdb=nprf;
617 results->nb_abb=nreg;
618 results->up_limit=rsi;
619 results->lo_limit=rso;
620 results->inc_x=dx;
621 results->inc_y=dy;
622 results->nfitpnts=r_size;
623 results->nb_iter=iter;
|
624 scholl 1.11 results->fldfr=fldfr;
|
625 scholl 1.7 results->ahi=lahi;
|
626 scholl 1.4 results->error1=ifail;
627 results->error2=ret_gsl;
|
628 scholl 1.9 results->cc=input->cc;
|
629 arta 1.1 if (ifail == 0)
630 results->quality=9;
|
631 scholl 1.6 //? if (cont>=nang && ret_gsl<0)
632 if (ret_gsl<0)
|
633 scholl 1.4 {
|
634 scholl 1.6 results->quality=1;
|
635 scholl 1.4 ret_code=ERR_LIMBFIT_FIT_FAILED;
|
636 scholl 1.13 if (results->debug)
|
637 scholl 1.4 {
638 sprintf(log_msg," ret_gsl<0 = %2d", ret_gsl);
|
639 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
640 scholl 1.4 }
641 }
|
642 arta 1.1 results->fits_ldfs_data=save_ldf;
|
643 scholl 1.6 results->fits_params=save_params;
644 results->fits_alpha_beta=save_alpha_beta;
|
645 scholl 1.11 if (fldf == 1) results->fits_fulldfs=save_full_ldf; else results->fits_fulldfs=0;
|
646 arta 1.1 free(D);
647 free(LDF);
|
648 scholl 1.4 } // end limb OK
|
649 arta 1.1 else
650 {
|
651 scholl 1.13 if (results->debug)
|
652 arta 1.1 {
|
653 scholl 1.4 sprintf(log_msg," limb.f routine returned ifail = %2d", ifail);
|
654 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
655 arta 1.1 }
656 // Update Returned Structure when process failed
657 results->numext=0;
|
658 scholl 1.4 results->error1=ifail;
659 results->error2=-1;
|
660 arta 1.1 results->quality=0;
661 ret_code=ERR_LIMBFIT_FAILED;
|
662 scholl 1.2 results->ann_wd=w;
663 results->mxszannv=S;
664 results->nb_ldf=nang;
665 results->nb_rdb=nprf;
666 results->nb_abb=nreg;
667 results->up_limit=rsi;
668 results->lo_limit=rso;
669 results->inc_x=dx;
670 results->inc_y=dy;
671 results->nfitpnts=r_size;
672 results->nb_iter=iter;
|
673 scholl 1.7 results->ahi=lahi;
|
674 scholl 1.2 results->max_limb=0;
675 results->cmean=0;
|
676 scholl 1.6 results->nb_fbins=0;
|
677 scholl 1.11 results->fldfr=4;
|
678 scholl 1.9 results->cc=input->cc;
|
679 scholl 1.4 } // end limb failed
|
680 scholl 1.11 // IS: do not free those (save_ldf,save_params,save_alpha_beta,save_full_ldf) passed from or to the structure !
|
681 scholl 1.14 free(rprf);
682 free(alph);
683 free(beta);
684 free(b0);
685 free(sb0);
686 free(lprf);
|
687 scholl 1.13 if (results->debug)
|
688 arta 1.1 {
|
689 scholl 1.4 sprintf(log_msg," >>>>end of limbfit with: %d", ret_code);
|
690 scholl 1.13 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
691 arta 1.1 }
|
692 scholl 1.4 sprintf(log_msg," end: RC: %d - limb:%d - fit:%d - quality: %d", ret_code, results->error1, results->error2, results->quality);
|
693 scholl 1.13 lf_logmsg("INFO", "APP", 0, 0, log_msg, log_msg_code, results->opf);
|
694 scholl 1.6
|
695 arta 1.1 return(ret_code);
696
697 } /* end of main */
698
699
700 /*--------------------------------------------------------------------------*/
|
701 scholl 1.13 int gaussfit(double y[], double t[],double sigma[], double A[], double erro[], long N, int degf, int debug, FILE *opf)
|
702 arta 1.1 /* Calculate a Least SqrareGaussian + Quadratic fit */
703 /* Uses the GNU Scientific Library */
704 /* Marcelo Emilio (c) v 1.0 Jan 2009 */
705 /* fits A[0]exp(-z^2/2)+A[3]+A[4]*x+A[5]*x^2 */
706 /* z=(t-A[1])/A[2] */
707 /* Need to add : */
708 /* #include <gsl/gsl_rng.h> */
709 /* #include <gsl/gsl_randist.h> */
710 /* #include <gsl/gsl_vector.h> */
711 /* #include <gsl/gsl_blas.h> */
712 /* #include <gsl/gsl_multifit_nlin.h> */
713 /* #include "expfit.c" */
714 /* compiles as gcc program.c --lm -lgsl -lgslcblas */
715 /* y --> f(x) - ordinate values */
716 /* t --> x(i) - abscissa values */
717 /* sigma --> independent gaussian erros */
718 /* N --> Number of points in the vector */
719 /* A --> Input a initial guess and output the result */
720 /* erro --> Chi-square of A */
721
722 {
|
723 scholl 1.4 static char *log_msg_code="gaussfit";
|
724 arta 1.1 int ret_code=0;
|
725 scholl 1.10 char log_msg[200];
|
726 scholl 1.13
|
727 arta 1.1 const gsl_multifit_fdfsolver_type *T;
728 gsl_multifit_fdfsolver *s;
729 int status;
730 unsigned int i, iter = 0;
731 const size_t n = N;
|
732 scholl 1.13 const size_t p = degf;
|
733 arta 1.1
734 gsl_matrix *covar = gsl_matrix_alloc (p, p);
735 /* double t[N], y[N], sigma[N]; */
736 struct dataF d = { n, t, y, sigma};
737 gsl_multifit_function_fdf f;
738
739 /* Initial Guess */
|
740 scholl 1.13 double x_init[degf];
741 for (i=0; i <degf; i++)
|
742 arta 1.1 x_init[i]=A[i];
743
744 gsl_vector_view x = gsl_vector_view_array (x_init, p);
745 const gsl_rng_type * type;
746 gsl_rng * r;
747
748 gsl_rng_env_setup();
749
750 type = gsl_rng_default;
751 r = gsl_rng_alloc (type);
752
753 f.f = &expb_f;
754 f.df = &expb_df;
755 f.fdf = &expb_fdf;
756 f.n = n;
757 f.p = p;
758 f.params = &d;
759
760 T = gsl_multifit_fdfsolver_lmsder;
761 s = gsl_multifit_fdfsolver_alloc (T, n, p);
762 gsl_multifit_fdfsolver_set (s, &f, &x.vector);
763 arta 1.1
764 do
765 {
766 iter++;
767 status = gsl_multifit_fdfsolver_iterate (s);
|
768 scholl 1.13
769 if (status==0 || status == GSL_ENOPROG)
770 status = gsl_multifit_test_delta (s->dx, s->x,1e-4, 1e-4); /* IS: change epsrel (last arg of gsl_multifit_test_delta) after the SDO Launch */
|
771 arta 1.1 else
|
772 scholl 1.4 {
|
773 scholl 1.13 ret_code=ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED;
|
774 scholl 1.4 if (debug)
775 {
|
776 scholl 1.13 sprintf(log_msg," iter: %u, status: %d (%s)\n", iter,status,gsl_strerror (status));
|
777 scholl 1.4 lf_logmsg("WARNING", "APP", ret_code, status, log_msg, log_msg_code, opf);
778 }
779 }
|
780 scholl 1.13
|
781 arta 1.1 }
782 while (status == GSL_CONTINUE && iter < 500);
|
783 scholl 1.13
|
784 arta 1.1
|
785 scholl 1.13 if (status == 0)
|
786 arta 1.1 {
|
787 scholl 1.13
|
788 arta 1.1 gsl_multifit_covar (s->J, 0.0, covar);
789
790 #define FIT(i) gsl_vector_get(s->x, i)
791 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
792
793 double chi = gsl_blas_dnrm2(s->f);
794 double dof = (double)(n - p);
795 double c = GSL_MAX_DBL(1, chi / sqrt(dof));
796
797
|
798 scholl 1.13 for (i=0; i <degf; i++) {
|
799 arta 1.1 A[i]=FIT(i);
800 erro[i]=c*ERR(i);
801 }
802 }
803 else
804 {
805 ret_code=ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED;
|
806 scholl 1.4 if (debug)
807 {
808 sprintf(log_msg," (gsl_multifit_fdfsolver_iterate) failed, (iter=%u) gsl_errno=%d", iter,status);
809 lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);
810 }
|
811 arta 1.1 }
812
813
814 if (debug==2)
815 {
|
816 scholl 1.4 sprintf(log_msg," Nonlinear LSQ fitting status = %s (%d)", gsl_strerror (status),status);
817 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);
|
818 arta 1.1 }
|
819 scholl 1.13 if (debug)
820 {
821 sprintf(log_msg," finished at iter# %u\n", iter);
822 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);
823 }
|
824 arta 1.1 gsl_multifit_fdfsolver_free (s);
825 gsl_matrix_free (covar);
826 gsl_rng_free (r);
827 return(ret_code);
828 }
829
830
831 /*--------------------------------------------------------------------------*/
832
|
833 scholl 1.16 double fin_min(double A[], double m, int range, int debug, FILE *opf)
|
834 arta 1.1 /* Calculate the maximum of the quadratic + gaussian function */
835 /* using Brent algorithm */
836 /* Marcelo Emilio (c) v 1.0 Jan 2009 */
837 /* A are the parameters of the gaussian */
838 /* m is the guess where the function is maximum */
839 /* Calculate the maximum around */
840 /* (-5 + m) < m > (5 - m) ( in pixel units ) */
841 /* Returns the m vaule where the function is maximum */
842 /* requires: #include "expmax.c" --> definition of the function */
843 /* #include <gsl/gsl_errno.h> */
844 /* #include <gsl/gsl_math.h> */
845 /* #include <gsl/gsl_min.h> */
846 {
|
847 scholl 1.4 static char *log_msg_code="fin_min";
|
848 arta 1.1 int status;
849 gsl_set_error_handler_off();
|
850 scholl 1.4 int ret_code=0;
|
851 scholl 1.10 char log_msg[200];
|
852 arta 1.1
853 int iter = 0, max_iter = 1000;
854 const gsl_min_fminimizer_type *T;
855 gsl_min_fminimizer *s;
856
857 gsl_function F;
858
859 //double m_expected = m+0.01; //1e-8;
860 //double a = m-5000, b = m+5000;
861 double a = m-range, b = m+range; // initially = 2
|
862 scholl 1.13 // if (degf == 4)
863 //struct exp_plus_quadratic_params params= { A[0], A[1], A[2], A[3] };
864 // else
|
865 arta 1.1 struct exp_plus_quadratic_params params= { A[0], A[1], A[2], A[3], A[4], A[5] };
866 F.function = &exp_plus_quadratic_function;
867 F.params = ¶ms;
868
869 T = gsl_min_fminimizer_brent;
870 s = gsl_min_fminimizer_alloc (T);
871
872 status=gsl_min_fminimizer_set (s, &F, m, a, b);
873 if (status)
874 {
|
875 scholl 1.4 if (status == GSL_EINVAL)
876 {
|
877 scholl 1.6 ret_code=ERR_GSL_FINMIN_NOMIN_FAILED;
|
878 scholl 1.4 if (debug)
879 {
880 sprintf(log_msg," (gsl_min_fminimizer_set) doesn't find a minimum, m=%f", m);
881 lf_logmsg("DEBUG", "APP", 0, status, log_msg, log_msg_code, opf);
882 }
883 }
884 else
885 {
886 ret_code=ERR_GSL_FINMIN_SET_FAILED;
|
887 scholl 1.6 if (debug)
888 {
889 sprintf(log_msg," (gsl_min_fminimizer_set) failed, gsl_errno=%d", status);
890 lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);
891 }
|
892 arta 1.1 }
|
893 scholl 1.11 gsl_min_fminimizer_free(s);
|
894 scholl 1.4 return ret_code;
|
895 arta 1.1 }
896
897 do
898 {
899 iter++;
900 status = gsl_min_fminimizer_iterate (s);
901
902 m = gsl_min_fminimizer_x_minimum (s);
903 a = gsl_min_fminimizer_x_lower (s);
904 b = gsl_min_fminimizer_x_upper (s);
905
|
906 scholl 1.13 status = gsl_min_test_interval (a, b,1E-3, 0.0);
|
907 arta 1.1 }
908 while (status == GSL_CONTINUE && iter < max_iter);
909
910 if (status) // regarder lequel plantait et verifier si on a autre chose que GSL_EINVAL comme erreur
911 { // ajouter debug pour les messages apres
|
912 scholl 1.4 if (status == GSL_EINVAL)
913 {
914 ret_code=ERR_GSL_FINMIN_PRO_FAILED;
915 if (debug)
916 {
917 sprintf(log_msg," invalid argument, n=%8.5f", a);
918 lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);
919 }
920 }
921 else
922 {
923 ret_code=ERR_GSL_FINMIN_PRO_FAILED;
924 if (debug)
925 {
926 sprintf(log_msg," failed, gsl_errno=%d", status);
927 lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);
928 }
|
929 arta 1.1 }
|
930 scholl 1.11 gsl_min_fminimizer_free(s);
|
931 scholl 1.4 return ret_code;
|
932 arta 1.1 }
933
934 if (debug==2)
935 {
|
936 scholl 1.4 sprintf(log_msg," a:%8.5f b:%8.5f m:%8.5f iter=%d, b-a:%f", a,b,m,iter,b-a);
937 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);
|
938 arta 1.1 }
939
940
941 gsl_min_fminimizer_free (s);
942
943 return m;
944
945 }
|
946 scholl 1.11 //---------------------------------------------------------------------------------------------------------------------
947 // Make full ldfs
948 //---------------------------------------------------------------------------------------------------------------------
949 float median(float * tmed, int siz)
950 {
951 int m=siz%2;
952 int s=siz/2;
953 return m==0?(tmed[s-1]+tmed[s])/2:(tmed[s]);
954 }
955
956 int mk_fldfs(float cmx, float cmy, double radius, int naxis_row, int naxis_col, long npixels,
957 float *data, float **save_full_ldf, int *bins1, int *bins2, FILE *opf, int debug)
958 {
959 static char *log_msg_code="mk_fldfs";
960
961 int status=0;
962 int retcode=0;
963 int fulldf_nrows, fulldf_ncols=2;
964 // compute radius array
965 float *data2 = (float *) malloc(sizeof(float) * npixels);
966 if(!data2)
967 scholl 1.11 {
968 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (data2)", log_msg_code, opf);
969 return ERR_MALLOC_FAILED;
970 }
971
972 unsigned long ti,tx,ty;
973 float tdx,tdy;
974 unsigned long cnt=0;
975 for (tx=0;tx<naxis_col;tx++)
976 {
977 for (ty=0;ty<naxis_row;ty++)
978 {
979 ti=naxis_col*ty+tx;
980 if (data[ti]>-2147483648.)
981 {
982 tdx=(float)fabs(tx-cmx);
983 tdy=(float)fabs(ty-cmy);
984 data2[ti]=(float)sqrt(tdx*tdx + tdy*tdy);
985 cnt++;
986 } else data2[ti]=900000.;
987 }
988 scholl 1.11 }
989 if (debug)
990 {
991 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "building array", log_msg_code, opf);
992 }
993 // index them
994 unsigned long *indx;
995 indx=lvector(1,npixels,&status);
996 if (status<0)
997 {
998 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "lvector(indx)", opf);
999 return ERR_MALLOC_FAILED;
1000 }
1001 if (debug)
1002 {
1003 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "vector", log_msg_code, opf);
1004 } retcode=indexx(npixels,data2,indx);
1005 if (retcode<0)
1006 {
1007 lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "indexx", opf);
1008 return ERR_NR_STACK_TOO_SMALL;
1009 scholl 1.11 }
1010 if (debug)
1011 {
1012 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "indexx", log_msg_code, opf);
1013 }
1014 // make bins
1015 int rc=4; // size in pixels of DeltaR of the last bin before the limb: 2*the distorsion
1016 float v1=(float)(1.-(rc/radius));
1017 float v2=1/(1-(v1*v1));
1018 int bins=(int)v2-1;
1019 unsigned int st=(unsigned int)(M_PI*radius*radius);
1020 int sb=st/bins;
1021 int ns=round(sb);
1022 int sr=st%bins;
1023 if (sr==0) *bins1=bins; else *bins1=bins+1;
1024 int mbins=bins+5; // to get what is over the limb: +1 = residual, then +4 outside
1025 fulldf_nrows=mbins;
1026 *bins2=mbins;
1027 // compute them
1028 float *ttmed1, *ttmed2;
1029 // printf("pixels %ld cnt: %lu BINS: %d, MBINS: %d rs: %f st=%u, sb=%d, sr=%d, ns=%d \n",npixels,cnt,bins,mbins,radius,st,sb,sr,ns);
1030 scholl 1.11
1031 ttmed1=vector(1,ns,&status);
1032 if (status<0)
1033 {
1034 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "vector(ttmed1)", opf);
1035 return ERR_MALLOC_FAILED;
1036 }
1037 ttmed2=vector(1,ns,&status);
1038 if (status<0)
1039 {
1040 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "vector(ttmed2)", opf);
1041 return ERR_MALLOC_FAILED;
1042 }
1043 float *t_med_int = (float *) malloc((mbins)*sizeof(float));
1044 if(!t_med_int)
1045 {
1046 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (t_med_int)", log_msg_code, opf);
1047 return ERR_MALLOC_FAILED;
1048 }
1049 float *t_med = (float *) malloc((mbins)*sizeof(float));
1050 if(!t_med)
1051 scholl 1.11 {
1052 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (t_med)", log_msg_code, opf);
1053 return ERR_MALLOC_FAILED;
1054 }
1055 int tk,next_i=0;
1056 for(tk=0;tk<mbins;tk++)
1057 {
1058 if (tk<bins || tk>bins)
1059 {
1060 for(ti=0;ti<ns;ti++)
1061 {
1062 ttmed1[ti]=data[indx[next_i+ti]];
1063 ttmed2[ti]=data2[indx[next_i+ti]];
1064 }
1065 retcode=sort(ns,ttmed1);
1066 if (retcode<0)
1067 {
1068 lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(1)", opf);
1069 return ERR_NR_STACK_TOO_SMALL;
1070 }
1071 /* if (debug)
1072 scholl 1.11 {
1073 lf_logmsg("DEBUG", "APP", NULL, status, "sort 1a", log_msg_code, opf);
1074 } retcode=sort(ns,ttmed2);
1075 */
1076 retcode=sort(ns,ttmed2);
1077 if (retcode<0)
1078 {
1079 lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(2)", opf);
1080 return ERR_NR_STACK_TOO_SMALL;
1081 }
1082 /* if (debug)
1083 {
1084 lf_logmsg("DEBUG", "APP", NULL, status, "sort 2a", log_msg_code, opf);
1085 }
1086 */
1087 t_med_int[tk]=median(ttmed1,ns);
1088 t_med[tk]=median(ttmed2,ns);
1089 next_i=(tk+1)*ns;
1090 }
1091 else if (tk == bins && sr != 0)
1092 {
1093 scholl 1.11 for(ti=0;ti<sr;ti++)
1094 {
1095 ttmed1[ti]=data[indx[next_i+ti]];
1096 ttmed2[ti]=data2[indx[next_i+ti]];
1097 }
1098 retcode=sort(sr,ttmed1);
1099 if (retcode<0)
1100 {
1101 lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(3)", opf);
1102 return ERR_NR_STACK_TOO_SMALL;
1103 }
|
1104 scholl 1.13 /* if (results->debug)
|
1105 scholl 1.11 {
|
1106 scholl 1.13 lf_logmsg("DEBUG", "APP", NULL, status, "sort 1b", log_msg_code, results->opf);
|
1107 scholl 1.11 }
1108 */
1109 retcode=sort(sr,ttmed2);
1110 if (retcode<0)
1111 {
1112 lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(4)", opf);
1113 return ERR_NR_STACK_TOO_SMALL;
1114 }
|
1115 scholl 1.13 /* if (results->debug)
|
1116 scholl 1.11 {
|
1117 scholl 1.13 lf_logmsg("DEBUG", "APP", NULL, status, "sort 2b", log_msg_code, results->opf);
|
1118 scholl 1.11 }
1119 */
1120 t_med_int[tk]=median(ttmed1,sr);
1121 t_med[tk]=median(ttmed2,sr);
1122 next_i=next_i+sr;
1123 }
1124 }
1125 if (debug)
1126 {
1127 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "after all sorts", log_msg_code, opf);
1128 }
1129 // combine the 2 arrays = > pas necessaire pour l'integration dans le vrai fits
1130 *save_full_ldf = (float *) malloc((fulldf_nrows*fulldf_ncols)*sizeof(float));
1131 if(!save_full_ldf)
1132 {
1133 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_full_ldf)", log_msg_code, opf);
1134 return ERR_MALLOC_FAILED;
1135 }
1136
1137
1138 float *p_full1=&t_med_int[0];
1139 scholl 1.11 float *pl_full1=&t_med_int[fulldf_nrows-1];
1140 float *p_full2=&t_med[0];
1141 float *pl_full2=&t_med[fulldf_nrows-1];
1142 float *p_save_full_ldf=save_full_ldf[0]; //&save_full_ldf2[0];
1143 while (p_full1 <= pl_full1) *(p_save_full_ldf++)=*(p_full1++);
1144 while (p_full2 <= pl_full2) *(p_save_full_ldf++)=*(p_full2++);
1145
1146 // free
1147 free(data2);
1148 free(t_med_int);
1149 free(t_med);
1150 // plus all others... vectors...
1151 free_lvector(indx,1,npixels);
1152 free_vector(ttmed1,1,ns);
1153 free_vector(ttmed2,1,ns);
1154 if (debug)
1155 {
1156 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "end flds", log_msg_code, opf);
1157 }
1158 return 0;
1159 }
|
Return to
limbfit.c
CVS log
Up to [Development] / JSOC / proj / limbfit / apps