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

Generated on Tue Nov 23 18:04:20 2004 for Snifs by doxygen 1.3.3