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

focus/source/focus_spectro.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00016 /* =========================================================== */
00017 
00018 #include <IFU_io.h>
00019 #include <IFU_math.h>
00020 #include <snifs.h>
00021 #include <focus.h>
00022 #include "../../extract/incl/extract.h" /* Probably not the right way... */
00023 #include <gsl/gsl_statistics.h>
00024 
00025 #ifdef HAVE_LIBDISLIN
00026 #include <dislin.h>
00027 #endif
00028 
00029 /* === Doxygen Comment ======================================= */
00085 /* =========================================================== */
00086 
00087 int main(int argc, char **argv)
00088 {
00089   int const NARCMAX=20;        /* Maximum number of arcs in input catalogue */
00090   int const ARCWIDTH=128; /* Full-width of arc raster used for correlations */
00091 
00092   IMAGE2D arcframe;
00093   char **argval, **arglabel;
00094   char *catname, *devname, arcname[lg_name+1];
00095   int CAT, MISSSPECFOC=0, ASCII, SEQ;
00096   int i,j;
00097   int status, narc, foc_start, foc_step;
00098   float specfoc, minfoc, maxfoc, range;
00099   double **arctab;
00100   double parpeak[6],eparpeak[6];
00101   double peakx[NARCMAX], *peaky, *peakdy, peakf[NARCMAX];
00102   double peakI[NARCMAX],peakS[NARCMAX],peakSx[NARCMAX],peakSy[NARCMAX];
00103   double peakdI[NARCMAX],peakdS[NARCMAX],peakdSx[NARCMAX],peakdSy[NARCMAX];
00104   double par[3],epar[3];
00105   double msig,emsig,chi2,gof,xnorm;
00106 #ifdef HAVE_LIBDISLIN
00107   Plot plot;
00108   char title[4*lg_ident+1];
00109   float *x,*y,*dy;
00110 #endif
00111 
00112   set_purpose("Compute spectrograph focus from a catalogue of arc exposures");
00113   set_arglist("-in|cat focus.cat -dev ASCII -seq NULL -range 3.");
00114 
00115   init_snifs("$Name:  $");
00116   init_session(argv,argc,&arglabel,&argval);
00117 
00118   if (DEBUG) {
00119     print_msg("$Id: focus_spectro.c,v 1.12 2004/11/12 16:24:37 ycopin Exp $");
00120     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00121   }
00122 
00123   /* ===== Parameters ============================== */
00124 
00125   catname = argval[0];                       /* Input file/catalogue */
00126 
00127   if ((CAT = (strstr(catname,".cat") != NULL))) {
00128     print_msg("o Arc catalog: %s",catname);
00129     RD_catalog(catname, arcname);
00130   } 
00131   else {
00132     strcpy(arcname, catname);
00133     print_msg("o Arc name: %s",arcname);
00134   }
00135 
00136   devname = argval[1];                       /* Output device */
00137   upper_strg(devname);
00138   ASCII = !strcmp(devname,"ASCII");
00139   if (ASCII) print_msg("o ASCII output");
00140   else {
00141 #ifdef HAVE_LIBDISLIN
00142     print_msg("o Output DISLIN device: %s",devname);
00143 #else
00144     print_error("DISLIN output not supported in current executable");
00145     exit_session(ERR_BAD_PARAM);
00146 #endif
00147   }
00148 
00149   if ((SEQ = is_set(argval[2])) &&           /* Focus sequence */
00150       sscanf(argval[2],"%d,%d",&foc_start,&foc_step) != 2) {
00151     print_error("Invalid option -seq start,step (has '%s')",argval[2]);
00152     exit_session(ERR_BAD_PARAM);
00153   }
00154   else print_msg("o Default linear focus sequence: start = %d, step = %d",
00155                  foc_start, foc_step);
00156 
00157   get_argval(3,"%f",&range);                 /* Focus-range quality level */
00158 
00159   /* ===== Loop on catalog (if needed) ============================== */
00160 
00161   alloc2d(&arctab,ARCWIDTH,ARCWIDTH,DOUBLE);
00162   narc = 0;                                  /* Effective nb of arcs */
00163 
00164   do {                                       /* Loop over arc frames */
00165 
00166     /* Open image */
00171     if ((status = open_image(&arcframe,arcname,"I",-1,"arc"))) 
00172       continue;
00173 
00174     /* Extract central raster */
00175     if (extract_arc(&arcframe, ARCWIDTH, arctab)) { 
00176       close_frame(&arcframe);
00177       continue;
00178     }
00179 
00180     /* Read spectrograph position */
00181     specfoc = read_focus(&arcframe);
00182 
00183     if (specfoc < 0) {                       /* Frame number */
00184       MISSSPECFOC = 1;
00185       specfoc = ABS(specfoc);
00186     } 
00187     else {                                   /* Focus position */
00188       if (SEQ && !narc) {
00189         print_warning("Cannot force focus position, option -seq discarded");
00190         SEQ = 0;
00191       }
00192       
00193       if (MISSSPECFOC) {
00194         print_error("Cannot mix focus position and frame number");
00195         close_frame(&arcframe);
00196         exit_session(UNKNOWN);
00197       }
00198     }
00199 
00200     /* Spectrograph position normalization by the first spectrograph position */
00201     if (!narc) {                             /* First spectrograph position */
00202       if (specfoc) {                         /* Non-zero position */
00203         xnorm = specfoc;
00204         if (DEBUG) print_msg("Spectrograph position normalization: %.2f",xnorm);
00205       }
00206       else {                                 /* Zero position */
00207         print_warning("Cannot properly normalize spectrograph position");
00208         xnorm = 1;
00209       }
00210     }
00211     peakx[narc] = specfoc/xnorm;             /* Normalization */
00212 
00213     /* Close arc frame */
00214     close_frame(&arcframe);
00215 
00216     /* Compute and fit auto-correlation peak */
00217     if (DEBUG) print_msg("Computing and fitting arc "
00218                          "auto-correlation peak (%d)",narc+1);
00219 
00220     if (Fit_Xpeak(arctab,arctab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,
00221                   parpeak,eparpeak))
00222       continue;
00223 
00224     /* Convert correlation quantities to intrinsic quantities */
00225     comp_mean_disp(parpeak, eparpeak, &msig, &emsig);
00226 
00227     if (DEBUG) print_parpeak(parpeak, eparpeak, msig, emsig);
00228 
00229     /* Store results in growing arrays */
00230     peakI[narc]  = parpeak[2]; peakdI[narc]  = eparpeak[2];
00231     peakSx[narc] = parpeak[3]; peakdSx[narc] = eparpeak[3];
00232     peakSy[narc] = parpeak[4]; peakdSy[narc] = eparpeak[4];
00233     peakS[narc]  = msig;       peakdS[narc]  = emsig;
00234 
00235     narc++;
00236 
00237     if (narc == NARCMAX) {
00238       print_warning("Cannot handle more than %d arcs, keep going",NARCMAX);
00239       break;
00240     }
00241 
00242   } while (CAT && RD_catalog(catname, arcname)); /* Loop over arc frames */
00243 
00244   free2d(&arctab,DOUBLE);
00245 
00246   /* ===== Focus Output ============================== */
00247 
00248   print_msg(" #   Focus    Intens     SigX     SigY    <Sig>   d<Sig>");
00249   for (i=0; i<narc; i++) 
00250     print_msg("%2d %7.0f %9.3g %8.3f %8.3f %8.3f %8.2g",
00251               i+1,peakx[i]*xnorm,
00252               peakI[i],peakSx[i],peakSy[i],peakS[i],peakdS[i]);
00253 
00254   /* Which quantity to be *minimized*? */
00255   peaky  = peakS; peakdy = peakdS;
00256   print_msg("Focus indicator: <Sig> = RMS(SigX,SigY)");
00257 
00258   if (!CAT) {
00259     /* Minimal output */
00260     if (!VERBOSE) {
00261       printf("%f\n",peaky[0]);
00262       fflush(stdout);
00263     }
00264     exit_session(OK);
00265   }  
00266 
00267 #ifdef HAVE_LIBDISLIN
00268   if (!ASCII && CAT) {
00269     plot_initialize(&plot,!ASCII,devname,NULL,NULL,0,0);
00270     
00271     /* Conversion from double to float */
00272     x = (float *)malloc(narc*sizeof(float));
00273     y = (float *)malloc(narc*sizeof(float));
00274     dy= (float *)malloc(narc*sizeof(float));
00275     for (i=0; i<narc; i++) {
00276       x[i] = peakx[i]*xnorm; 
00277       if (SEQ) x[i] = foc_start + x[i]*foc_step;
00278       y[i]  = peaky[i]; 
00279       dy[i] = peakdy[i];
00280     }
00281 
00282     plot_set_limits(x,y,narc,&(plot.Limits));
00283     plot_increase_limits(&(plot.Limits));
00284 
00285     metafl(plot.devname);                    /* Output device */
00286     disini();
00287     if (DEBUG) errmod("ALL","ON");
00288     else errmod("ALL","OFF");
00289     complx(); texmod("ON");                  /* TeX support */
00290     labdig(2,"Y");                   /* 2 digits on Y-labels (instead of 1) */
00291     if (VERBOSE) paghdr("$Name:  $", "",1,0); /* Technical header */
00292     plot_axes(&plot,"Focus","<\\sigma >");
00293     sprintf(title,"Spectrograph Focus\nCatalogue: %s\n",catname);
00294     chncrv("BOTH");                 /* Automatic change of curve attributes */
00295     plot_array_err(&(plot.Symbol),x,y,narc,NULL,dy);
00296   }
00297 #endif
00298 
00299   /* ===== Focus fit ============================== */
00300 
00301   /* Initial guesses */
00302   if ( peaky[1] < peaky[0] ) {               /* Concave parabola */
00303     i = gsl_stats_min_index(peaky,1,narc);
00304     j = gsl_stats_max_index(peaky,1,narc);
00305   }
00306   else {                                     /* Convex parabola */
00307     i = gsl_stats_max_index(peaky,1,narc);
00308     j = gsl_stats_min_index(peaky,1,narc);
00309   }
00310   if (peakx[i] == peakx[j]) {
00311     print_error("All the frames seem to have the same focus, aborting");
00312     exit_session(UNKNOWN);
00313   }
00314   par[1] = peakx[i];                         /* Extremum position */
00315   par[2] = peaky[i];                         /* Extremum height */
00316   par[0] = (peaky[j]-par[2])/SQ(peakx[j] - par[1]); /* MinMax parabola */
00317 
00318   if (DEBUG) {
00319     print_msg("Focus guess: %g px at focus = %f [curv=%g]",
00320               par[2],par[1],par[0]);
00321 #ifdef HAVE_LIBDISLIN
00322     if (!ASCII && CAT) {
00323       for (i=0; i<narc; i++) y[i] = par[0]*SQ(x[i]/xnorm-par[1]) + par[2];
00324       curve(x,y,narc);                       /* Estimated parabola */
00325     }
00326 #endif
00327   }
00328 
00329   epar[0] = 0; /* Eastern egg: Print-out: 1, Derivatives: 2 */
00330 
00331   status = nllsqfit_bnd(peakx, peaky, peakdy, narc, par, NULL, NULL, 3,
00332                         *nllsqfit_quad, peakf, NULL, epar, &chi2,&gof);
00333 
00334   par[0] /= SQ(xnorm); epar[0] /= SQ(xnorm); /* Renormalize */
00335   par[1] *= xnorm;     epar[1] *= xnorm;
00336 
00337   minfoc = par[1] - sqrt(range/100.*par[2]/par[0]); /* Focus range */
00338   maxfoc = par[1] + sqrt(range/100.*par[2]/par[0]);
00339 
00340   if (SEQ) {
00341     print_msg("Converting frame number to focus position (%d,%d)",
00342               foc_start,foc_step);
00343     par[1] =  par[1]*foc_step + foc_start;
00344     epar[1]= epar[1]*foc_step;
00345     minfoc =  minfoc*foc_step + foc_start;
00346     maxfoc =  maxfoc*foc_step + foc_start;
00347   }
00348 
00349   print_msg("Focus fit: status=%d, Chi2=%g, GoF=%g",status,chi2,gof);
00350   print_msg("Focus position: %.2f +- %.2f",par[1],epar[1]);
00351   print_msg("Best focus: %.2f +- %.2f px",par[2],epar[2]);
00352   print_msg("Focus range [%.0f%%]: %.2f - %.2f",range,minfoc,maxfoc);
00353   if (DEBUG) print_msg("Focus curvature: %g +- %g",par[0],epar[0]);
00354 
00355   /* Minimal output */
00356   if (!VERBOSE) {
00357     printf("%f\n",par[1]);
00358     fflush(stdout);
00359   }
00360 
00361 #ifdef HAVE_LIBDISLIN
00362   if (!ASCII && CAT) {
00363     float x0=par[1],y0=par[2],dx0=epar[1],dy0=epar[2];
00364 
00365     for (i=0; i<narc; i++) y[i] = peakf[i]; /* Conversion from double to float */
00366 
00367     curve(x,y,narc);                         /* Adjusted parabola */
00368 
00369     plot.Symbol.shape = 21;                  /* Best focus */
00370     plot_array_err(&(plot.Symbol),&x0,&y0,1,&dx0,&dy0);
00371 
00372     color("BLUE");
00373     rline(minfoc,par[2]*(1+range/100.),    /* Focus range */
00374           maxfoc,par[2]*(1+range/100.));
00375     color("FORE");
00376 
00377     free(x); free(y); free(dy);
00378 
00379     sprintf(title,"%sBest focus: %.2f \\pm  %.2f (%.2f \\pm  %.2f)\n",
00380             title,x0,dx0,y0,dy0);
00381     plot_title(title);
00382 
00383     if (strcmp(devname,"CONS")) print_msg("Output file: %s",getfil());
00384 
00385     disfin();
00386   }
00387 #endif
00388 
00389   /* Test if optimal focus is inside the probbed range */
00390   if (par[1]<minfoc || par[1]>maxfoc) {
00391     print_error("Focus position (%.2f) outside %.0%%-range [%.0f,%.0f]",
00392                 par[1],range,minfoc,maxfoc);
00393     status = 5;                              /* Unused by nllsqfit_bnd */
00394   }
00395 
00396   exit_session(status);
00397   return(status);
00398 }

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