JSOC_Documentation: proj/limbfit/apps/limbfit

00001 /* 
00002     I.Scholl 
00003 
00004         limbfit(): adapted from Solar Astrometry Program (solarSDO.c)
00005                                     Marcelo Emilio - Jan 2010 
00006 
00007 
00008     #define CODE_NAME       "limbfit_tas"
00009     #define CODE_VERSION    "V6.01" 
00010     #define CODE_DATE       "Fri Mar  6 22:30:42 PST 2015" 
00011 */
00012 
00013 #include "limbfit_tas.h"
00014 
00015 #define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))
00016 
00017 void sav_b0(float *pf_sb0, float *pl_sb0, float *pf_b0)
00018 {
00019     float *p_sb0=pf_sb0;
00020     float *p_b0=pf_b0;
00021     while(p_sb0<=pl_sb0) *(p_sb0++)=*(p_b0++);
00022 }
00023 
00024 void sum_b0(float *beta, float *pf_b0, float *pl_b0)
00025 {
00026     float *p_b0=pf_b0;
00027     float *p_beta=beta;
00028     while(p_b0 <= pl_b0) *(p_b0++)=*(p_b0)+*(p_beta++);             
00029 }
00030 
00031 int limbfit(LIMBFIT_INPUT *input, LIMBFIT_OUTPUT *results, LIMBFIT_IO_PUT *ios)
00032 {
00033 static char *log_msg_code="limbfit";
00034 char log_msg[200];
00035 
00036 if (results->debug) lf_logmsg("DEBUG", "APP", 0, 0, "", log_msg_code, results->opf);
00037 
00038 /************************************************************************
00039                     INIT VARS coming from build_lev1.c
00040 ************************************************************************/   
00041     float *data=input->data;
00042 //  float MIN_VALUE=BAD_PIXEL_VALUE;
00043 //  float MAX_VALUE=(float)fabs(BAD_PIXEL_VALUE);
00044 
00045 
00046 int     ret_code = 0;
00047 
00048 int     w        = ANNULUS_WIDTH;
00049 int     S        = MAX_SIZE_ANN_VARS;
00050 int     nang     = NUM_LDF; //180
00051 int     nprf     = NUM_RADIAL_BINS; //64
00052 int     nreg     = NUM_AB_BINS;
00053 float   rsi      = LO_LIMIT;
00054 float   rso      = UP_LIMIT;
00055 int     r_size   = NUM_FITPNTS;
00056 int     grange   = GUESS_RANGE;
00057 float   dx       = INC_X;
00058 float   dy       = INC_Y;
00059 int     naxis_row= input->img_sz0;
00060 int     naxis_col= input->img_sz1;
00061 float   lahi     = AHI;
00062 int     fldf;
00063 //int       skip     = SKIPGC;
00064 //float     flag     = BAD_PIXEL_VALUE; 
00065 int     iter     = NB_ITER;
00066 if (input->iter!=NB_ITER) iter=input->iter; 
00067 fldf=input->fldf; 
00068 
00069 if (input->spe==1)
00070 {
00071     iter=NB_ITER2;
00072     lahi=AHI2;
00073     rsi=LO_LIMIT2;
00074     rso=UP_LIMIT2;
00075 }
00076     
00077 /************************************************************************/
00078 /*                        set parameters                               */
00079 /************************************************************************/
00080 long npixels=naxis_row*naxis_col;
00081 long ii, jj, i,j;
00082 int jk; 
00083 float cmx, cmy,r;//, guess_cx, guess_cy, guess_r;
00084 int nitr=0, ncut=0;
00085 
00086 r  = (float)naxis_row/2;
00087     
00088 /* Initial guess estimate of center position from Richard's code*/
00089     cmx=(float)input->ix;
00090     cmy=(float)input->iy;
00091     r=(float)input->ir;
00092 
00093 
00094 /************************************************************************/
00095 /*                        Make annulus data                             */
00096 /************************************************************************/
00097     int nbc=3;
00098     float *p_anls=ios->pf_anls;     
00099 
00100     if (ios->is_firstobs == 0) 
00101     {
00102         ios->is_firstobs=1;
00103         float d;
00104         float w2p=r+(float)w/2;
00105         float w2m=r-(float)w/2;
00106         double iimcmy2,jjmcmx;
00107         /* Select Points */
00108         jk=0;//-1;
00109         
00110         for(ii = 0; ii < naxis_row; ii++)
00111         {
00112             iimcmy2=(ii-cmy)*(ii-cmy);
00113             for(jj = 0; jj < naxis_col; jj++) 
00114             { 
00115                 jjmcmx=jj-cmx;
00116                 d=(float)sqrt(iimcmy2+(jjmcmx)*(jjmcmx));
00117                 if (d<=w2p && d>=w2m) 
00118                 {
00119                     jk++;
00120                     *(p_anls++)=(float)jj;
00121                     *(p_anls++)=(float)ii;
00122                     *(p_anls++)=data[ii*naxis_col+jj];                   
00123                 }
00124             }
00125         }
00126         ios->anls_nbpix=jk;
00127         ios->pl_anls=p_anls-1;
00128     
00129         if (results->debug)
00130         {
00131             sprintf(log_msg," jk = %d", jk);
00132             lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00133         }
00134         if ((jk*3) >= S) 
00135         {
00136             lf_logmsg("ERROR", "APP", ERR_SIZE_ANN_TOO_BIG, 0,"nbc>S", log_msg_code, results->opf);
00137             return ERR_SIZE_ANN_TOO_BIG;
00138         }
00139     }
00140     else 
00141     {
00142         jk=ios->anls_nbpix;
00143         ii=0;
00144         p_anls=p_anls+2;
00145         while(p_anls<=ios->pl_anls) 
00146         {
00147             *(p_anls)=data[(int)ios->anls[nbc*ii+1]*naxis_col+(int)ios->anls[nbc*ii]];
00148             ii++;
00149             p_anls=p_anls+3;
00150         }       
00151     }
00152 
00153 /************************************************************************/
00154 /*                  Call Fortran Code Subrotine limb.f                  */
00155 /************************************************************************/
00156     float *rprf, *lprf, *alph, *beta, *b0; //, *sb0, *beta1, *beta2, *beta3;
00157     long ab_nrow=nreg, ab_ncol=2;
00158 
00159     // contains the LDFs AXIS
00160     rprf = (float *) malloc(sizeof(float)*(nprf));
00161         if(!rprf) 
00162         {
00163             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (rprf)", log_msg_code, results->opf);
00164             return ERR_MALLOC_FAILED;
00165         }
00166     float *p_rprf, *pl_rprf;
00167     // contains LDFs inc. average ldf
00168     lprf = (float *) malloc(sizeof(float)*((nang+1)*nprf));
00169         if(!lprf) 
00170         {
00171             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (lprf)", log_msg_code, results->opf);
00172             return ERR_MALLOC_FAILED;
00173         }
00174     float *p_lprf, *pl_lprf;
00175     
00176     alph = (float *) malloc(sizeof(float)*(nreg));
00177         if(!alph) 
00178         {
00179             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (alph)", log_msg_code, results->opf);
00180             return ERR_MALLOC_FAILED;
00181         }
00182     float *pf_alph=&alph[0];
00183     float *p_alph;
00184     float *pl_alph=&alph[ab_nrow-1];
00185     
00186     beta = (float *) malloc(sizeof(float)*(nreg));
00187         if(!beta) 
00188         {
00189             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (xbeta)", log_msg_code, results->opf);
00190             return ERR_MALLOC_FAILED;
00191         }
00192     float *pf_beta=&beta[0];
00193     float *p_beta;
00194     
00195     b0 = (float *) malloc(sizeof(float)*(nreg));
00196         if(!b0) 
00197         {
00198             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (b0)", log_msg_code, results->opf);
00199             return ERR_MALLOC_FAILED;
00200         }
00201     float *pf_b0=&b0[0];
00202     float *p_b0;
00203     float *pl_b0=&b0[ab_nrow-1]; 
00204     
00205     /*
00206     sb0 = (float *) malloc(sizeof(float)*(nreg));
00207         if(!sb0) 
00208         {
00209             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (sb0)", log_msg_code, results->opf);
00210             return ERR_MALLOC_FAILED;
00211         }
00212     float *pf_sb0=&sb0[0];
00213     float *p_sb0;
00214     float *pl_sb0=&sb0[ab_nrow-1]; 
00215     */
00216     
00217     p_b0=pf_b0; 
00218     while(p_b0<=pl_b0) *(p_b0++)=0.;
00219 
00220     // fortran call:
00221     if (results->debug) lf_logmsg("DEBUG", "APP", 0, 0, "entering limb", log_msg_code, results->opf);
00222 
00223     int centyp=0; //=0 do center calculation, =1 skip center calculation
00224     //#1
00225     if (input->cc == 0) centyp=1;
00226     
00227     //in limb_ call: beta contains the output beta, b0 is initialized with 0.0
00228     int it=1, ifail=0;  
00229     
00230     do
00231     {
00232         // init everything that is passed to fortran
00233         //
00234         p_alph=pf_alph;
00235         p_beta=pf_beta;
00236         while(p_alph<=pl_alph) 
00237         {
00238             *(p_beta++)=0.;
00239             *(p_alph++)=0.;
00240         }
00241         p_rprf=&rprf[0];
00242         pl_rprf=&rprf[nprf-1];
00243         while(p_rprf<=pl_rprf) *(p_rprf++)=0.;  
00244         p_lprf=&lprf[0];
00245         pl_lprf=&lprf[((nang+1)*nprf)-1];
00246         while(p_lprf<=pl_lprf) *(p_lprf++)=0.;
00247 
00248             if (results->debug)
00249             {
00250                 sprintf(log_msg,"entering limb# %d /%d iters to do)", it,iter);
00251                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00252             }
00253         ifail=0;
00254 /*
00255     printf("before: %d,%f,%f,%f,%f,%f,\n",centyp,lahi,rsi,rso,dx,dy);
00256     printf("--:%d,%d,%d,%f,%f,%f,\n",ifail,ncut,nitr,cmx,cmy,r);
00257     printf("--:%d,%d\n",S,jk);
00258     printf(">>: %f,%f,%f,%f\n",rprf[0],rprf[63],lprf[0],lprf[11583]);
00259     printf(">00>: %f,%f,%f\n",ios->anls[0],ios->anls[1],ios->anls[2]);
00260     printf(">lp>: %f,%f,%f\n",ios->anls[7289973],ios->anls[7289974],ios->anls[7289975]);
00261     printf(">lv>: %f\n",*(ios->pl_anls));
00262     printf(">-2>: %f,%f,%f\n",ios->anls[7289967],ios->anls[7289968],ios->anls[7289969]);
00263     printf(">-1>: %f,%f,%f\n",ios->anls[7289970],ios->anls[7289971],ios->anls[7289972]);
00264     printf(">+1>: %f,%f,%f\n",ios->anls[7289976],ios->anls[7289977],ios->anls[7289978]);
00265     printf(">>: %f,%f,%f,%f,%f,%f\n",alph[0],alph[255],beta[0],beta[255],b0[0],b0[255]);
00266 */
00267         limb_(&ios->anls[0],&jk, &cmx, &cmy, &r, &nitr, &ncut, &rprf[0], &lprf[0],  &rsi, &rso, 
00268                 &dx, &dy, &alph[0], &beta[0], &ifail, &b0[0], &centyp, &lahi); 
00269             if (results->debug)
00270             {
00271                 sprintf(log_msg,"exiting limb ifail= %d", ifail);
00272                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00273             }
00274 /*
00275     printf("after: %d,%f,%f,%f,%f,%f,\n",centyp,lahi,rsi,rso,dx,dy);
00276     printf("--:%d,%d,%d,%f,%f,%f,\n",ifail,ncut,nitr,cmx,cmy,r);
00277     printf("--:%d,%d\n",S,jk);
00278     printf(">>: %f,%f,%f,%f\n",rprf[0],rprf[63],lprf[0],lprf[11583]);
00279     printf(">00>: %f,%f,%f\n",ios->anls[0],ios->anls[1],ios->anls[2]);
00280     printf(">lp>: %f,%f,%f\n",ios->anls[7289973],ios->anls[7289974],ios->anls[7289975]);
00281     printf(">lv>: %f\n",*(ios->pl_anls));
00282     printf(">-2>: %f,%f,%f\n",ios->anls[7289967],ios->anls[7289968],ios->anls[7289969]);
00283     printf(">-1>: %f,%f,%f\n",ios->anls[7289970],ios->anls[7289971],ios->anls[7289972]);
00284     printf(">+1>: %f,%f,%f\n",ios->anls[7289976],ios->anls[7289977],ios->anls[7289978]);
00285     printf(">>: %f,%f,%f,%f,%f,%f\n",alph[0],alph[255],beta[0],beta[255],b0[0],b0[255]);
00286 */
00287         if(ifail==0)    
00288         {
00289             centyp=1;
00290             //if(it==iter && it>1) sav_b0(pf_sb0,pl_sb0,pf_b0);
00291             sum_b0(&beta[0],pf_b0,pl_b0);
00292         }
00293         it++;
00294     }
00295     while(it<=iter && ifail==0);
00296 
00297 
00298 
00299     if (results->debug)
00300     {
00301         sprintf(log_msg," nitr = %6d", nitr);
00302         lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00303         sprintf(log_msg," cmx = %6.2f, cmy = %6.2f", cmx, cmy);
00304         lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00305     }
00306     
00307     if(ifail == 0) // || ifail == 7) 
00308     {       
00309         /************************************************************************/
00310         /*          Compute the mean of the center of the image                 */
00311         /************************************************************************/
00312         double cmean;
00313         if (input->src <2)
00314         {
00315             double ctot=0.0;
00316             long nbp=0;
00317             float limx_m=cmx-250;
00318             float limx_p=cmx+250;
00319             float limy_m=cmy-250;
00320             float limy_p=cmy+250;
00321             for (i=limx_m;i<limx_p;i++)
00322             {
00323                 for (j=limy_m;j<limy_p;j++)
00324                 {
00325                     ctot=ctot+data[i*naxis_col+j]; 
00326                     nbp++;
00327                 }
00328             }
00329             cmean=ctot/nbp;
00330             if (results->debug)
00331             {
00332                 sprintf(log_msg," cmean = %6.4f (ctot= %6.4f , nbp=%ld)", cmean,ctot,nbp);
00333                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00334             }
00335         }
00336         /************************************************************************/
00337         /*                  Compute the Inflection Point                      */
00338         /************************************************************************/
00339         float *LDF, *D, *t_ip, *p_t_ip;
00340         LDF = (float *) malloc(sizeof(float)*(nprf));
00341             if(!LDF) 
00342             {
00343                 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (LDF)", log_msg_code, results->opf);
00344                 return ERR_MALLOC_FAILED;
00345             }
00346 
00347         D = (float *) malloc(sizeof(float)*(nprf));
00348             if(!D)  
00349             {
00350                 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (D)", log_msg_code, results->opf);
00351                 return ERR_MALLOC_FAILED;
00352             }
00353         t_ip    = (float *) malloc((nang)*sizeof(float));
00354             if(!t_ip) 
00355             {
00356                 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_ip)", log_msg_code, results->opf);
00357                 return ERR_MALLOC_FAILED;
00358             }
00359         p_t_ip=&t_ip[0];
00360 
00361         double maxim; 
00362         double radius = 0.;
00363         int cont, c, ret_gsl;
00364         float h;
00365             
00366         //if (iter > 1) ab_ncol=3; else ab_ncol=2;
00367         int degf=6;
00368         //if (iter > 2) degf=4; else degf=6; // problem with degf = 4...
00369         double A[degf];
00370         double erro[degf];
00371 
00372 
00373          /* copy the last pixel of the mean ldf in the prev one to correct a fortran issue */
00374         lprf[(nprf*nang)+63]=lprf[(nprf*nang)+62];
00375 
00376         /* -----------------------------------------------------------------------------
00377         
00378             3 radius:
00379                 * A[1]          = maximum of the Gaussian
00380                 * radius        = maximum of the Gaussian + quadratic fit -> R_LFS
00381                 * ip, t_ip[*]   = inflection point 
00382         
00383         ------------------------------------------------------------------------------- */          
00384             
00385         for (cont=0; cont<=nang; cont++)
00386         { 
00387          //cont=nang;
00388             ret_gsl=0;
00389              /* dx */
00390              h=(float)rprf[1]-rprf[0];
00391             
00392              /* Take the last (mean) LDF from lprf vector */
00393              for(ii = 0; ii < nprf ; ii++) LDF[ii]=lprf[(nprf*cont)+ii];
00394 
00395              /* Calculate the Derivative */
00396              D[0]=(-3*LDF[4]+16*LDF[3]-36*LDF[2]+48*LDF[1]-25*LDF[0])/(12*h);
00397              D[1]=(LDF[4]-6*LDF[3]+18*LDF[2]-10*LDF[1]-3*LDF[0])/(12*h);
00398              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);
00399              D[nprf-2]=(-LDF[nprf-5]+6*LDF[nprf-4]-18*LDF[nprf-3]+10*LDF[nprf-2]+3*LDF[nprf-1])/(12*h);
00400              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);
00401 
00402              /* square the result */
00403             // for(ii = 0; ii < nprf ; ii++) 
00404             //   D[ii]=D[ii]*D[ii]; //pointers here!
00405             float* p_D=&D[0];
00406             float* pl_D=&D[nprf-1];
00407             while(p_D<=pl_D) *(p_D++)=(*(p_D) * *(p_D));
00408 
00409              /* find the maximum */
00410              jj=0;
00411              maxim=-1;
00412              for(ii = 1; ii < nprf-1 ; ii++) 
00413              {
00414                 if (D[ii] > maxim) 
00415                 {
00416                     maxim=(double)D[ii];
00417                     jj=ii;
00418                 }
00419              }
00420             
00421             /* improve the maximum estimation looking at the 2 neibors points */
00422             t_ip[cont]=rprf[jj]-(float)h/2*(D[jj-1]-D[jj+1])/(2*D[jj]-D[jj-1]-D[jj+1]);
00423             if (results->debug)
00424             {
00425                 if (cont==nang)
00426                 {
00427                     sprintf(log_msg," Inflection Point 1: %ld %d %8.5f %8.5f", jj, nprf, rprf[jj], t_ip[cont]);
00428                     lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00429                 }
00430             }
00431         } // endfor-cont
00432         cont--;
00433         // now only the mean ldf
00434         /* FIT A GAUSSIAN PLUS QUADRATIC FUNCTION */
00435 
00436         /* Select Region */
00437         int m_in, m_ex, N;
00438          /* Gaussian + quadratic function */
00439          /* After Launch verify if this number is OK */
00440          /* Should contain the pixels around the 3*FWHM level */
00441         m_in=jj-r_size;
00442         m_ex=jj+r_size;
00443         
00444         if (m_in < 0) m_in = 0;
00445         if (m_ex > nprf-1) m_ex =nprf-1;
00446         N=m_ex-m_in+1;
00447 
00448         /* parameters:- initial guess */
00449         A[0]= maxim;
00450         A[1]= t_ip[cont];
00451         A[2]= 0.5;
00452         A[3]= 12*maxim;
00453         for (c=0;c<degf;c++) erro[c]=0.;
00454     
00455         if (N >= degf)  // degree of freedom 6 parameters, 6 values minimum
00456         {
00457             double px[N], der[N] , sigma[N];
00458             
00459             jj=-1;
00460             for(ii = m_in; ii <= m_ex ; ii++)
00461             {
00462                    jj++;
00463                    px[jj]=rprf[ii];
00464                    der[jj]=D[ii];
00465                    sigma[jj] = 1.;
00466             }
00467             //A[4]= der[N-1]-der[0];
00468             //if (debug==2) fprintf(opf,"%s DEBUG_INFO in limbfit: #: %d\n",LOGMSG1,cont);
00469 
00470             ret_gsl=gaussfit(der, px, sigma, A, erro, N, degf,results->debug,results->opf);
00471 
00472             if (ret_gsl < 0)
00473             {
00474                 if (cont==nang)
00475                 {
00476                     sprintf(log_msg," gaussfit failed for the averaged LDF %d err:%d", cont,ret_gsl);
00477                     lf_logmsg("ERROR", "APP", 0, ret_gsl, log_msg, log_msg_code, results->opf);
00478                 }
00479                 for (c=0;c<degf;c++) 
00480                 {   
00481                     A[c]=0.;
00482                     erro[c]=0.;
00483                 }
00484                 radius=0.;
00485             }
00486             else
00487             {           
00488                 /* FIND THE MAXIMUM OF THE GAUSSIAN PLUS QUADRATIC FUNCTION */
00489                 radius = A[1]; /* Initial Guess */
00490                 radius = fin_min(A, radius, grange, results->debug,results->opf);
00491                 if (radius < 0)
00492                 {
00493                     ret_gsl=(int)radius;
00494                     if (cont==nang)
00495                     {
00496                         sprintf(log_msg," fin_min failed for the averaged LDF %d err:%d", cont, ret_gsl);
00497                         lf_logmsg("ERROR", "APP", 0, ret_gsl, log_msg, log_msg_code, results->opf);
00498                     }
00499                     for (c=0;c<degf;c++) 
00500                     {   
00501                         A[c]=0.;
00502                         erro[c]=0.;
00503                     }
00504                     radius=0.;
00505                 }
00506             }
00507         } 
00508         else 
00509         { 
00510             sprintf(log_msg," Inflection point too close to annulus data border %d", cont);
00511             lf_logmsg("WARNING", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00512             for (c=0;c<degf;c++) erro[c]=0.;
00513             // in this case: A[1]=radius=t_ip[last] 
00514             radius=t_ip[cont];
00515         }
00516         // save them
00517         for (c=0;c<degf;c++)
00518         {
00519             results->fits_as[c]=(float)A[c];
00520             results->fits_es[c]=(float)erro[c];
00521         }
00522         if (results->debug)
00523         {   
00524             sprintf(log_msg," Inflection Point 2: %8.5f %8.5f %8.5f", A[1], erro[1], radius);
00525             lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00526             if (degf==4)
00527             {
00528                 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f ", A[0],A[1],A[2],A[3]);
00529                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00530                 sprintf(log_msg," -----: %8.5f %8.5f %8.5f %8.5f ", erro[0],erro[1],erro[2],erro[3]);
00531                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00532             }
00533             else 
00534             {
00535                 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]);
00536                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00537                 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]);
00538                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00539             }
00540         }
00541     
00542         //  Full ldfs                     
00543         int fldfr=0;    //proc result: 0=ok; 1=failed; 2=cannot be processed; 3=malloc pb; 4=not processed; 
00544         int fulldf_nrows, fulldf_ncols=2;
00545         int bins1=0, bins2=0;
00546         int retcode=0;
00547         float *save_full_ldf;
00548         if (fldf==1)
00549         {
00550             // full LDF
00551 
00552             //if (ret_gsl<0 || cmx < 0.01 || cmy < 0.01 || radius < 0.01) // that I don't remember what for 
00553             if (cmx < 0.01 || cmy < 0.01 || radius < 0.01)
00554             {
00555                 fulldf_nrows=1;
00556                 bins1=0;
00557                 save_full_ldf  = (float *) malloc((2)*sizeof(float));
00558                     if(!save_full_ldf) 
00559                     {
00560                         lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_full_ldf)", log_msg_code, results->opf);
00561                         return ERR_MALLOC_FAILED;
00562                     }               
00563                 save_full_ldf[0]=0;
00564                 save_full_ldf[1]=0;
00565                 fldfr=3;
00566             }
00567             else
00568             {
00569                 retcode=mk_fldfs(cmx, cmy, radius, naxis_row, naxis_col, npixels, data, &save_full_ldf, &bins1, &bins2, results->opf, results->debug);
00570                 if (retcode == ERR_MALLOC_FAILED)
00571                 {
00572                     fldfr=2;
00573                     return ERR_MALLOC_FAILED;
00574                 }
00575                 else
00576                     if (retcode == ERR_NR_STACK_TOO_SMALL)
00577                     {
00578                         ret_code=ERR_LIMBFIT_FLDF_FAILED;
00579                         fldfr=1;
00580                         //shouldnt exist?
00581                     }
00582                 fulldf_nrows=bins2;
00583             }
00584         }           
00585 
00586 
00587         /* ---------------------Save LDFs & AB data in a FITS file--------------
00588         
00589          LDF = lprf(ldfs) & rprf (axis) come from limbfit.f         
00590         
00591                      to be read in IDL as: FLOAT     = Array[65, 182]
00592 
00593                     NAXIS1          FLOAT               65
00594                     NAXIS2          FLOAT              182
00595         
00596             -----------------
00597             | axis   |  fsn |
00598             -----------------
00599             | lfds   |  ip  |
00600             | ...    |  ... |
00601             | ....   |  ... |
00602             -----------------
00603             | avg ldf|  ip  |
00604             -----------------
00605         
00606         --------------------------------------------------------------------- */
00607         // for saving them in FITS file
00608         float *save_ldf, *save_alpha_beta; 
00609         long ldf_nrow=nang+2, ldf_ncol=nprf+1;  //ii -nbr of ldfs, jj -nbr of points for each ldf
00610         save_ldf    = (float *) malloc((ldf_nrow*ldf_ncol)*sizeof(float));
00611             if(!save_ldf) 
00612             {
00613                 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_ldf)", log_msg_code, results->opf);
00614                 return ERR_MALLOC_FAILED;
00615             }
00616         save_alpha_beta = (float *) malloc((ab_nrow*ab_ncol)*sizeof(float));
00617             if(!save_alpha_beta) 
00618             {
00619                 lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_alpha_beta)", log_msg_code, results->opf);
00620                 return ERR_MALLOC_FAILED;
00621             }       
00622 
00623         float *p_sldf=&save_ldf[0];
00624         //float *pl_sldf=&save_ldf[(ldf_nrow*ldf_ncol)-1];
00625         float *plc_sldf;
00626 
00627         // axis
00628         p_rprf=&rprf[0];
00629         pl_rprf=&rprf[nprf-1];
00630         while (p_rprf <= pl_rprf)
00631             *(p_sldf++)=*(p_rprf++);
00632 
00633 //      *(p_sldf++)=(float)input->fsn;
00634         *(p_sldf++)=0.0;
00635 
00636 /*old       
00637             
00638         //ldfs
00639         p_lprf=&lprf[0];
00640         p_sldf=&save_ldf[ldf_ncol]; //      shouldn't be? pl_lprf=&lprf[ldf_ncol]; ???
00641         while (p_sldf <= pl_sldf)
00642             *(p_sldf++)=*(p_lprf++);
00643 */          
00644         int cpt_row=0;
00645         p_lprf=&lprf[0];
00646         pl_lprf=&lprf[((nang+1)*nprf)-1];
00647 
00648         while (p_lprf <= pl_lprf)
00649         {
00650             plc_sldf=&lprf[((cpt_row+1)*nprf)-1];
00651             while (p_lprf <= plc_sldf)
00652                 *(p_sldf++)=*(p_lprf++);
00653             *(p_sldf++)=*(p_t_ip++);
00654             cpt_row++;
00655         }
00656 
00657         /* ---------------------------------------------------------------------
00658         
00659          AB         write first alpha then beta
00660                     to be read in IDL as:  Array[2, 256]
00661         
00662                     NAXIS1          FLOAT                2
00663                     NAXIS2          FLOAT              256
00664  
00665         --------------------------------------------------------------------- */
00666         p_alph=pf_alph;
00667         p_b0=pf_b0;
00668         float *p_save_alpha_beta=&save_alpha_beta[0];
00669         while (p_alph <= pl_alph) 
00670         {
00671             *(p_save_alpha_beta++)=*(p_alph++);
00672             *(p_save_alpha_beta++)=*(p_b0++);
00673         }
00674         /*
00675         if (iter>1) 
00676         {
00677             p_sb0=pf_sb0;
00678             while (p_sb0 <= pl_sb0) *(p_save_alpha_beta++)=*(p_sb0++);
00679         }
00680         */ 
00681 
00682         // Update Returned Structure when process succeeded                    
00683         results->cenx=cmx;
00684         results->ceny=cmy;
00685         results->radius=radius;
00686         results->cmean=cmean;
00687         results->numext=3;      
00688         results->fits_ldfs_naxis1=ldf_ncol; 
00689         results->fits_ldfs_naxis2=ldf_nrow;
00690         results->fits_fldfs_tfields=fulldf_ncols;
00691         results->fits_fldfs_nrows=fulldf_nrows;
00692         results->fits_ab_naxis1=ab_ncol;
00693         results->fits_ab_naxis2=ab_nrow;
00694         results->nb_fbins=bins1;
00695         results->ann_wd=w;
00696         results->mxszannv=S;
00697         results->nb_ldf=nang;
00698         results->nb_rdb=nprf;
00699         results->nb_abb=nreg;
00700         results->up_limit=rsi;
00701         results->lo_limit=rso;
00702         results->inc_x=dx;
00703         results->inc_y=dy;
00704         results->nfitpnts=r_size;
00705         results->nb_iter=iter;
00706         results->fldfr=fldfr;
00707         results->ahi=lahi;
00708         results->error1=ifail;      
00709         results->error2=ret_gsl;        
00710         if (ifail == 0) 
00711             results->quality=9; 
00712         //? if (cont>=nang && ret_gsl<0)
00713         if (ret_gsl<0)
00714         {
00715             results->quality=1; 
00716             ret_code=ERR_LIMBFIT_FIT_FAILED;
00717             if (results->debug)
00718             {   
00719                 sprintf(log_msg,"  ret_gsl<0 = %2d", ret_gsl);
00720                 lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00721             }       
00722         }
00723         results->fits_ldfs_data=save_ldf; 
00724         results->fits_alpha_beta=save_alpha_beta; 
00725         if (fldf == 1) results->fits_fulldfs=save_full_ldf; else results->fits_fulldfs=0;
00726         free(D);
00727         free(LDF);
00728         free(t_ip);
00729     } // end limb OK
00730     else 
00731     {
00732         if (results->debug)
00733         {   
00734             sprintf(log_msg," limb.f routine returned ifail = %2d", ifail);
00735             lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00736         }       
00737         // Update Returned Structure when process failed                    
00738         results->numext=0;      
00739         results->error1=ifail;
00740         results->error2=-1;
00741         results->quality=0;
00742         ret_code=ERR_LIMBFIT_FAILED;
00743         results->ann_wd=w;
00744         results->mxszannv=S;
00745         results->nb_ldf=nang;
00746         results->nb_rdb=nprf;
00747         results->nb_abb=nreg;
00748         results->up_limit=rsi;
00749         results->lo_limit=rso;
00750         results->inc_x=dx;
00751         results->inc_y=dy;
00752         results->nfitpnts=r_size;
00753         results->nb_iter=iter;
00754         results->ahi=lahi;
00755         results->cmean=0;
00756         results->nb_fbins=0;
00757         results->fldfr=4;
00758     } // end limb failed
00759     // IS: do not free those (save_ldf,save_params,save_alpha_beta,save_full_ldf) passed from or to the structure !
00760 
00761     free(rprf);
00762     free(alph);
00763     free(beta);
00764     free(b0);
00765     //free(sb0);
00766 //  free(lprf);
00767 
00768     if (results->debug)
00769     {   
00770         sprintf(log_msg," >>>>end of limbfit with: %d", ret_code);
00771         lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00772     }
00773     sprintf(log_msg," end: RC: %d - limb:%d - fit:%d - quality: %d", ret_code, results->error1, results->error2, results->quality);
00774     lf_logmsg("INFO", "APP", 0, 0, log_msg, log_msg_code, results->opf);
00775 
00776 return(ret_code);
00777 
00778 } /* end of main */
00779 
00780 
00781 /*--------------------------------------------------------------------------*/
00782 int gaussfit(double y[], double t[],double sigma[], double A[], double erro[], long N, int degf, int debug, FILE *opf)
00783 /* Calculate a Least SqrareGaussian + Quadratic fit              */
00784 /*      Uses the GNU Scientific Library                          */
00785 /* Marcelo Emilio (c) v 1.0 Jan 2009                             */
00786 /* fits A[0]exp(-z^2/2)+A[3]+A[4]*x+A[5]*x^2                     */
00787 /* z=(t-A[1])/A[2]                                               */
00788 /* Need to add :                                                 */
00789 /* #include <gsl/gsl_rng.h>                                      */
00790 /* #include <gsl/gsl_randist.h>                                  */
00791 /* #include <gsl/gsl_vector.h>                                   */
00792 /* #include <gsl/gsl_blas.h>                                     */
00793 /* #include <gsl/gsl_multifit_nlin.h>                            */
00794 /* #include "expfit.c"                                           */
00795 /* compiles as  gcc program.c --lm -lgsl -lgslcblas              */
00796 /* y --> f(x) - ordinate values                                  */
00797 /* t --> x(i) - abscissa values                                  */
00798 /* sigma --> independent gaussian erros                          */
00799 /* N --> Number of points in the vector                          */
00800 /* A --> Input a initial guess and output the  result            */
00801 /* erro --> Chi-square of A                                      */
00802 
00803 {
00804     static char *log_msg_code="gaussfit";
00805     int ret_code=0;
00806     char log_msg[200];
00807 
00808     const gsl_multifit_fdfsolver_type *T;
00809     gsl_multifit_fdfsolver *s;
00810     int status;
00811     unsigned int i, iter = 0;
00812     const size_t n = N;
00813     const size_t p = degf;
00814          
00815     gsl_matrix *covar = gsl_matrix_alloc (p, p);
00816     /* double t[N], y[N], sigma[N]; */
00817     struct dataF d = { n, t, y, sigma};
00818     gsl_multifit_function_fdf f;
00819     
00820     /* Initial Guess */
00821     double x_init[degf];
00822     for (i=0; i <degf; i++) 
00823        x_init[i]=A[i];
00824     
00825     gsl_vector_view x = gsl_vector_view_array (x_init, p);
00826     const gsl_rng_type * type;
00827     gsl_rng * r;
00828          
00829     gsl_rng_env_setup();
00830          
00831     type = gsl_rng_default;
00832     r = gsl_rng_alloc (type);
00833          
00834     f.f = &expb_f;
00835     f.df = &expb_df;
00836     f.fdf = &expb_fdf;
00837     f.n = n;
00838     f.p = p;
00839     f.params = &d;
00840 
00841     T = gsl_multifit_fdfsolver_lmsder;
00842     s = gsl_multifit_fdfsolver_alloc (T, n, p);
00843     gsl_multifit_fdfsolver_set (s, &f, &x.vector);  
00844 
00845     do
00846     {
00847         iter++;
00848         status = gsl_multifit_fdfsolver_iterate (s);     
00849 
00850         if (status==0 || status == GSL_ENOPROG) 
00851             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 */ 
00852         else 
00853         {
00854             ret_code=ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED; 
00855             if (debug) 
00856             {
00857                 sprintf(log_msg," iter: %u, status: %d (%s)\n", iter,status,gsl_strerror (status));
00858                 lf_logmsg("WARNING", "APP", ret_code, status, log_msg, log_msg_code, opf);          
00859             }
00860         }
00861 
00862     }
00863     while (status == GSL_CONTINUE && iter < 500);
00864 
00865             
00866     if (status == 0)  
00867     {
00868 
00869         gsl_multifit_covar (s->J, 0.0, covar);
00870              
00871         #define FIT(i) gsl_vector_get(s->x, i)
00872         #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
00873              
00874         double chi = gsl_blas_dnrm2(s->f);
00875         double dof = (double)(n - p);
00876         double c = GSL_MAX_DBL(1, chi / sqrt(dof)); 
00877           
00878         
00879         for (i=0; i <degf; i++) {
00880             A[i]=FIT(i);
00881             erro[i]=c*ERR(i);
00882         }  
00883     }
00884     else
00885     {
00886         ret_code=ERR_GSL_GAUSSFIT_FDFSOLVER_FAILED;         
00887         if (debug)
00888         {
00889             sprintf(log_msg," (gsl_multifit_fdfsolver_iterate) failed, (iter=%u) gsl_errno=%d", iter,status);
00890             lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);
00891         }
00892     }
00893 
00894     
00895     if (debug==2)
00896     {                                                                       
00897         sprintf(log_msg," Nonlinear LSQ fitting status = %s (%d)", gsl_strerror (status),status);
00898         lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);
00899     }
00900     if (debug) 
00901     {
00902         sprintf(log_msg," finished at iter# %u\n", iter);
00903         lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);            
00904     }
00905     gsl_multifit_fdfsolver_free (s);
00906     gsl_matrix_free (covar);
00907     gsl_rng_free (r);
00908     return(ret_code);
00909 }
00910 
00911 
00912 /*--------------------------------------------------------------------------*/
00913 
00914 double fin_min(double A[], double m, int range, int debug, FILE *opf)
00915 /* Calculate the maximum of the quadratic + gaussian function    */
00916 /*                 using Brent algorithm                         */
00917 /* Marcelo Emilio (c) v 1.0 Jan 2009                             */
00918 /* A are the parameters of the gaussian                          */
00919 /* m is the guess where the function is maximum                  */
00920 /* Calculate the maximum around                                  */
00921 /*                 (-5 + m) < m > (5 - m)  ( in pixel units )    */
00922 /* Returns the m vaule where the function is maximum             */
00923 /* requires: #include "expmax.c" --> definition of the function  */
00924 /*           #include <gsl/gsl_errno.h>                          */
00925 /*           #include <gsl/gsl_math.h>                           */
00926 /*           #include <gsl/gsl_min.h>                            */
00927 {
00928     static char *log_msg_code="fin_min";
00929     int status;
00930     gsl_set_error_handler_off();
00931     int ret_code=0;
00932     char log_msg[200];
00933     
00934      int iter = 0, max_iter = 1000;
00935      const gsl_min_fminimizer_type *T;
00936      gsl_min_fminimizer *s;
00937      
00938      gsl_function F;
00939     
00940     //double m_expected = m+0.01; //1e-8;
00941     //double a = m-5000, b = m+5000;
00942     double a = m-range, b = m+range;      // initially = 2
00943 //  if (degf == 4)
00944     //struct exp_plus_quadratic_params params= { A[0], A[1], A[2], A[3] }; 
00945 //  else
00946     struct exp_plus_quadratic_params params= { A[0], A[1], A[2], A[3], A[4], A[5] }; 
00947       F.function = &exp_plus_quadratic_function;
00948       F.params = &params;
00949          
00950       T = gsl_min_fminimizer_brent;
00951       s = gsl_min_fminimizer_alloc (T);
00952     
00953       status=gsl_min_fminimizer_set (s, &F, m, a, b);     
00954       if (status)
00955       {
00956             if (status == GSL_EINVAL) 
00957             {
00958                 ret_code=ERR_GSL_FINMIN_NOMIN_FAILED;
00959                 if (debug)
00960                 {
00961                     sprintf(log_msg," (gsl_min_fminimizer_set) doesn't find a minimum, m=%f", m);
00962                     lf_logmsg("DEBUG", "APP", 0, status, log_msg, log_msg_code, opf);
00963                 }
00964             } 
00965             else 
00966             {
00967                 ret_code=ERR_GSL_FINMIN_SET_FAILED;
00968                 if (debug)
00969                 {
00970                     sprintf(log_msg," (gsl_min_fminimizer_set) failed, gsl_errno=%d", status);
00971                     lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);            
00972                 }
00973             }
00974             gsl_min_fminimizer_free(s);
00975             return ret_code;
00976       }
00977       
00978       do
00979       {
00980        iter++;
00981        status = gsl_min_fminimizer_iterate (s);
00982          
00983        m = gsl_min_fminimizer_x_minimum (s);
00984        a = gsl_min_fminimizer_x_lower (s);
00985        b = gsl_min_fminimizer_x_upper (s);
00986          
00987        status = gsl_min_test_interval (a, b,1E-3, 0.0);
00988       }
00989       while (status == GSL_CONTINUE && iter < max_iter);
00990       
00991       if (status)       // regarder lequel plantait et verifier si on a autre chose que GSL_EINVAL comme erreur
00992         {               // ajouter debug pour les messages apres
00993             if (status == GSL_EINVAL) 
00994             {
00995                 ret_code=ERR_GSL_FINMIN_PRO_FAILED;
00996                 if (debug)
00997                 {
00998                     sprintf(log_msg," invalid argument, n=%8.5f", a);
00999                     lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);            
01000                 }
01001             } 
01002             else 
01003             {           
01004                 ret_code=ERR_GSL_FINMIN_PRO_FAILED;
01005                 if (debug)
01006                 {
01007                     sprintf(log_msg," failed, gsl_errno=%d", status);
01008                     lf_logmsg("ERROR", "APP", ret_code, status, log_msg, log_msg_code, opf);            
01009                 }
01010             }
01011             gsl_min_fminimizer_free(s);
01012             return ret_code; 
01013         }
01014       
01015       if (debug==2)
01016       {
01017             sprintf(log_msg," a:%8.5f b:%8.5f m:%8.5f iter=%d, b-a:%f", a,b,m,iter,b-a);
01018             lf_logmsg("DEBUG", "APP", 0, 0, log_msg, log_msg_code, opf);
01019       }
01020 
01021 
01022       gsl_min_fminimizer_free (s);
01023          
01024       return m;
01025 
01026 }
01027 //---------------------------------------------------------------------------------------------------------------------
01028 // Make full ldfs
01029 //---------------------------------------------------------------------------------------------------------------------
01030 float median(float * tmed, int siz)
01031 {
01032     int m=siz%2;
01033     int s=siz/2;
01034     return m==0?(tmed[s-1]+tmed[s])/2:(tmed[s]);
01035 }
01036 
01037 int mk_fldfs(float cmx, float cmy, double radius, int naxis_row, int naxis_col, long npixels, 
01038                     float *data, float **save_full_ldf, int *bins1, int *bins2, FILE *opf, int debug)
01039 {
01040 static char *log_msg_code="mk_fldfs";
01041 
01042     int status=0;
01043     int retcode=0;
01044     int fulldf_nrows, fulldf_ncols=2;
01045     // compute  radius array
01046     float *data2 = (float *) malloc(sizeof(float) * npixels);
01047         if(!data2) 
01048         {
01049             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (data2)", log_msg_code, opf);
01050             return ERR_MALLOC_FAILED;
01051         }
01052 
01053     unsigned long ti,tx,ty;
01054     float tdx,tdy;
01055     unsigned long cnt=0;
01056     for (tx=0;tx<naxis_col;tx++)
01057     {
01058         for (ty=0;ty<naxis_row;ty++)
01059         {
01060             ti=naxis_col*ty+tx;
01061             if (data[ti]>-2147483648.)
01062             {   
01063                 tdx=(float)fabs(tx-cmx);
01064                 tdy=(float)fabs(ty-cmy);
01065                 data2[ti]=(float)sqrt(tdx*tdx + tdy*tdy);
01066                 cnt++;
01067             } else data2[ti]=900000.;   
01068         }
01069     }
01070         if (debug) 
01071             {
01072                 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "building array", log_msg_code, opf);          
01073             }
01074     // index them
01075     unsigned long *indx;
01076     indx=lvector(1,npixels,&status); 
01077         if (status<0) 
01078         {
01079             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "lvector(indx)", opf);
01080             return ERR_MALLOC_FAILED;
01081         }
01082         if (debug) 
01083             {
01084                 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "vector", log_msg_code, opf);          
01085             }   retcode=indexx(npixels,data2,indx); 
01086         if (retcode<0) 
01087         {
01088             lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "indexx", opf);
01089             return ERR_NR_STACK_TOO_SMALL;
01090         }
01091         if (debug) 
01092             {
01093                 lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "indexx", log_msg_code, opf);          
01094             }
01095     // make bins            
01096     int rc=4; // size in pixels of DeltaR of the last bin before the limb: 2*the distorsion
01097     float v1=(float)(1.-(rc/radius));
01098     float v2=1/(1-(v1*v1));
01099     int bins=(int)v2-1;     
01100     unsigned int st=(unsigned int)(M_PI*radius*radius);
01101     int sb=st/bins;
01102     int ns=round(sb);
01103     int sr=st%bins;
01104     if (sr==0) *bins1=bins; else *bins1=bins+1;
01105     int mbins=bins+5; // to get what is over the limb: +1 = residual, then +4 outside
01106     fulldf_nrows=mbins;
01107     *bins2=mbins;
01108     // compute them 
01109     float *ttmed1, *ttmed2;
01110     // 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);
01111 
01112     ttmed1=vector(1,ns,&status); 
01113         if (status<0) 
01114         {
01115             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "vector(ttmed1)", opf);
01116             return ERR_MALLOC_FAILED;
01117         }
01118     ttmed2=vector(1,ns,&status); 
01119         if (status<0) 
01120         {
01121             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed", "vector(ttmed2)", opf);
01122             return ERR_MALLOC_FAILED;
01123         }
01124     float *t_med_int    = (float *) malloc((mbins)*sizeof(float));
01125         if(!t_med_int) 
01126         {
01127             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (t_med_int)", log_msg_code, opf);
01128             return ERR_MALLOC_FAILED;
01129         }
01130     float *t_med        = (float *) malloc((mbins)*sizeof(float));
01131         if(!t_med) 
01132         {
01133             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (t_med)", log_msg_code, opf);
01134             return ERR_MALLOC_FAILED;
01135         }
01136     int tk,next_i=0;
01137     for(tk=0;tk<mbins;tk++)
01138     {
01139         if (tk<bins || tk>bins)
01140         {
01141             for(ti=0;ti<ns;ti++) 
01142             {               
01143                 ttmed1[ti]=data[indx[next_i+ti]];
01144                 ttmed2[ti]=data2[indx[next_i+ti]];
01145             }
01146             retcode=sort(ns,ttmed1);
01147                 if (retcode<0) 
01148                 {
01149                     lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(1)", opf);
01150                     return ERR_NR_STACK_TOO_SMALL;
01151                 }
01152 /*              if (debug) 
01153                 {
01154                     lf_logmsg("DEBUG", "APP", NULL, status, "sort 1a", log_msg_code, opf);          
01155                 }           retcode=sort(ns,ttmed2);
01156 */
01157             retcode=sort(ns,ttmed2);
01158                 if (retcode<0) 
01159                 {
01160                     lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(2)", opf);
01161                     return ERR_NR_STACK_TOO_SMALL;
01162                 }
01163 /*              if (debug) 
01164                 {
01165                     lf_logmsg("DEBUG", "APP", NULL, status, "sort 2a", log_msg_code, opf);          
01166                 }
01167 */
01168             t_med_int[tk]=median(ttmed1,ns);
01169             t_med[tk]=median(ttmed2,ns);
01170             next_i=(tk+1)*ns;
01171         } 
01172         else if (tk == bins && sr != 0)
01173         {
01174             for(ti=0;ti<sr;ti++) 
01175             {               
01176                 ttmed1[ti]=data[indx[next_i+ti]];
01177                 ttmed2[ti]=data2[indx[next_i+ti]];
01178             }
01179             retcode=sort(sr,ttmed1);
01180                 if (retcode<0) 
01181                 {
01182                     lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(3)", opf);
01183                     return ERR_NR_STACK_TOO_SMALL;
01184                 }
01185 /*              if (results->debug) 
01186                 {
01187                     lf_logmsg("DEBUG", "APP", NULL, status, "sort 1b", log_msg_code, results->opf);         
01188                 }
01189 */
01190             retcode=sort(sr,ttmed2);
01191                 if (retcode<0) 
01192                 {
01193                     lf_logmsg("ERROR", "APP", ERR_NR_STACK_TOO_SMALL, 0,"stack too small", "sort(4)", opf);
01194                     return ERR_NR_STACK_TOO_SMALL;
01195                 }
01196 /*              if (results->debug) 
01197                 {
01198                     lf_logmsg("DEBUG", "APP", NULL, status, "sort 2b", log_msg_code, results->opf);         
01199                 }
01200 */
01201             t_med_int[tk]=median(ttmed1,sr);
01202             t_med[tk]=median(ttmed2,sr);
01203             next_i=next_i+sr;
01204         }
01205     }
01206     if (debug) 
01207         {
01208             lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "after all sorts", log_msg_code, opf);         
01209         }
01210     // combine the 2 arrays = > pas necessaire pour l'integration dans le vrai fits 
01211     *save_full_ldf  = (float *) malloc((fulldf_nrows*fulldf_ncols)*sizeof(float));
01212         if(!save_full_ldf) 
01213         {
01214             lf_logmsg("ERROR", "APP", ERR_MALLOC_FAILED, 0,"malloc failed (save_full_ldf)", log_msg_code, opf);
01215             return ERR_MALLOC_FAILED;
01216         }
01217 
01218 
01219     float *p_full1=&t_med_int[0];
01220     float *pl_full1=&t_med_int[fulldf_nrows-1];
01221     float *p_full2=&t_med[0];
01222     float *pl_full2=&t_med[fulldf_nrows-1]; 
01223     float *p_save_full_ldf=save_full_ldf[0]; //&save_full_ldf2[0];
01224     while (p_full1 <= pl_full1) *(p_save_full_ldf++)=*(p_full1++);
01225     while (p_full2 <= pl_full2) *(p_save_full_ldf++)=*(p_full2++);  
01226     
01227     // free
01228     free(data2);
01229     free(t_med_int);
01230     free(t_med);
01231     // plus all others...  vectors...
01232     free_lvector(indx,1,npixels);
01233     free_vector(ttmed1,1,ns);
01234     free_vector(ttmed2,1,ns);
01235                 if (debug) 
01236                 {
01237                     lf_logmsg("DEBUG", "APP", DEBUG_MSG, status, "end flds", log_msg_code, opf);            
01238                 }
01239 return 0;
01240 }
proj/limbfit/apps/limbfit_tas.c
