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

extract/source/prepare_mask.c

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

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