JSOC_Documentation: proj/limbfit/apps/limbfit.c Source File

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