Main Page | Modules | Data Structures | File List | Data Fields | Globals | Related Pages

extract/source/find_mpup.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00015 /* =========================================================== */
00016 
00017 #include <IFU_io.h>
00018 #include <IFU_math.h>
00019 #include <snifs.h>
00020 #include <extract.h>
00021 
00022 /* =========================================================== */
00023 
00025 #define DBG_PREFIX "dbg_fmp"
00026 
00032 #define FITNEIGHBOR 1.6
00033 //#define FITNEIGHBOR -8.
00034 
00038 #define THRESHORDER 0.3
00039 
00041 #define USE_PROFERR 1
00042 
00043 /* Global variables needed by nllsqfit_distMLA and model_distMLA */
00044 int *glmlaml, *glmlanl;        /* Mpup position indices */      
00045 double *glmlaxpup, *glmlaypup; /* Mpup positions        */
00046 double glmlaxlc, glmlaylc;     /* Central lens position */
00047 
00048 /* Global variables describing the geometric PSF */
00049 int glnggeo;               /* Nb of gaussian components */
00050 double *glIgeo, *glXgeo, *glSgeo;
00051 
00052 /* Global variables describing the global PSF */
00053 int glngpsf;               /* Nb of gaussian components */
00054 
00055 /* Global variables describing the core PSF = geom. PSF (x) global PSF */
00056 double **glIpsf, **glSpsf; /* Global variables needed by nllsqfit_globalPSF */
00057 
00058 /* Global variables needed by nllsqfit_globalPSF */
00059 int *glpcen;               /* Length of the central mpup profiles */      
00060 
00061 /* Global variables describing the local PSF */
00062 int glngloc;             /* Nb of gaussian components */
00063 
00064 
00065 /* === Doxygen Comment ======================================= */
00108 /* =========================================================== */
00109 
00110 int main(int argc, char **argv)
00111 {
00112   char **argval, **arglabel;
00113 
00114   TABLE mask,intable,table;
00115   IMAGE2D mpupima;
00116   TIGERfile xcube, ycube;
00117   SPECTRUM *histospec;
00118   SPECTRUM spec;
00119   Channel channel;
00120   
00121   char *pupname, *tabname, *maskname;
00122   char ananame[lg_name+1], tmpname[lg_name+1], dident[lg_ident+1];
00123   int *iadj, *iend, *iwork;
00124   int *ppup, *opup, *ipup, *jpup, *indx, *indx2, *nopup, *lpup;
00125   int *icen, *pcen;
00126   int ANALYSIS, DISCARD_GLOBPSF, TABLE, OFFSET;
00127   int i,j,k,i0,j0,i1,j1,i2,j2, status;
00128   int nbin,npup,nvpup,ntri,ncen,npts,nggeo,ngpsf,ngloc;
00129   int xhalfw,xfullw, yhalfw,yfullw, no;
00130   int colx,coly,colf,coln,colo,colxf,colyf,colixf,coliyf;
00131   int colsxf,colsyf,colrxf,colryf;
00132   int colno, coli, colj, colxin, colyin;
00133   float *xpup, *ypup, *xthe, *ythe, *dpup, *spup;
00134   float *i1g,*x1g,*s1g;
00135   float mindist, maxdist, fracthresh, lbdaref;
00136   float x0,y0;
00137   float ftmp, scale, fmax, lens_px;
00138   double xgau[3], ygau[3], dxgau[3], dygau[3], fl[3],dfl[3];
00139   double fracdist[2], offset[2];
00140   double *xloc, *dxloc;
00141   double *fpup, *rpup;
00142   double *xarr, *xfit, *yarr, *yfit, *xsig;
00143   double *xtri, *ytri, *stri, *ltri, *atri, *dwork;
00144   double *scen;
00145   double lenssize, angle, distor;
00146   double mean,sigma,median,sigmed;
00147   double min,max,threshold;
00148   double rmsx, rmsy, chi2, gof, norm, geonorm, RoN2;
00149   /* Distortion parameters */
00150 #ifdef USE_R4_DIST
00151   int ndpar = 6; /* lenssize,angle,xc,yc,alpha,beta */
00152 #else
00153   int ndpar = 5; /* lenssize,angle,xc,yc,alpha */
00154 #endif
00155   double dpar[ndpar], edpar[ndpar];
00156   /* Global PSF parameters */
00157   int ncpsfpar;
00158   double *cpsfpar, *ecpsfpar, *lcpsfpar, *ucpsfpar;
00159   double *psf_param;
00160   
00161   /* ##### INTRODUCTION ############################## */
00162 
00163   set_purpose("Analysis of micro-pupil frame");
00164   set_arglist("-in none -nbin 100 -threshold 0.99 "
00165               "-fracdist 0.1,0.1 -analysis NULL -mask mask "
00166               "-nggeo 3 -ngpsf 3 -ngloc 1 -lambda none "
00167               "-table NULL -offset NULL");
00168 
00169   init_snifs("$Name:  $","$Revision: 1.30 $");
00170   init_session(argv,argc,&arglabel,&argval);
00171 
00172   if (DEBUG) {
00173     print_msg("$Id: find_mpup.c,v 1.30 2005/09/15 21:43:14 ycopin Exp $");
00174     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00175   }
00176 
00177   /* ===== Parameters ============================== */
00178 
00179   /* ----- Input ------------------------------ */
00180 
00181   pupname = argval[0];          /* mpup frame                                */
00182   get_argval(1,"%d",&nbin);     /* nbin in mpup frame histogram              */
00183   get_argval(2,"%f",&fracthresh); /* frac of pixel population to be studied  */
00184                         /* 1/2-size (x & y) of the analysis area around mpup */
00185   if (sscanf(argval[3], "%lf,%lf",&fracdist[0],&fracdist[1]) != 2) { 
00186     print_error("Invalid option -fracdist fx,fy (got '%s')", argval[3]);
00187     exit_session(ERR_BAD_PARAM);
00188   }
00189   if ((ANALYSIS = is_set(argval[4]))) 
00190     get_argval(4,"%s",ananame);
00191   else if (DEBUG) 
00192     sprintf(ananame,"%s_mpups",DBG_PREFIX);
00193   maskname = argval[5];
00194   get_argval(6,"%d",&nggeo);    /* Gaussian comp. in the geom. PSF  */
00195   get_argval(7,"%d",&ngpsf);    /* Gaussian comp. in the global PSF */
00196   get_argval(8,"%d",&ngloc);    /* Gaussian comp. in the local PSF  */
00197   get_argval(9,"%f",&lbdaref);  /* Lambda ref for the mpup          */
00198   if ((TABLE = is_set(argval[10]))) 
00199     tabname = argval[10];       /* Input mpup-position table        */
00200   if ((OFFSET = is_set(argval[11])) &&       /* User-defined offset */
00201       sscanf(argval[11], "%lf,%lf",&offset[0],&offset[1]) != 2) { 
00202     print_error("Invalid option -offset dx,dy (got '%s')", argval[11]);
00203     exit_session(ERR_BAD_PARAM);
00204   }
00205 
00206   /* ----- Output ------------------------------ */
00207 
00208   print_msg("o Input mpup frame: %s",pupname);
00209   print_msg("o Nb of bins in histogram: %d",nbin);
00210   print_msg("o Threshold cut-off: %.2f%%",fracthresh*100);
00211   print_msg("o Mpup window half-width: %.2f x %.2f interlens",
00212             fracdist[0],fracdist[1]);
00213   if (ANALYSIS) {
00214     print_msg("o ANALYSIS: Output mpups table: %s",ananame);
00215     print_msg("o Geometric PSF: %d-gaussian description", nggeo);
00216     print_msg("o Global PSF: %d-gaussian description", ngpsf);
00217     if (ngloc < 0) {
00218       ngloc = -ngloc;
00219       DISCARD_GLOBPSF = 1;
00220       print_msg("o Core PSF = Geometric PSF");
00221     }
00222     else {
00223       DISCARD_GLOBPSF = 0;
00224       print_msg("o Core PSF = Geometric PSF (x) Global PSF");
00225     }
00226     print_msg("o Local PSF: %d-gaussian description", ngloc);
00227   }
00228 
00229   /* We do not handle local PSF with more than 1 gaussian... */
00230   if (ABS(ngloc) > 1) {
00231     print_error("Local PSF with more than *one* gaussian not implemented yet!");
00232     exit_session(ERR_NOIMPL);
00233   }
00234 
00235   print_msg("o Reference wavelength: %.1f", lbdaref);
00236 
00237   if (TABLE) print_msg("o Input 1st-order mpup-position table: %s", tabname);
00238   if (TABLE && OFFSET) print_msg("o User-defined offset: %g x %g px", 
00239                                  offset[0],offset[1]);
00240 
00241   /* ===== Micropupils frame ============================== */
00242 
00243   if ((status = open_image(&mpupima,pupname,"I",-1,"micropupils")))
00244     exit_session(status);
00245 
00246   if (!(channel = read_channel(&mpupima))) exit_session(channel);
00247 
00248   if (ANALYSIS && USE_PROFERR) {
00249     RoN2 = read_mean_RoN(&mpupima, channel);
00250     RoN2 *= RoN2;
00251   }
00252 
00253   /* ##### LOCATING MPUP ############################## */
00254 
00255   if (!TABLE) {
00256 
00257     /* ===== Frame statistics ============================== */
00258 
00259     frame_stats(&mpupima,&mean,&sigma,&median,&sigmed);
00260     if (DEBUG) {
00261       print_msg("   Min=%.5f, Max=%.5f",mpupima.min,mpupima.max);
00262       print_msg("   Mean=%.5f (%.5f), Median=%.5f (%.5f)",
00263                 mean,sigma,median,sigmed);
00264     }
00265 
00266     /* ===== Frame histogram ============================== */
00267 
00268     min = 1e-4; /* instead of mpupima.min */
00269     max = MIN(mpupima.max,median+50*sigmed); /* instead of mpupima.max */
00270 
00271     print_msg("Building frame histogram:\n"
00272               "   %d bins in [%.2f,%.2f]",nbin,min,max);
00273     histospec = frame_histo(&mpupima,min,max,nbin);
00274 
00275     if (DEBUG) {
00276       sprintf(tmpname,"%s_frame_histo",DBG_PREFIX);
00277       print_msg("DEBUG: storing frame histogram in %s",tmpname);
00278       sprintf(dident,"Histo of %s",mpupima.name);
00279       disable_erase_flag();
00280       dump_spec_mem(histospec,tmpname,dident,"N",LOSVD_RAW);
00281       restore_erase_flag();
00282     }
00283 
00284     /* ===== Compute threshold ============================== */
00285 
00286     threshold = histo_threshold(histospec,fracthresh);  
00287     print_msg("Threshold for a %.2f%% cut-off: %.2f",fracthresh*100,threshold);
00288 
00289     free_spec_mem(histospec); 
00290     free(histospec);
00291 
00292     /* ===== Detection ============================== */
00293 
00294     npup = detect_mpups(mpupima, threshold, &ppup, &fpup, &xpup, &ypup);
00295 
00296     minmax(fpup,npup,&min,&max,&i,&j);
00297     print_msg("   => %d classes found [%.2f-%.2f]",npup,min,max);
00298 
00299     /* mean_sig(fpup,npup,&mean,&sigma);
00300        print_msg("   => Classe flux: mean = %.2f +- %.2f",mean,sigma); */
00301 
00302     /* ===== Mpups orders ============================== */
00303 
00304     /* Flux threshold between the different orders */
00305     threshold = min + (max - min)*THRESHORDER;
00306     y0 = mpupima.starty + mpupima.ny*mpupima.stepy/2;  /* Frame y-middle */
00307 
00308     print_msg("Threshold between the different orders: "
00309               "flux=%.2f (%.1f), y=%.2f",threshold,THRESHORDER,y0);
00310 
00311     opup = (int *)malloc(npup*sizeof(int));    /* order */
00312     for (i=0; i<npup; i++) {
00313       if (fpup[i] > threshold) opup[i] = 1;  /* 1st order mpups are bright */
00314       else if (ypup[i] > y0)   opup[i] = 0;  /* Oth order mpups are above  */
00315       else                     opup[i] = 2;  /* Others are 2nd order mpups */
00316     }
00317     mpup_orders(fpup,opup,npup,fl,dfl);
00318 
00319     /* ===== Mpups histograms ============================== */
00320 
00321     if (DEBUG) {
00322       disable_erase_flag();
00323       sprintf(tmpname,"%s_mpups_histo",DBG_PREFIX);
00324       print_msg("DEBUG: Building 1st-order mpups histogram %s:\n"
00325                 "   %d bins in ]%.2f,%.2f[",tmpname, 
00326                 nbin/2,fl[1]-5*dfl[1],fl[1]+5*dfl[1]);
00327       histospec = array_histo(fpup,npup,fl[1]-5*dfl[1],fl[1]+5*dfl[1],nbin/2);
00328       sprintf(dident,"Histo of 1st-order mpups (%s)",mpupima.name);
00329       dump_spec_mem(histospec,tmpname,dident,"N",LOSVD_RAW);
00330       free_spec_mem(histospec); 
00331       free(histospec);
00332 
00333       sprintf(tmpname,"%s_mpups_histo02",DBG_PREFIX);
00334       print_msg("DEBUG: Building non-1st-order mpups histogram %s:\n"
00335                 "   %d bins in ]%.2f,%.2f[",tmpname, 
00336                 nbin/2,0.,fl[0]+5*dfl[0]);
00337       histospec = array_histo(fpup,npup,0,fl[0]+5*dfl[0],nbin/2);
00338       sprintf(dident,"Histo of non-1st-order mpups (%s)",mpupima.name);
00339       dump_spec_mem(histospec,tmpname,dident,"N",LOSVD_RAW);
00340       free_spec_mem(histospec); 
00341       free(histospec);
00342       restore_erase_flag();
00343     }
00344 
00345   } /* !TABLE */
00346 
00347   else {
00348 
00349     /* ----- Open ------------------------------ */
00350 
00351     print_msg("Reading 1st-order mpup positions from table %s...", tabname);
00352     if (open_table(&intable,tabname,"I") < 0) {
00353       print_error("Cannot open table %s",tabname);
00354       exit_session(ERR_OPEN);
00355     }
00356     handle_select_flag(&intable,'W',NULL);
00357     npup = intable.row;
00358     print_msg("   => %d positions found",npup);
00359 
00360     /* open :X1 and :Y1 columns */
00361     colxin = get_col_ref(&intable,"X1");
00362     colyin = get_col_ref(&intable,"Y1");
00363     if (colxin<0 || colyin<0) {
00364       print_error("Cannot find columns X1,Y1 in table %s",intable.name);
00365       exit_session(ERR_NOCOL);
00366     }
00367 
00368     /* ----- Read ------------------------------ */
00369 
00370     xpup = (float *)malloc(npup*sizeof(float));
00371     ypup = (float *)malloc(npup*sizeof(float));
00372 
00373     RD_col(&intable,colxin,xpup);
00374     RD_col(&intable,colyin,ypup);
00375 
00376     opup = (int *)malloc(npup*sizeof(int)); /* order */
00377     for (i=0; i<npup; i++) opup[i] = 1;
00378 
00379     if (OFFSET) {
00380       for (i=0; i<npup; i++) {
00381         xpup[i] -= offset[0];
00382         ypup[i] -= offset[1];
00383       }
00384     }
00385 
00386   } /* TABLE */
00387 
00388   /* ===== Mpups table ============================== */
00389 
00390   if (ANALYSIS || DEBUG) {
00391     print_msg("ANALYSIS: Creating mpups table %s...",ananame);
00392 
00393     if (create_table(&table,ananame,-1,-1,'Q',"Lenses") < 0) {
00394       print_error("Cannot create table %s",ananame);
00395       exit_session(ERR_CREAT);
00396     }
00397 
00398     colno = create_col(&table,LAB_COL_NO,INT,'R',"I4", "");   /* NO       */
00399     colx  = create_col(&table,"XCEN", FLOAT,'R',"F9.6","px"); /* centroid */
00400     coly  = create_col(&table,"YCEN", FLOAT,'R',"F9.6","px");
00401     colo  = create_col(&table,"ORDER",INT,  'R',"I3",  "");   /* order    */
00402 
00403     for (i=0; i<npup; i++) {
00404       no = i+1;
00405       WR_tbl(&table,i,colno,&no);
00406       WR_tbl(&table,i,colx,xpup+i);
00407       WR_tbl(&table,i,coly,ypup+i);
00408       WR_tbl(&table,i,colo,opup+i);
00409     }
00410 
00411     if (!TABLE) {
00412       colf = create_col(&table,"FLUX", FLOAT,'R',"F9.6","");  /* flux     */
00413       coln = create_col(&table,"NPIX", INT,  'R',"I9",  "");  /* npix     */
00414 
00415       for (i=0; i<npup; i++) {
00416         ftmp = fpup[i]; WR_tbl(&table,i,colf,&ftmp);
00417         WR_tbl(&table,i,coln,ppup+i);
00418       }
00419     }
00420 
00421     colxf = create_col(&table,"XG",   FLOAT,'R',"F9.6","px"); /* gaussian */
00422     colyf = create_col(&table,"YG",   FLOAT,'R',"F9.6","px"); /* centroid */
00423     colixf= create_col(&table,"IXG",  FLOAT,'R',"F9.6","");   /* intensity*/
00424     coliyf= create_col(&table,"IYG",  FLOAT,'R',"F9.6","");
00425     colsxf= create_col(&table,"SXG",  FLOAT,'R',"F9.6","px"); /* sigma    */
00426     colsyf= create_col(&table,"SYG",  FLOAT,'R',"F9.6","px");
00427     colrxf= create_col(&table,"RMSXG",FLOAT,'R',"F9.6","");   /* RMS      */
00428     colryf= create_col(&table,"RMSYG",FLOAT,'R',"F9.6","");
00429 
00430   }
00431 
00432   if (!TABLE) { free(ppup); free(fpup); }
00433 
00434   /* ===== Mpup profile gaussian-fit ============================== */
00435 
00436   lens_px = (channel==BLUE_CHANNEL)? SNIFSB_LENS_PX:SNIFSR_LENS_PX;
00437   xhalfw = (int)lens_px*fracdist[0];
00438   yhalfw = (int)lens_px*fracdist[1];
00439   xfullw = 2*xhalfw + 1;
00440   yfullw = 2*yhalfw + 1;
00441   print_msg("Micropupils Y-integration for analysis over "
00442             "2x%d+1=%d x 2x%d+1=%d px [%.2f,%.2f*%.1f]",
00443             xhalfw,xfullw,yhalfw,yfullw,fracdist[0],fracdist[1],lens_px);
00444 
00445   xarr = (double *)malloc(xfullw*sizeof(double));
00446   xfit = (double *)malloc(xfullw*sizeof(double));
00447   yarr = (double *)malloc(yfullw*sizeof(double));
00448   yfit = (double *)malloc(yfullw*sizeof(double));
00449 
00450   i1g = (float *)calloc(npup,sizeof(float));
00451   x1g = (float *)calloc(npup,sizeof(float));
00452   s1g = (float *)calloc(npup,sizeof(float));
00453 
00454   if (ANALYSIS || DEBUG) {
00455     print_msg("ANALYSIS: cubes %s_[x|y]prof with raw and gaussian profiles",
00456               ananame);
00457     sprintf(tmpname,"%s_xprof",ananame);
00458     if (create_tiger_frame(&xcube,tmpname,-1,-1,mpupima.stepx,FLOAT,
00459                            ananame,"Mpup X-profile (gaussian)",NULL) < 0) {
00460       print_error("Cannot create cube %s",tmpname);
00461       exit_session(ERR_CREAT);
00462     }
00463     sprintf(tmpname,"%s_yprof",ananame);
00464     if (create_tiger_frame(&ycube,tmpname,-1,-1,mpupima.stepy,FLOAT,
00465                            ananame,"Mpup Y-profile (gaussian)",NULL) < 0) {
00466       print_error("Cannot create cube %s",tmpname);
00467       exit_session(ERR_CREAT);
00468     }
00469 
00470 #ifdef MEGAVERBOSE
00471     print_msg("Mpup [O]:     Ix   (dIx)       X    (dX)    Sx   (dSx)     RMSx"
00472               "     Iy   (dIy)       Y    (dY)    Sy   (dSy)     RMSy");
00473 #endif
00474   }
00475 
00476   reset_print_progress();
00477   nvpup = 0;               /* Number of valid 1st-order mpups */
00478   for (k=0; k<npup; k++) { /* Loop over pupils */
00479     if (!DEBUG) print_progress("Computing gaussian center",100*(k+1)/npup, 5.);
00480     /* Area of analysis */
00481     pixel_frame(&mpupima,xpup[k],ypup[k],&i0,&j0);
00482     i1 = MAX(0,i0-xhalfw);  i2 = MIN(mpupima.nx-1,i0+xhalfw);
00483     j1 = MAX(0,j0-yhalfw);  j2 = MIN(mpupima.ny-1,j0+yhalfw);
00484     coord_frame(&mpupima,i1,j1,&x0,&y0);
00485 
00486     /* X-center fitting: integration along Y ---------- */ 
00487     for (i=i1; i<=i2; i++) {
00488       xarr[i-i1]=0;
00489       for (j=j1; j<=j2; j++) xarr[i-i1] += RD_frame(&mpupima,i,j);
00490       xarr[i-i1] /= (j2-j1+1);               /* Mean instead of sum */
00491     }
00492 
00493     rmsx = gaussian_mpup(xarr,i2-i1+1,x0,mpupima.stepx,xgau,dxgau,xfit);
00494 
00495     /* Store results */
00496     if (rmsx >= 0) {
00497       xpup[nvpup] = xgau[1];
00498       i1g[k]      = xgau[0]; /* Storage for latter on (ANALYSIS) */
00499       x1g[k]      = xgau[1];
00500       s1g[k]      = xgau[2];
00501     }
00502 
00503     if (ANALYSIS || DEBUG) { /* Store profile and fit (zero-ed if error) */
00504       put_tiger_arrays(&xcube,k+1,xarr,xfit,DOUBLE,i2-i1+1,x0);
00505     }
00506 
00507     /* Y-center fitting: integration along X ---------- */ 
00508     for (j=j1; j<=j2; j++) {
00509       yarr[j-j1]=0;
00510       for (i=i1; i<=i2; i++) yarr[j-j1] += RD_frame(&mpupima,i,j);
00511       yarr[j-j1] /= (i2-i1+1);               /* Mean instead of sum */
00512     }
00513     
00514     rmsy = gaussian_mpup(yarr,j2-j1+1,y0,mpupima.stepy,ygau,dygau,yfit);
00515 
00516     /* Store results */
00517     if (rmsy >= 0) ypup[nvpup] = ygau[1];
00518 
00519     if (ANALYSIS || DEBUG) { /* Store profile and fit (zero-ed if error) */
00520       put_tiger_arrays(&ycube,k+1,yarr,yfit,DOUBLE,j2-j1+1,y0);
00521     }
00522 
00523     /* Store results in table */
00524     if (ANALYSIS || DEBUG) {
00525       ftmp = (float)(xgau[0]); WR_tbl(&table,k,colixf,&ftmp);
00526       ftmp = (float)(xgau[1]); WR_tbl(&table,k, colxf,&ftmp);
00527       ftmp = (float)(xgau[2]); WR_tbl(&table,k,colsxf,&ftmp);
00528       ftmp = (float)(rmsx);    WR_tbl(&table,k,colrxf,&ftmp);
00529       ftmp = (float)(ygau[0]); WR_tbl(&table,k,coliyf,&ftmp);
00530       ftmp = (float)(ygau[1]); WR_tbl(&table,k, colyf,&ftmp);
00531       ftmp = (float)(ygau[2]); WR_tbl(&table,k,colsyf,&ftmp);
00532       ftmp = (float)(rmsy);    WR_tbl(&table,k,colryf,&ftmp);
00533     }
00534     
00535     /* Keep only 1st order mpups that have a good shape in X *and* Y */
00536     if (rmsx >= 0 && rmsy >= 0) {
00537 #ifdef MEGAVERBOSE
00538       if (DEBUG) {
00539         print_msg("a%3d [%d]:"
00540                   " %6.3f (%5.0e) %7.2f (%5.0e) %.3f (%5.0e) %7.2e"
00541                   " %6.3f (%5.0e) %7.2f (%5.0e) %.3f (%5.0e) %7.2e",
00542                   k+1,opup[k],
00543                   xgau[0],dxgau[0],xgau[1],dxgau[1],xgau[2],dxgau[2],rmsx,
00544                   ygau[0],dygau[0],ygau[1],dygau[1],ygau[2],dygau[2],rmsy);
00545       }
00546 #endif
00547       if (opup[k] == 1) nvpup++; /* Keep valid 1st order mpups only */
00548     }
00549     else {
00550       if (DEBUG) print_msg("a%3d [%d]: invalid fit (status %c:%d)", 
00551                            k+1,opup[k],(rmsx<0)?'X':'Y',
00552                            (rmsx<0)?ABS((int)rmsx):ABS((int)rmsy));
00553     }
00554 
00555   } /* Loop over pupils */
00556 
00557   /* Memory management */
00558   close_frame(&mpupima);
00559 
00560   free(xarr); free(xfit); free(yarr); free(yfit);
00561   free(opup);
00562 
00563   if (ANALYSIS || DEBUG) {
00564     write_file_class(&xcube, DEBUG_FILE); close_tiger_frame(&xcube);
00565     write_file_class(&ycube, DEBUG_FILE); close_tiger_frame(&ycube);
00566     close_table(&table);
00567   }
00568 
00569   /* Valid 1st-order mpups */
00570   npup = nvpup; 
00571   print_msg("   => Working on %d valid 1st-order mpups",npup);
00572 
00573   /* ##### DISTORTIONS ############################## */
00574 
00575   /* ===== Ordering and triangulating ============================== */
00576 
00577   print_msg("Triangulating...");
00578 
00579   iadj = (int *)malloc((6*npup-9)*sizeof(int));
00580   iend = (int *)malloc(npup*sizeof(int));
00581 
00582   status = 2; /* xpup and ypup are reordered */
00583   reordr(&npup,&status,xpup,ypup,NULL,iend); /* use iend as a temp.array */
00584   trmesh(&npup, xpup, ypup, iadj, iend, &status); 
00585   if (status != 0) {
00586     print_error("Status %d while triangulating", status);
00587     exit_session(UNKNOWN);
00588   }
00589 
00590   /* ===== Triangle definitions ============================== */
00591 
00592   xtri = (double *)malloc(npup*6 *sizeof(double)); /* barycenter */
00593   ytri = (double *)malloc(npup*6 *sizeof(double));
00594   stri = (double *)malloc(npup*6 *sizeof(double)); /* surface */
00595   atri = (double *)malloc(npup*6 *sizeof(double)); /* angle */
00596   ltri = (double *)malloc(npup*12*sizeof(double)); /* short sides */
00597 
00598   ntri = study_triangulation(mpupima, xpup, ypup, npup, iadj, iend, 
00599                              xtri, ytri, stri, ltri, atri);
00600 
00601   print_msg("   => %d unique regular triangles",ntri);
00602 
00603   free(iadj); free(iend);
00604 
00605   /* Median size & angle */
00606 
00607   iwork =    (int *)malloc(2*ntri*sizeof(int));
00608   dwork = (double *)malloc(2*ntri*sizeof(double));
00609 
00610   lenssize = median_sig(ltri,2*ntri,&sigmed,dwork,iwork);
00611   
00612   print_msg("   Lens size (median): %.2f +- %.2f px",lenssize,sigmed);
00613   
00614   angle = median_sig(atri,ntri,&sigmed,dwork,iwork);
00615 
00616   print_msg("   Angle (median):     %.2f +- %.2f deg",
00617             angle/M_PI*180,sigmed/M_PI*180);
00618 
00619   free(iwork); free(dwork);
00620   free(xtri); free(ytri); free(stri); free(atri); free(ltri);
00621 
00622   /* ===== First look at distortion ============================== */
00623 
00624   /* Central mpup */
00625   x0 = median_f(xpup,npup);
00626   y0 = median_f(ypup,npup);
00627   print_msg("   xc,yc (median):     %.2f x %.2f",x0,y0);  
00628 
00629 //  x0 = FMP_XC;
00630 //  y0 = FMP_YC;
00631 //  print_warning("Lens array center is hard-coded to %.2f x %.2f",x0,y0);
00632 
00633   for (mindist=MAXFLOAT, i=0; i<npup; i++) {
00634     ftmp = hypot(xpup[i]-x0,ypup[i]-y0);
00635     if (ftmp < mindist) {
00636       i0 = i;
00637       mindist = ftmp;
00638     }
00639   }
00640   print_msg("   Central lens: m%d (%.2f x %.2f) at %.2f px",
00641             i0+1,xpup[i0],ypup[i0],mindist);
00642 
00643   /* Lens array (i,j)-coordinates */
00644   ipup =   (int *)malloc(npup*sizeof(int));    /* mpup indices */
00645   jpup =   (int *)malloc(npup*sizeof(int));
00646   xthe = (float *)malloc(npup*sizeof(float));  /* theoritical positions */
00647   ythe = (float *)malloc(npup*sizeof(float));
00648   dpup = (float *)malloc(npup*sizeof(float));  /* distance obs./theor. */
00649 
00650   rotate_coord(xpup,ypup,npup,xpup[i0],ypup[i0],-angle,xthe,ythe);
00651   for (i=0; i<npup; i++) {
00652     ipup[i] = rint((xthe[i] - xpup[i0])/lenssize);
00653     jpup[i] = rint((ythe[i] - ypup[i0])/lenssize);
00654   }
00655 
00656   /* ===== Distortion parameters fit ============================== */
00657 
00658   print_msg("Distortion parameter fit (%d param, %d mpups)...",ndpar,npup);
00659 
00660   /* Initial guess */
00661   dpar[0] = lenssize; /* MLA lenssize */
00662   dpar[1] = angle;    /* MLA angle */
00663   dpar[2] = x0;       /* Center of distortion */
00664   dpar[3] = y0;
00665   distor = (channel==BLUE_CHANNEL)?
00666     (SNIFSB_CAM_E + SNIFSB_COLL_E)/SQ(SNIFSB_CAM_F/(CCDB_PIX*1e-4)):
00667     (SNIFSR_CAM_E + SNIFSR_COLL_E)/SQ(SNIFSR_CAM_F/(CCDR_PIX*1e-4));
00668   dpar[4] = (channel==BLUE_CHANNEL)?
00669     SNIFSB_DISTORTION:SNIFSR_DISTORTION;
00670   if (DEBUG) {
00671     print_msg("   Theoritical distortion parameter: %g", distor);
00672     print_msg("   Distortion parameter fixed to: %g", dpar[4]);
00673   }
00674 
00675   /* Normalize fitting parameters */
00676   dpar[0] /= FMP_LENSSNORM;
00677   dpar[1] /= FMP_ANGLENORM;
00678   dpar[2] /= FMP_POSITNORM;
00679   dpar[3] /= FMP_POSITNORM;
00680   dpar[4] /= FMP_ALPHANORM;
00681 
00682 #ifdef USE_R4_DIST
00683   print_warning("Distortion law includes an extra r^4 component (%d)",ndpar);
00684   dpar[5] = SQ(dpar[4]*FMP_ALPHANORM); /* Distortion: beta */
00685   dpar[5] /= SQ(FMP_ALPHANORM);
00686 #endif
00687 
00688   /* Global variables */
00689   glmlaxpup = (double *)malloc(npup*sizeof(double));
00690   glmlaypup = (double *)malloc(npup*sizeof(double));
00691   glmlaml   =    (int *)malloc(npup*sizeof(int));
00692   glmlanl   =    (int *)malloc(npup*sizeof(int));
00693   for (i=0; i<npup; i++) {
00694     glmlaxpup[i] = xpup[i]; /* cannot memcpy because of types float/double */
00695     glmlaypup[i] = ypup[i];
00696   }
00697   memcpy((int *)glmlaml,(int *)ipup,npup*sizeof(int));
00698   memcpy((int *)glmlanl,(int *)jpup,npup*sizeof(int));
00699   glmlaxlc = xpup[i0];
00700   glmlaylc = ypup[i0];
00701 
00702   /* Unbound Least-square fit */
00703   xarr = (double *)calloc(npup,sizeof(double)); /* Minimizing distances: 0 */
00704   xfit = (double *)malloc(npup*sizeof(double));
00705 
00706   edpar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00707 
00708   status = nllsqfit_bnd(NULL, xarr, NULL, npup, dpar, NULL, NULL, ndpar,
00709                         *nllsqfit_distMLA, xfit, NULL, edpar, &chi2,&gof);
00710   
00711   free(xarr); free(xfit);
00712 
00713   /* Renormalize fitting parameters */
00714   dpar[0] *= FMP_LENSSNORM;   edpar[0] *= FMP_LENSSNORM;
00715   dpar[1] *= FMP_ANGLENORM;   edpar[1] *= FMP_ANGLENORM;
00716   dpar[2] *= FMP_POSITNORM;   edpar[2] *= FMP_POSITNORM;
00717   dpar[3] *= FMP_POSITNORM;   edpar[3] *= FMP_POSITNORM;
00718   dpar[4] *= FMP_ALPHANORM;   edpar[4] *= FMP_ALPHANORM;
00719 
00720   print_msg("Status = %d, chi2 = %f",status,chi2);
00721   print_msg("   Lens size: %.2f +- %.2f px",dpar[0],edpar[0]);
00722   print_msg("   Angle:     %.2f +- %.2f deg",
00723             dpar[1]/M_PI*180,edpar[1]/M_PI*180);
00724   print_msg("   xc,yc:     %.2f x %.2f +-(%.2fx%.2f)",
00725             dpar[2],dpar[3],edpar[2],edpar[3]);
00726   print_msg("   alpha:     %.2e +- %.0e", dpar[4],edpar[4]);
00727 
00728 #ifdef USE_R4_DIST
00729   dpar[5] *= SQ(FMP_ALPHANORM);   edpar[5] *= SQ(FMP_ALPHANORM);
00730   print_msg("   beta:      %.2e +- %.0e", dpar[5],edpar[5]);
00731 #endif
00732 
00733   /* ##### MASK ############################## */
00734 
00735   lenssize = dpar[0]; /* MLA lenssize */
00736   angle    = dpar[1]; /* MLA angle */
00737   x0       = dpar[2]; /* Center of distortion */
00738   y0       = dpar[3];
00739 
00740   /* ===== Observed & Undistorted theoritical positions ==================== */
00741 
00742   rotate_coord(xpup,ypup,npup,xpup[i0],ypup[i0],-angle,xthe,ythe);
00743   for (i=0; i<npup; i++) {
00744     /* Update (i,j) coordinates according to new central lens */
00745     ipup[i] = rint((xthe[i] - xpup[i0])/lenssize);
00746     jpup[i] = rint((ythe[i] - ypup[i0])/lenssize);
00747     xthe[i] = xpup[i0] + ipup[i]*lenssize;
00748     ythe[i] = ypup[i0] + jpup[i]*lenssize;
00749   }
00750   rotate_coord(xthe,ythe,npup,xpup[i0],ypup[i0],angle,xthe,ythe);
00751 
00752   /* ===== Distorted theoritical positions ============================== */
00753 
00754   /* Global variables needed by model_distMLA */
00755   memcpy((int *)glmlaml,(int *)ipup,npup*sizeof(int));
00756   memcpy((int *)glmlanl,(int *)jpup,npup*sizeof(int));
00757   glmlaxlc = xpup[i0];
00758   glmlaylc = ypup[i0];
00759 
00760   model_distMLA(npup,dpar,ndpar, xthe,ythe,dpup); /* make use of glob.var. */
00761 
00762   /* ===== Residuals ============================== */
00763 
00764   median = median_f(dpup,npup);
00765   mean_sig_f(dpup,npup,&mean,&sigma);
00766   print_msg("   => Residuals: Med=%.2f px, Mean=%.2f px sig=%.2f px",
00767             median,mean,sigma);
00768 
00769   if (DEBUG) {                               /* Residuals histogram */
00770     dwork = (double *)malloc(npup*sizeof(double));
00771     indx2 =    (int *)malloc(npup*sizeof(int));
00772     for (i=0; i<npup; i++) dwork[i] = dpup[i];
00773     indexx(npup,dwork,indx2);
00774     sprintf(tmpname,"%s_res_histo",DBG_PREFIX);
00775     print_msg("DEBUG: Building residual histogram %s:\n"
00776               "   %d bins in ]%.2f,%.2f[",tmpname,
00777               nbin/4,dwork[indx2[0]],dwork[indx2[npup-1]]);
00778     histospec = array_histo(dwork,npup,
00779                             dwork[indx2[0]],dwork[indx2[npup-1]],nbin/4);
00780     sprintf(dident,"Histo of residuals (%s)",mpupima.name);
00781     disable_erase_flag();
00782     dump_spec_mem(histospec,tmpname,dident,"N",LOSVD_RAW);
00783     restore_erase_flag();
00784     free_spec_mem(histospec); 
00785     free(histospec); free(dwork); free(indx2);
00786   }
00787 
00788   /* ===== Sorting ============================== */
00789 
00790   indx = (int *)malloc(npup*sizeof(int));     /* sorting index */
00791 
00792   /* Sort mpups according to... */
00793   rpup = (double *)malloc(npup*sizeof(double));
00794   for (i=0; i<npup; i++) {
00795     //    rpup[i] = hypot(xpup[i]-x0,ypup[i]-y0);  /* radial distance */
00796     rpup[i] = ipup[i] + jpup[i]*100;         /* Bottom to top, left to right */
00797   }
00798   indexx(npup,rpup,indx);
00799   i0 = indx[0];                              /* Lens #1 */
00800   invert_indexx(indx, npup);
00801   //  print_msg("Sorting from center of distortion: m%d (%.2f x %.2f) at %.2f px -> lens #%d",
00802   //            i0+1,xpup[i0],ypup[i0],rpup[i0],indx[i0]+1);
00803   print_msg("Sorting from bottom-left corner: m%d (%.2f x %.2f) -> lens #%d",
00804             i0+1,xpup[i0],ypup[i0],indx[i0]+1);
00805   free(rpup);
00806   
00807   /* ===== Mask creation ============================== */
00808 
00809   print_msg("Creating mask table '%s'...",maskname);
00810 
00811   if (create_table(&mask,maskname,-1,-1,'Q',"Lenses") < 0) {
00812     print_error("Cannot create table %s",maskname);
00813     exit_session(ERR_CREAT);
00814   }
00815 
00816   colno = create_col(&mask,LAB_COL_NO,     INT,  'R',"I4",  "");
00817   colx  = create_col(&mask,LAB_COL_XLD,    FLOAT,'R',"F9.6","px");
00818   coly  = create_col(&mask,LAB_COL_YLD,    FLOAT,'R',"F9.6","px");
00819   colxf = create_col(&mask,LAB_COL_XTH_LND,FLOAT,'R',"F9.6","px");
00820   colyf = create_col(&mask,LAB_COL_YTH_LND,FLOAT,'R',"F9.6","px");
00821   colf  = create_col(&mask,"RESND",        FLOAT,'R',"F9.6","px");
00822   coli  = create_col(&mask,LAB_COL_I,      INT,  'R',"I3",  "");
00823   colj  = create_col(&mask,LAB_COL_J,      INT,  'R',"I3",  "");
00824 
00825   colxf = create_col(&mask,LAB_COL_XTH_LD,FLOAT,'R',"F9.6","px");
00826   colyf = create_col(&mask,LAB_COL_YTH_LD,FLOAT,'R',"F9.6","px");
00827   colf  = create_col(&mask,"RESD",        FLOAT,'R',"F9.6","px");
00828 
00829   for (i=0; i<npup; i++) {
00830     no = indx[i]+1;
00831     WR_tbl(&mask,i,colno,&no);    /* Lens number            */
00832     WR_tbl(&mask,i,colx, xpup+i); /* Obs. position          */
00833     WR_tbl(&mask,i,coly, ypup+i);
00834     WR_tbl(&mask,i,colxf,xthe+i); /* Theor. pos. (no dist.) */
00835     WR_tbl(&mask,i,colyf,ythe+i);
00836     WR_tbl(&mask,i,colf, dpup+i); /* Residual (px)          */
00837     WR_tbl(&mask,i,coli, ipup+i); /* (i,j)-coordinates      */
00838     WR_tbl(&mask,i,colj, jpup+i);
00839 
00840     WR_tbl(&mask,i,colxf,xthe+i); /* Theor. pos. (w/ dist.) */
00841     WR_tbl(&mask,i,colyf,ythe+i);
00842     WR_tbl(&mask,i,colf, dpup+i); /* Residual (px)          */
00843   }
00844 
00845   /* Mask descriptors */
00846   WR_desc(&mask, CHANNEL, CHAR, lg_ident+1, channel_names[channel]);
00847 
00848   print_msg("Writing distortion parameters (%d) in "
00849             "mask descriptors %s and %s",ndpar,MPUP_NDISTPAR,MPUP_DISTPAR);
00850   if (DEBUG) print_msg("   DISTOR: lenssize,angle,x0,y0,alpha,etc.");
00851   WR_desc(&mask,MPUP_NDISTPAR,INT,1,&ndpar);
00852   WR_desc(&mask,MPUP_DISTPAR, DOUBLE,ndpar,dpar);
00853 
00854   free(glmlaxpup); free(glmlaypup); free(glmlaml); free(glmlanl);
00855 
00856   free(indx);
00857   free(xpup); free(ypup);
00858   free(dpup); free(xthe); free(ythe);
00859   free(ipup); free(jpup); 
00860 
00861 
00862   /* ##### PROFILE ANALYSIS ############################## */
00863 
00864   if (ANALYSIS) {
00865 
00866     print_msg("========== ANALYSIS ==========");
00867 
00868     /* ===== Profile datacube ============================== */
00869 
00870     /* WARNING: the X-profile datacube contains all the profiles that
00871        were correctly adjusted by a single gaussian, including non-1st
00872        order mpups. The associated table contains *all* the mpups. */
00873 
00874     /* Opening X-profile datacube */
00875     sprintf(tmpname,"%s_xprof",ananame);
00876     print_msg("Opening %s X-profile cube",tmpname);
00877     if (open_tiger_frame(&xcube,tmpname,"IO") < 0) {
00878       print_error("Cannot open %s",tmpname);
00879       exit_session(ERR_OPEN);
00880     }
00881 
00882     /* Getting (gaussian-fit) coordinates */
00883     npup = xcube.nbspec;
00884     xpup = (float *)malloc(npup*sizeof(float));  /* Position               */
00885     ypup = (float *)malloc(npup*sizeof(float));
00886     nopup=   (int *)malloc(npup*sizeof(int));    /* Lens number (!= mask)  */
00887     lpup =   (int *)malloc(npup*sizeof(int));    /* Line in analysis table 
00888                                                     nopup[i] = lpup[i]+1   */
00889     opup =   (int *)malloc(npup*sizeof(int));    /* Order */
00890     spup = (float *)malloc(npup*sizeof(float));  /* Dispersion (X-1D gauss) */
00891 
00892     get_lenses_coord(&xcube,"XG","YG",nopup,xpup,ypup,lpup);
00893     get_lenses_col(&xcube,"SXG",spup,lpup);
00894     get_lenses_col(&xcube,"ORDER",opup,lpup); 
00895 
00896     /* ===== Geometric PSF ============================== */
00897 
00898     switch (channel) {
00899     case BLUE_CHANNEL:
00900       scale = (SNIFSB_CAM_F * SNIFS_LENS_F) / 
00901         (SNIFSB_COLL_F * SNIFS_GAMMA * UHT_FOCR * CCDB_PIX*1e-3);
00902       break;
00903     case RED_CHANNEL:
00904       scale = (SNIFSR_CAM_F * SNIFS_LENS_F) / 
00905         (SNIFSR_COLL_F * SNIFS_GAMMA * UHT_FOCR * CCDR_PIX*1e-3);
00906       break;
00907     default:
00908       print_error("Unknown SNIFS spectroscopic channel (%d)",channel);
00909       exit_session(ERR_NOIMPL);
00910       break;
00911     }
00912     print_msg("Geometric PSF scale = %.2f pix", scale);
00913 
00914     glnggeo = nggeo;
00915     glIgeo = (double *)malloc(nggeo*sizeof(double));
00916     glXgeo = (double *)malloc(nggeo*sizeof(double));
00917     glSgeo = (double *)malloc(nggeo*sizeof(double));
00918     geonorm = set_geoPSF(scale); /* Initialize the global variables */
00919     print_msg("Geometric PSF norm = %.2f",geonorm);
00920 
00921     /* ===== Fit of the global PSF ============================== */
00922 
00923     /* ----- Mpups selection ------------------------------ */
00924 
00925     icen = (int *)malloc(npup*sizeof(int)); /* Arr. of central spectra index */
00926 
00927     /* Mpup selection */
00928 
00929     /* The central mpups (i.e. closest to the center of distortion x0,y0) are
00930        not the sharpest ones, and therefore cannot be used to estimate the
00931        core PSF. therefore, we look for the sharpest mpup after median spatial
00932        filtering. */
00933 
00934     /* Sharpest filtered mpup */
00935 
00936     maxdist = lenssize*ABS(FITNEIGHBOR);
00937 
00938     ncen = select_mpup_sharpest_smooth(maxdist,npup,xpup,ypup,opup,spup, icen);
00939 
00940     /* Array of selected spectra length and start */
00941     pcen =    (int *)malloc(ncen*sizeof(int));
00942     scen = (double *)malloc(ncen*sizeof(double));
00943 
00944     /* Array w/ all the concatenated profiles */
00945     xarr = (double *)malloc(ncen*xfullw*sizeof(double));
00946     xfit = (double *)malloc(ncen*xfullw*sizeof(double));
00947 
00948     for (npts=k=i=0; i<ncen; i++) {
00949       get_tiger_spec(&xcube, &spec, NULL, nopup[icen[i]]);
00950       pcen[i] = spec.npts;
00951       scen[i] = spec.start;
00952       npts   += pcen[i];     /* Total nb of points up to now */
00953       for (j=0; j<pcen[i]; j++, k++) xarr[k] = RD_spec(&spec,j);
00954       free_spec_mem(&spec);
00955       if (DEBUG) {
00956         print_msg("   %2d: Lens a%3d (%7.2fx%7.2f), %d px, S1G=%.3f",
00957                   i+1,lpup[icen[i]]+1,xpup[icen[i]],ypup[icen[i]],
00958                   pcen[i],s1g[lpup[icen[i]]]);
00959       }
00960     }
00961     print_msg("   => %d selected 1st-order mpups (total: %d pixels)",
00962               ncen,npts);
00963 
00964     /* ----- Initialization ------------------------------ */
00965 
00966     /* Parameters: 1 (S0) + 2*(ngpsf-1) (Ii+Si) + 2*ncen (In+Xn) */
00967     ncpsfpar = 2*ngpsf-1 + 2*ncen;
00968     cpsfpar = (double *)malloc(ncpsfpar*sizeof(double)); /* parameters */
00969     ecpsfpar= (double *)malloc(ncpsfpar*sizeof(double)); /* errors     */
00970 
00971     /* Initial guesses */
00972 
00973     /* Global PSF */
00974     cpsfpar[0] = 0.4;                /* Inner disp. (Intensity normed to 1) */
00975     for (k=i=1; i<ngpsf; i++, k+=2) {
00976       cpsfpar[k+1] = 3*cpsfpar[k-1];          /* Dispersion */
00977       cpsfpar[k]   = cpsfpar[0]/cpsfpar[k+1]; /* Intensity */
00978     }
00979     /* Individual mpups */
00980     for (i=0; i<ncen; i++) {
00981       cpsfpar[k++] = i1g[lpup[icen[i]]]/2.;        /* Intensity */
00982       cpsfpar[k++] = x1g[lpup[icen[i]]] - scen[i]; /* Position  */
00983     }
00984 
00985     /* Errors on profiles if any */
00986     /* One could use directly variance from preprocessing, but it would require to
00987        store the variance spectra in the xcube, instead of the gaussian-fit. Not sure 
00988        it's worth it... */
00989     if (USE_PROFERR) {
00990       /* The photon noise assumption could be OK for real data. Tested. */
00991       print_warning("Global PSF fit: assuming photon noise (discarding preproc variance)");
00992       xsig = (double *)malloc(ncen*xfullw*sizeof(double));
00993       for (k=0; k<npts; k++) xsig[k] = sqrt(ABS(xarr[k]) + RoN2);
00994     }
00995     else xsig = NULL;
00996 
00997     /* Lower boundaries */
00998     lcpsfpar = (double *)calloc(ncpsfpar,sizeof(double));  /* Everything >0 */
00999     lcpsfpar[0] = 1e-3;                  /* Force central gaussian non-null */
01000     for (j=2*ngpsf, i=0; i<ncen; i++, j+=2) {
01001       lcpsfpar[j-1] = 1e-3;              /* Icen > 0 */
01002       lcpsfpar[j]   = -MAXDOUBLE;        /* Xcen */
01003     }
01004     
01005 
01006     /* Upper boundaries */
01007     ucpsfpar = (double *)calloc(ncpsfpar,sizeof(double));
01008     for (i=0; i<ncpsfpar; i++) ucpsfpar[i] = MAXDOUBLE;
01009 //    for (k=i=1; i<ngpsf; i++, k+=2) {
01010 //      ucpsfpar[k+1] = 3*cpsfpar[2*(ngpsf-1)]; /* Dispersion < 3*max(disp) */
01011 //    }
01012     for (k=i=1; i<ngpsf; i++, k+=2) {
01013       ucpsfpar[k+1] = 5*cpsfpar[k+1];        /* Dispersion < 3*est(disp) */
01014     }
01015 
01016     //#ifdef MEGAVERBOSE
01017     if (DEBUG) {
01018       /* Don't use k++ in print_msg: strange results */
01019       print_msg("Global PSF fit: Initial guesses");
01020       print_msg("   Global PSF:  S0 = %4.2f",cpsfpar[0]);
01021       for (k=i=1; i<ngpsf; i++, k+=2) 
01022         print_msg("      I%d = %4.2f S%d = %4.2f [%4.2f-%4.2f]",
01023                   i,cpsfpar[k],i,cpsfpar[k+1],
01024                   lcpsfpar[k+1],ucpsfpar[k+1]);
01025       print_msg("   Selected mpups:");
01026       for (i=0; i<ncen; i++, k+=2)
01027         print_msg("      a%3d: I = %4.2f X = %.2f",
01028                   lpup[icen[i]]+1,cpsfpar[k],cpsfpar[k+1]);
01029       
01030     }
01031     //#endif
01032 
01033     /* Normalization */
01034     for (j=2*ngpsf, i=0; i<ncen; i++, j+=2) {
01035       cpsfpar[j-1] /= FMP_IPSFNORM;
01036       cpsfpar[j  ]  = (cpsfpar[j]-pcen[i]/2.)/FMP_XPSFNORM;
01037     }
01038 
01039     /* ----- nllsqfit_bnd ------------------------------ */
01040 
01041     /* Global variables */
01042     glngpsf = ngpsf;
01043     glpcen  = (int *)malloc(ncen*sizeof(int));
01044     memcpy((int *)glpcen,(int *)pcen,ncen*sizeof(int));
01045 
01046     ecpsfpar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
01047 
01048     status = nllsqfit_bnd(NULL,xarr,xsig,npts,cpsfpar,lcpsfpar,ucpsfpar,ncpsfpar,
01049                           *nllsqfit_globalPSF, xfit, NULL, ecpsfpar, 
01050                           &chi2,&gof);
01051 
01052     /* Renormalization */
01053     for (j=2*ngpsf, i=0; i<ncen; i++, j+=2) {
01054       cpsfpar[j-1] *= FMP_IPSFNORM;
01055       cpsfpar[j  ]  = cpsfpar[j]*FMP_XPSFNORM + pcen[i]/2.;
01056       ecpsfpar[j-1]*= FMP_IPSFNORM;
01057       ecpsfpar[j  ]*= FMP_XPSFNORM;
01058     }
01059 
01060     /* ----- Results ------------------------------ */
01061 
01062     print_msg("Global PSF fit: status=%d, chi2=%f",status,chi2);
01063     print_msg("   Global PSF:            S0 = %.3f +-%.3f",
01064               cpsfpar[0],ecpsfpar[0]);
01065     for (k=i=1; i<ngpsf; i++, k+=2) 
01066       print_msg("      I%d = %.3f +-%.3f, S%d = %.3f +-%.3f",
01067                 i,cpsfpar[k],ecpsfpar[k],i,cpsfpar[k+1],ecpsfpar[k+1]);
01068     print_msg("   Selected mpups:");
01069     for (i=0; i<ncen; i++, k+=2)
01070       print_msg("      a%3d: I = %.2f +-%.0e, X = %.2f +-%.0e",
01071                 lpup[icen[i]]+1,
01072                 cpsfpar[k],ecpsfpar[k],cpsfpar[k+1],ecpsfpar[k+1]);
01073 
01074     /* Global PSF storage in an easy-to-handle array */
01075     psf_param = (double *)malloc((1+2*ngpsf)*sizeof(double));
01076     psf_param[0] = scale;
01077     for (i=0; i<ngpsf; i++) {
01078       psf_param[2*i+1] = (i==0)? 1.:cpsfpar[2*i-1];
01079       psf_param[2*i+2] = cpsfpar[2*i];
01080     }
01081 
01082     if (DEBUG) {
01083       /* Raw and adjusted profiles */
01084       sprintf(tmpname,"%s_gpsf_%d%d",ananame,nggeo,ngpsf);
01085       print_msg("ANALYSIS: cube %s with raw/adjusted central profiles",
01086                 tmpname);
01087       if (create_tiger_frame(&ycube,tmpname,-1,-1,mpupima.stepx,FLOAT,
01088                              ananame,"Selected mpup X-profile",NULL) < 0) {
01089         print_error("Cannot create cube %s",tmpname);
01090         exit_session(ERR_CREAT);
01091       }
01092 //      for (k=0; k<npts; k++) { /* LOG */
01093 //        xarr[k] = log10(xarr[k]); 
01094 //        xfit[k] = log10(xfit[k]);
01095 //      }
01096       for (k=i=0; i<ncen; k+=pcen[i],i++) {
01097         put_tiger_arrays(&ycube,nopup[icen[i]],xarr+k,xfit+k,
01098                          DOUBLE,pcen[i],scen[i]);
01099       }
01100 
01101       write_file_class(&ycube, DEBUG_FILE); 
01102       close_tiger_frame(&ycube);
01103 
01104       disable_erase_flag();
01105 
01106       /* Geometric PSF - normalized (maximum) to 1 */
01107       sprintf(tmpname,"%s_geopup",DBG_PREFIX);
01108       sprintf(dident,"Geom. PSF %s %d%d",mpupima.name,nggeo,ngpsf);
01109       print_msg("DEBUG: spectrum %s with %d-gaussian geometric PSF",
01110                 tmpname,nggeo);
01111       create_spec(&spec, tmpname, xfullw, -xhalfw, 1., FLOAT, dident, "");
01112       for (fmax=-MAXFLOAT, j=0; j<xfullw; j++) {
01113         ftmp = geoPSF(-xhalfw+j);
01114         WR_spec(&spec,j,ftmp);
01115         fmax = MAX(fmax,ftmp);
01116       }
01117       for (j=0; j<xfullw; j++) WR_spec(&spec,j,RD_spec(&spec,j)/fmax);
01118       write_file_class(&spec, DEBUG_FILE); 
01119       close_spec(&spec);
01120 
01121       /* Global PSF - normalized (maximum) to 1 */
01122       sprintf(tmpname,"%s_psf",DBG_PREFIX);
01123       sprintf(dident,"Global PSF %s %d%d",mpupima.name,nggeo,ngpsf);
01124       print_msg("DEBUG: spectrum %s with %d-gaussian global PSF",
01125                 tmpname,ngpsf);
01126       create_spec(&spec, tmpname, xfullw, -xhalfw, 1., FLOAT, dident, "");
01127       for (fmax=-MAXFLOAT, j=0; j<xfullw; j++) {
01128         ftmp = globalPSF(-xhalfw+j,psf_param);
01129         WR_spec(&spec,j,ftmp);
01130         fmax = MAX(fmax,ftmp);
01131       }
01132       for (j=0; j<xfullw; j++) WR_spec(&spec,j,RD_spec(&spec,j)/fmax);
01133       write_file_class(&spec, DEBUG_FILE); 
01134       close_spec(&spec);
01135 
01136       restore_erase_flag();
01137     }
01138 
01139     /* ----- Partial deconvolution ------------------------------ */
01140 
01141     print_msg("Partial deconvolution of the global PSF");
01142     norm = deconv_globalPSF(ngpsf, psf_param, 0.75);
01143     print_msg("Global PSF norm = %f",norm);
01144 
01145     /* ----- Memory management ------------------------------ */
01146 
01147     free(xarr); free(xfit); 
01148     if (xsig != NULL) free(xsig);
01149     free(glpcen);
01150     free(icen); free(pcen); free(scen);
01151     free(cpsfpar); free(ecpsfpar); free(lcpsfpar); free(ucpsfpar);
01152 
01153     /* ----- Mask descriptors ------------------------------ */
01154 
01155     j = (DISCARD_GLOBPSF)? -ngloc:ngloc; /* temporary */
01156     print_msg("Writing descriptors in mask '%s': \n"
01157               "%s=%f, %s=%d, %s=%d, %s=%d",maskname,MPUP_LBDAPUP,lbdaref,
01158               MPUP_NGGEO,nggeo,MPUP_NGPSF,ngpsf,MPUP_NGLOC,j);
01159     WR_desc(&mask,MPUP_LBDAPUP,FLOAT,1,&lbdaref);
01160     WR_desc(&mask,MPUP_NGGEO,  INT,1,&nggeo);
01161     WR_desc(&mask,MPUP_NGPSF,  INT,1,&ngpsf);
01162     WR_desc(&mask,MPUP_NGLOC,  INT,1,&j);
01163 
01164     /* Temporary: do not send address of global variable as parameter */
01165     xarr = (double *)malloc(3*nggeo*sizeof(double));
01166     for (j=0; j<nggeo; j++) {
01167       xarr[3*j  ] = glIgeo[j]/geonorm; /* Normalization */
01168       xarr[3*j+1] = glXgeo[j]; 
01169       xarr[3*j+2] = glSgeo[j];
01170     }
01171     print_msg("Writing geometric normalized PSF model in descriptor %s "
01172               "(%d elements)",MPUP_GEOPAR,3*nggeo);
01173     if (DEBUG) print_msg("   %s: I0,X0,S0,I1,X1,S1,etc. (including scale)",
01174                          MPUP_GEOPAR);
01175     WR_desc(&mask,MPUP_GEOPAR,DOUBLE,3*nggeo,xarr);
01176     free(xarr);
01177 
01178     print_msg("Writing global PSF model in descriptor %s (%d elements)",
01179               MPUP_PSFPAR,1+2*ngpsf);
01180     if (DEBUG) print_msg("   %s: scale,I0,S0,I1,S1,etc.",MPUP_PSFPAR);
01181     WR_desc(&mask,MPUP_PSFPAR,DOUBLE,(1+2*ngpsf),psf_param);
01182 
01183     /* ===== Global pupil analysis ============================== */
01184 
01185     /* ----- Core PSF ------------------------------ */
01186 
01187     if (DISCARD_GLOBPSF) { /* Do not use the global PSF  */
01188       ngpsf   = 0;         /* Initial ngpsf in argval[7] */
01189       glngpsf = ngpsf;
01190     }
01191 
01192     alloc2d(&glIpsf,nggeo,MAX(ngpsf,1),DOUBLE);
01193     alloc2d(&glSpsf,nggeo,MAX(ngpsf,1),DOUBLE);
01194 
01195     /* If ngpsf=0, psf_param is not needed, since core PSF = geom. PSF then */
01196     norm = set_corePSF(psf_param);
01197     print_msg("Core PSF norm = %f",norm);
01198 
01199     /* ----- Local PSF ------------------------------ */
01200 
01201     glngloc = ngloc;
01202     xloc  = (double *)malloc((1+2*ngloc)*sizeof(double));
01203     dxloc = (double *)malloc((1+2*ngloc)*sizeof(double));
01204 
01205     /* ----- Analysis table ------------------------------ */
01206 
01207     open_table(&table,ananame,"IO");
01208 
01209     colxf = create_col(&table,"XP",  FLOAT,'R',"F9.6","px"); /* local PSF*/
01210     colixf= create_col(&table,"IXP", FLOAT,'R',"F9.6","");   /* intensity*/
01211     colsxf= create_col(&table,"SXP", FLOAT,'R',"F9.6","px"); /* sigma    */
01212     colrxf= create_col(&table,"RMSXP",FLOAT,'R',"F9.6","");  /* RMS      */
01213 
01214     /* ----- Output cube ------------------------------ */
01215 
01216     sprintf(tmpname,"%s_xprof2",ananame);
01217     print_msg("ANALYSIS: cubes %s with raw and fully adjusted profiles",
01218               tmpname);
01219     if (create_tiger_frame(&ycube,tmpname,-1,-1,mpupima.stepx,FLOAT,
01220                            ananame,"Mpup X-profile (final)",NULL) < 0) {
01221       print_error("Cannot create cube %s",tmpname);
01222       exit_session(ERR_CREAT);
01223     }
01224 
01225     /* ----- Loop on mpups ------------------------------ */
01226 
01227 #ifdef MEGAVERBOSE
01228     if (DEBUG)
01229       print_msg("Mpup [S]:"
01230                 "     Ix   (dIx)       X    (dX)    Sx   (dSx)     RMSx");
01231 #endif
01232 
01233     reset_print_progress();
01234     for (i=0; i<npup; i++) {
01235       if (!DEBUG) print_progress("Final pupil analysis", 100*(i+1)/npup, 5.);
01236 
01237       /* Get the raw profile */
01238       get_tiger_spec(&xcube, &spec, NULL, nopup[i]);
01239 
01240       xarr = (double *)malloc(spec.npts*sizeof(double));
01241       xfit = (double *)malloc(spec.npts*sizeof(double));
01242 
01243       for (j=0; j<spec.npts; j++) xarr[j] = RD_spec(&spec,j);
01244       free_spec_mem(&spec);
01245 
01246       /* Initial guess - CAUTION: 0:position, 1:intensity, 2:dispersion */
01247       xloc[0] = x1g[lpup[i]] - spec.start;              /* Position   */
01248       xloc[2] = 0.1;                                    /* Dispersion */
01249       xloc[1] = i1g[lpup[i]]*s1g[lpup[i]]/norm/xloc[2]; /* Intensity  */
01250       for (j=1; j<ngloc; j++) { /* Next gaussians smaller and wider */
01251         xloc[2*j+1] = xloc[2*j-1]/10.;
01252         xloc[2*j+2] = xloc[2*j]  *10.;
01253       }
01254 
01255       rmsx = final_mpup(xarr,spec.npts,spec.start,spec.step,xloc,dxloc,xfit);
01256 
01257       /* Write in table, CAUTION: 0:position, 1:intensity, 2:dispersion */
01258       ftmp = (float)(xloc[0]); WR_tbl(&table,lpup[i], colxf,&ftmp);
01259       ftmp = (float)(xloc[1]); WR_tbl(&table,lpup[i],colixf,&ftmp);
01260       ftmp = (float)(xloc[2]); WR_tbl(&table,lpup[i],colsxf,&ftmp);
01261       ftmp = (float)(rmsx);    WR_tbl(&table,lpup[i],colrxf,&ftmp);
01262 
01263       if (rmsx >= 0) {
01264 #ifdef MEGAVERBOSE
01265         if (DEBUG) {
01266           print_msg("a%3d [%d]:"
01267                     " %8.2f (%5.0e) %7.2f (%5.0e) %.3f (%5.0e) %7.2e",
01268                     i+1,opup[i],
01269                     xloc[1],dxloc[1],xloc[0],dxloc[0],xloc[2],dxloc[2],rmsx);
01270           for (j=1; j<ngloc; j++) {
01271             print_msg("          %6.3f (%5.0e)    -            %.3f (%5.0e)",
01272                       xloc[2*j+1],dxloc[2*j+1],xloc[2*j+2],dxloc[2*j+2]);
01273           }
01274         }
01275 #endif
01276       }
01277       else {
01278         if (DEBUG) print_msg("a%3d [%d]: invalid fit (status %d)",
01279                              i+1,opup[i],(int)-rmsx);
01280       }
01281       
01282       /* Store in output datacube */
01283       put_tiger_arrays(&ycube,nopup[i],xarr,xfit,
01284                        DOUBLE,spec.npts,spec.start);
01285 
01286       free(xarr); free(xfit);
01287     }
01288 
01289     /* ===== Memory management ============================== */
01290 
01291     write_file_class(&ycube, DEBUG_FILE); close_tiger_frame(&ycube);
01292     close_table(&table);
01293 
01294     free2d(&glIpsf,DOUBLE); free2d(&glSpsf,DOUBLE);
01295     free(xloc); free(dxloc);
01296 
01297     free(psf_param);
01298     free(glIgeo); free(glXgeo); free(glSgeo);
01299 
01300     close_tiger_frame(&xcube);
01301     free(xpup); free(ypup); free(lpup); free(nopup); free(opup); free(spup);
01302 
01303   } /* ANALYSIS */
01304 
01305   /* ##### CONCLUSIONS ############################## */
01306 
01307   write_file_class(&mask, MSK_FILE); 
01308   close_table(&mask);
01309 
01310   free(i1g); free(x1g); free(s1g); 
01311 
01312   exit_session(OK);
01313   return(OK);
01314 }
01315 

Generated on Wed Oct 26 23:59:39 2005 for Snifs by doxygen 1.3.5