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

extract/lib/proc_plot.c

Go to the documentation of this file.
00001 /* === Doxygen Comment ======================================= */
00013 /* =========================================================== */
00014 
00015 #include <IFU_io.h>
00016 #include <IFU_math.h>
00017 #include <snifs.h>
00018 #include <extract.h>
00019 #include <gsl/gsl_statistics.h>
00020 
00021 #ifdef HAVE_LIBDISLIN
00022 #include <dislin.h>
00023 #endif
00024 
00026 #define NCOORD 100
00027 
00029 #define MAXNBIN 100
00030 
00032 #define NTICKS  6
00033 
00034 /* === Doxygen Comment ======================================= */
00049 /* =========================================================== */
00050 
00051 void plot_set_limits(const float x[],const float y[], int n, PlotLimits *limits)
00052 {
00053   float xmin,xmax,ymin,ymax;
00054 
00055   gsl_stats_float_minmax(&xmin,&xmax,x,1,n);
00056   gsl_stats_float_minmax(&ymin,&ymax,y,1,n);
00057 
00058   if (!limits->xmin && !limits->xmax && !limits->ymin && !limits->ymax) {
00059     /* Limits have never been set: set them */
00060     limits->xmin = xmin;
00061     limits->xmax = xmax;
00062     limits->ymin = ymin;
00063     limits->ymax = ymax;
00064   }
00065   else {
00066     /* Limits have *already* been set: update them */
00067     limits->xmin = MIN(limits->xmin,xmin);
00068     limits->xmax = MAX(limits->xmax,xmax);
00069     limits->ymin = MIN(limits->ymin,ymin);
00070     limits->ymax = MAX(limits->ymax,ymax);
00071   }
00072 }
00073 
00074 /* === Doxygen Comment ======================================= */
00078 /* =========================================================== */
00079 
00080 void plot_increase_limits(PlotLimits *limits)
00081 {
00082   double delta;
00083   
00084   if (limits->ymin == limits->ymax) delta = 1;
00085   else delta = (limits->ymax - limits->ymin) * 0.1;
00086   limits->ymin -= delta;
00087   limits->ymax += delta;
00088 
00089   if (limits->xmin == limits->xmax) delta = 1;
00090   else delta = (limits->xmax - limits->xmin) * 0.1;
00091   limits->xmin -= delta;
00092   limits->xmax += delta;
00093 }
00094 
00095 /* === Doxygen Comment ======================================= */
00104 /* =========================================================== */
00105 
00106 int inside_plot_limits(float x,float y, PlotLimits limits)
00107 {
00108   return( (x>=limits.xmin && x<=limits.xmax && 
00109            y>=limits.ymin && y<=limits.ymax)? TRUE:FALSE );
00110 }
00111 
00112 /* === Doxygen Comment ======================================= */
00124 /* =========================================================== */
00125 
00126 int plot_compticks(int nticks, float min, float max, float *first, float *step)
00127 {
00128   float delta, fact;
00129   int expo;
00130   
00131   delta = (max - min)/nticks;
00132   /* delta = a * 10^expo = a * fact, w/ 1<=a<10 */
00133   expo = floor(log10(delta));
00134   fact = pow(10,expo);  
00135 
00136   switch ((int)rint(3*log10(delta/fact))) {
00137   case 0: *step = 1.; break;
00138   case 1: *step = 2.; break;
00139   case 2: *step = 5.; break;
00140   case 3: *step = 10; break;
00141   }
00142   *step *= fact;
00143 
00144   *first = (*step)*(floor(min/(*step))+1);
00145   return(floor((max-(*first))/(*step)));
00146 }
00147 
00148 /* === Doxygen Comment ======================================= */
00162 /* =========================================================== */
00163 
00164 void plot_start(Plot *plot, float scale)
00165 {
00166 #ifdef HAVE_LIBDISLIN
00167   float xsize,ysize, dx,dy, axscale;
00168 
00169   metafl(plot->devname);                     /* Output device */
00170   filmod("delete");                          /* Overwrite output file */
00171   setpag("DA4P");                            /* Portrait */
00172   xsize = ABS(scale)*2100;
00173   ysize = ABS(scale)*2970;                   /* Window size [DISLIN units] */
00174 
00175   disini();                                  /* DISLIN initialization: level 0 -> 1 */
00176 
00177   if (DEBUG) errmod("CHECK","OFF");          /* No check on boundaries */
00178   else       errmod("ALL","OFF");
00179 
00180   if (VERBOSE) paghdr("$Name:  $", "",1,0);  /* Technical header */
00181 
00182   /* Axes lengthes: by default, resolution is 100/cm (e.g. 2100 x 2970 for
00183      DA4P), and axe lengthes is 2/3 of page (e.g. 1400 x 1980) */
00184   if (scale<0) {
00185     axscale = MIN(xsize/(dx = plot->Limits.xmax - plot->Limits.xmin),  /* X-scale */
00186                   ysize/(dy = plot->Limits.ymax - plot->Limits.ymin)); /* Y-scale */
00187     xsize = dx*axscale;
00188     ysize = dy*axscale;
00189   }
00190   axslen(xsize,ysize);
00191   
00192   center();                                  /* Use axspos(250,2700); to set position */
00193   labtyp("VERT","Y");
00194   intax();                                   /* Integer labels */
00195   /* ticpos("REVERS", "XY"); */              /* Ticks inward */
00196 
00197   setvlt("SPEC");                            /* A LUT brighter than RAIN */
00198   /* Fonts: disalf,simplx,complx,duplx,triplx,gothic,serif,helve,helves */
00199   complx(); 
00200 
00201   texmod("ON");                              /* TeX support */
00202 
00203 #else
00204   print_error("Function plot_start requires DISLIN support");
00205   /* C99 __func__ (gcc 3.0) */
00206   //  print_error("Function %s requires DISLIN support",__func__);
00207   exit_session(UNKNOWN);
00208 #endif
00209 }
00210 
00211 /* === Doxygen Comment ======================================= */
00222 /* =========================================================== */
00223 
00224 void plot_axes(const Plot *plot, char *xlabel,char *ylabel)
00225 {
00226 #ifdef HAVE_LIBDISLIN
00227 
00228   float xfirst,xstep,yfirst,ystep;
00229   
00230   name(xlabel,"x"); name(ylabel,"y");
00231 
00232   plot_compticks(NTICKS, plot->Limits.xmin,plot->Limits.xmax, &xfirst,&xstep);
00233   plot_compticks(NTICKS, plot->Limits.ymin,plot->Limits.ymax, &yfirst,&ystep);
00234 
00235   if (plot->mincut || plot->maxcut) {        /* 3D plot */
00236     graf3(plot->Limits.xmin,plot->Limits.xmax,xfirst,xstep,
00237           plot->Limits.ymin,plot->Limits.ymax,yfirst,ystep,
00238           plot->mincut,plot->maxcut,plot->mincut,
00239           (plot->maxcut-plot->mincut)/10.);
00240   }
00241   else {                                     /* Pure 2D plot */
00242     graf(plot->Limits.xmin,plot->Limits.xmax,xfirst,xstep,
00243          plot->Limits.ymin,plot->Limits.ymax,yfirst,ystep);
00244   }
00245 #else
00246   print_error("Function plot_axes requires DISLIN support");
00247   exit_session(UNKNOWN);
00248 #endif
00249 }
00250 
00251 /* === Doxygen Comment ======================================= */
00260 /* =========================================================== */
00261 
00262 void plot_title(char *string)
00263 {
00264 #ifdef HAVE_LIBDISLIN
00265   char *copy, *chunk;
00266   int i;
00267   
00268   copy = strdup(string);                     /* Copy the input string */
00269   for (i=0; i<4; i++)                        /* Split the input title */
00270     if ((chunk = strsep(&copy,"\n"))!=NULL && strlen(chunk))
00271       titlin(chunk,(i? (i+1):-1));           /* 1st title underlined */
00272 
00273   title();                                   /* Write the title */
00274   free(copy);
00275 #else
00276   print_error("Function plot_title requires DISLIN support");
00277   exit_session(UNKNOWN);
00278 #endif
00279 }
00280 
00281 /* === Doxygen Comment ======================================= */
00295 /* =========================================================== */
00296 
00297 void plot_set_bkgnd(IMAGE2D *bkgndima, Plot *plot, float mincut,float maxcut,
00298                     int *nx,int *ny, float ***zmat)
00299 {
00300   SPECTRUM *histospec;
00301   PlotLimits lim;
00302 
00303   int i,j;
00304   int imin,imax,jmax,jmin, nbin;
00305   float min,max;
00306 
00307   /* Update plot limits */
00308 
00309   lim = plot->Limits;                        /* Just to shorten the code */
00310 
00311   /* Restrict graph to background frame */
00312   lim.xmin = MAX(bkgndima->startx,lim.xmin);
00313   lim.xmax = MIN(bkgndima->endx  ,lim.xmax);
00314   lim.ymin = MAX(bkgndima->starty,lim.ymin);
00315   lim.ymax = MIN(bkgndima->endy  ,lim.ymax);
00316 
00317   pixel_frame(bkgndima,lim.xmin,lim.ymin,&imin,&jmin);
00318   pixel_frame(bkgndima,lim.xmax,lim.ymax,&imax,&jmax);
00319   coord_frame(bkgndima,imin,jmin,&(lim.xmin),&(lim.ymin));
00320   coord_frame(bkgndima,imax,jmax,&(lim.xmax),&(lim.ymax));
00321 
00322 //#ifdef HAVE_LIBDISLIN
00323 //  /* Set DISLIN coordinates... */
00324 //  sursze(lim.xmin,lim.xmax,lim.ymin,lim.ymax);
00325 //#endif
00326 //  /* ... but add half a pixel on each side to cover whole peripheric pixels */
00327 //  lim.xmin -= 0.5*bkgndima->stepx;
00328 //  lim.xmax += 0.5*bkgndima->stepx;
00329 //  lim.ymin -= 0.5*bkgndima->stepy;
00330 //  lim.ymax += 0.5*bkgndima->stepy;
00331 
00332   plot->Limits = lim;
00333 
00334   *nx = imax-imin+1;
00335   *ny = jmax-jmin+1;
00336 
00337   print_msg("Updating limits according to %s: [%.1f,%.1f:%.1f,%.1f] (%dx%d)",
00338             bkgndima->name, lim.xmin,lim.ymin,lim.xmax,lim.ymax, *nx,*ny);
00339 
00340   /* Dump frame in 2D-matrix */
00341 
00342   alloc2d(zmat,*nx,*ny,FLOAT);
00343   for (min=MAXFLOAT, max=-MAXFLOAT, i=0; i<*nx; i++) {
00344     for (j=0; j<*ny; j++) {
00345       (*zmat)[i][j] = RD_frame(bkgndima,i+imin,j+jmin);
00346       min = MIN(min,(*zmat)[i][j]);
00347       max = MAX(max,(*zmat)[i][j]);
00348     }
00349   }
00350   
00351   /* Set cuts */
00352 
00353   nbin = MIN(MAXNBIN,(int)((*nx)*(*ny)/10));
00354 
00355   print_msg("Building %s histogram: %d bins in [%g,%g]",
00356             bkgndima->name,nbin,min,max);
00357   histospec = subframe_histo(bkgndima,imin,imax,jmin,jmax,min,max,nbin);
00358 
00359   min = histo_threshold(histospec,mincut/100.);
00360   max = histo_threshold(histospec,maxcut/100.);
00361   print_msg("   [%.1f%%,%.1f%%] cut-off: [%g,%g],",mincut,maxcut,min,max);
00362 
00363   plot->mincut = min;
00364   plot->maxcut = max;
00365 
00366   free_spec_mem(histospec); 
00367   free(histospec);
00368 }
00369 
00370 
00371 /* === Doxygen Comment ======================================= */
00379 /* =========================================================== */
00380 
00381 void plot_max(const Maxima_Set *maxset, const Plot *plot, const double shift[])
00382 {
00383   char *xcoord, *ycoord;
00384   int i,j,counter,ASCII;
00385   float x,y;
00386 
00387   print_msg("Plotting max...");
00388 
00389   ASCII = !plot->Graph;
00390 
00391   if (ASCII) {
00392 
00393     /* ===== Initialization (ASCII) ============================== */
00394   
00395     xcoord = (char *)malloc(10*(NCOORD+1)*sizeof(char));
00396     ycoord = (char *)malloc(10*(NCOORD+1)*sizeof(char));
00397 
00398     strcpy(xcoord,"-xdata {");
00399     strcpy(ycoord,"-ydata {"); 
00400   }
00401 #ifdef HAVE_LIBDISLIN
00402   else {
00403     hsymbl(plot->Symbol.size);
00404     setclr(plot->Symbol.color);
00405   }
00406 #endif
00407     
00408   /* ===== Plotting max ============================== */
00409 
00410   for (counter=i=0; i<maxset->nb_ycoords; i++) { /* loop over y */
00411     y = maxset->line[i].ycoord;
00412     if (shift != NULL) y -= shift[1];        /* Apply offset */
00413     for (j=0; j<maxset->line[i].nb_max; j++) { /* loop over x */
00414       x = maxset->line[i].maxima[j].xcoord;
00415       if (shift != NULL) x -= shift[0];      /* Apply offset */
00416 
00417       if (inside_plot_limits(x,y,plot->Limits)) {
00418 
00419         if (ASCII) {
00420 
00421           /* ----- Plain ASCII output ------------------------------ */
00422 
00423           sprintf(xcoord,"%s %g",xcoord,x);
00424           sprintf(ycoord,"%s %g",ycoord,y);
00425           counter ++;
00426 
00427           if (counter%NCOORD == 0) { /* Display & reinitialize strings */
00428             printf("%s} %s}\n",xcoord,ycoord);
00429             fflush(stdout);
00430 
00431             strcpy(xcoord,"-xdata {");
00432             strcpy(ycoord,"-ydata {");
00433           }
00434         }
00435 #ifdef HAVE_LIBDISLIN
00436         else {
00437 
00438           /* ----- Graphical output ------------------------------ */
00439 
00440           rlsymb(plot->Symbol.shape,x,y);
00441         }
00442 #endif
00443       } 
00444         
00445     } /* x */
00446   } /* y */
00447 
00448   if (ASCII) {
00449 
00450     /* ===== Conclusion (ASCII) ============================== */
00451 
00452     if (counter%NCOORD != 0) { /* Display & reinitialize strings */
00453       printf("%s} %s}\n",xcoord,ycoord);
00454       fflush(stdout);
00455     }
00456 
00457     free(xcoord); free(ycoord);
00458   }
00459 }
00460 
00461 /* === Doxygen Comment ======================================= */
00482 /* =========================================================== */
00483 
00484 void plot_lens(const SnifsOptics *optics, int nopup,float xppup,float yppup, 
00485                int nlbda,const float lbda[], int order,
00486                const Plot *plot, float lbdaref,float pixsize,
00487                int dxnc, double *dxcoeff, const double shift[],
00488                double *xmla,double *ymla)
00489 {
00490   char *xcoord, *ycoord;
00491   int j,counter,ASCII, status;
00492 
00493   double xccd,yccd;
00494   
00495   ASCII = !plot->Graph;
00496 
00497   if (ASCII) {
00498 
00499     /* ===== Initialization (ASCII) ============================== */
00500   
00501     xcoord = (char *)malloc(10*(NCOORD+1)*sizeof(char));
00502     ycoord = (char *)malloc(10*(NCOORD+1)*sizeof(char));
00503 
00504     strcpy(xcoord,"-xdata {");
00505     strcpy(ycoord,"-ydata {"); 
00506   }
00507 #ifdef HAVE_LIBDISLIN
00508   else {
00509     hsymbl(plot->Symbol.size);
00510     setclr(plot->Symbol.color);
00511   }
00512 #endif
00513 
00514   /* Compute the image position on MLA if needed */
00515   if (*xmla == MAXDOUBLE || *ymla == MAXDOUBLE)
00516     /* [CCD/px] -> [MLA/mm] */
00517     snifs_optics_CCD2MLA(optics, xppup,yppup, xmla,ymla, pixsize,lbdaref);
00518 
00519   /* Loop over wavelengthes */
00520   for (counter=j=0; j<nlbda; j++) {
00521     /* [MLA/mm] -> [CCD/px] */
00522     snifs_optics_MLA2CCD(optics, *xmla,*ymla, &xccd,&yccd, 
00523                          pixsize,lbda[j],order);
00524 
00525     if (dxnc && order==1) /* Apply local adjustment to mask position */
00526       xccd += val_poly_nag(lbda[j], dxcoeff, dxnc,
00527                            optics->filter.inf_util, optics->filter.sup_util,
00528                            &status);
00529     if (shift != NULL) {  /* Apply offset */
00530       xccd -= shift[0];
00531       yccd -= shift[1];
00532     }
00533 
00534     if (inside_plot_limits(xccd,yccd,plot->Limits)) {
00535 
00536       if (ASCII) {
00537 
00538         /* ----- Plain ASCII output ------------------------------ */
00539 
00540         sprintf(xcoord,"%s %g",xcoord,xccd);
00541         sprintf(ycoord,"%s %g",ycoord,yccd);
00542         counter ++;
00543               
00544         if (counter%NCOORD == 0) { /* Display & reinitialize strings */
00545           printf("%s} %s}\n",xcoord,ycoord);
00546           fflush(stdout);
00547 
00548           strcpy(xcoord,"-xdata {");
00549           strcpy(ycoord,"-ydata {");
00550         }
00551       }
00552 #ifdef HAVE_LIBDISLIN
00553       else {
00554         char label[256];
00555 
00556         /* ----- Graphical output ------------------------------ */
00557 
00558         if (order!=0 && !counter) {          /* Label w/ lens_number[order] */
00559           sprintf(label," %d[%d]",nopup,order);
00560           rlmess(label,xccd,yccd);
00561           counter++;
00562         }
00563 
00564         if (plot->Symbol.size == 20) { /* Hard-coded: see plot_optics */
00565           if (order == 1)         
00566             hsymbl((int)plot->Symbol.size*
00567                    blaze_function(&(optics->grism),order,lbda[j]));
00568           else
00569             hsymbl((int)plot->Symbol.size*10*
00570                    blaze_function(&(optics->grism),order,lbda[j]));
00571         }
00572 
00573         rlsymb(plot->Symbol.shape,xccd,yccd);
00574       }
00575 #endif
00576     } /* inside_plot_limits */
00577 
00578   } /* lbda */
00579 
00580   if (ASCII) {
00581 
00582     /* ===== Conclusion (ASCII) ============================== */
00583 
00584     if (counter%NCOORD != 0) { /* Display strings */
00585       printf("%s} %s}\n",xcoord,ycoord);
00586       fflush(stdout);
00587     }
00588 
00589     free(xcoord); free(ycoord);
00590   }
00591 }
00592 
00593 /* === Doxygen Comment ======================================= */
00616 /* =========================================================== */
00617 
00618 void plot_lens_tab(const SnifsOptics *optics, float xppup,float yppup,
00619                    int nlbda,const float lbda[], int order,
00620                    TABLE *table, const int colid[], int row,
00621                    float lbdaref,float pixsize,
00622                    int dxnc, double *dxcoeff, const double shift[],
00623                    double *xmla,double *ymla)
00624 {
00625   int j, status;
00626   float xf, yf;
00627   double xccd,yccd;
00628 
00629   /* Compute the image position on MLA if needed */
00630   if (*xmla == MAXDOUBLE || *ymla == MAXDOUBLE)
00631     /* [CCD/px] -> [MLA/mm] */
00632     snifs_optics_CCD2MLA(optics, xppup,yppup, xmla,ymla, pixsize,lbdaref);
00633 
00634   /* Loop over wavelengthes */
00635   for (j=0; j<nlbda; j++) {
00636     /* [MLA/mm] -> [CCD/px] */
00637     snifs_optics_MLA2CCD(optics, *xmla,*ymla, &xccd,&yccd, 
00638                          pixsize,lbda[j],order);
00639 
00640     if (dxnc && order==1) /* Apply local adjustment to mask position */
00641       xccd += val_poly_nag(lbda[j], dxcoeff, dxnc,
00642                            optics->filter.inf_util, optics->filter.sup_util,
00643                            &status);
00644     if (shift != NULL) {  /* Apply offset */
00645       xccd -= shift[0];
00646       yccd -= shift[1];
00647     }
00648 
00649     /* Storage in table */
00650     xf = xccd; yf = yccd;
00651     WR_tbl(table,row,colid[2*j],  &xf);
00652     WR_tbl(table,row,colid[2*j+1],&yf);
00653 
00654   } /* lbda */
00655 
00656 }
00657 
00658 /* === Doxygen Comment ======================================= */
00669 /* =========================================================== */
00670 
00671 void plot_array_err(const PlotSymbol *symbol, float x[], float y[], int n,
00672                     float dx[], float dy[])
00673 {
00674 #ifdef HAVE_LIBDISLIN
00675   incmrk(-1); /* Plot every symbols without curve */
00676   marker(symbol->shape); /* Symbol */
00677   hsymbl(symbol->size);  /* Size   */
00678   setclr(symbol->color); /* Color  */
00679 
00680   /* Symbols */
00681   curve(x,y,n);
00682   
00683   /* Error-bars (if any) */
00684   if (dx != NULL) {
00685     bartyp("HORI");
00686     errbar(x,y,dx,dx,n);
00687   }
00688   if (dy != NULL) {
00689     bartyp("VERT");
00690     errbar(x,y,dy,dy,n);
00691   }
00692 
00693   incmrk(0); /* Default */
00694 
00695 #else
00696   print_error("Function plot_array_err requires DISLIN support");
00697   exit_session(UNKNOWN);  
00698 #endif
00699 }
00700 
00701 /* === Doxygen Comment ======================================= */
00713 /* =========================================================== */
00714 
00715 void plot_initialize(Plot *plot, int GRAPH, const char *devname, 
00716                      const PlotLimits *limits, const PlotSymbol *symbol, 
00717                      float mincut, float maxcut)
00718 {
00719   plot->Graph = GRAPH;
00720   strcpy(plot->devname,devname);
00721 
00722   if (limits != NULL) {
00723     plot->Limits.xmin = limits->xmin;
00724     plot->Limits.xmax = limits->xmax;
00725     plot->Limits.ymin = limits->ymin;
00726     plot->Limits.ymax = limits->ymax;
00727   }
00728   else {
00729     plot->Limits.xmin = 0;
00730     plot->Limits.xmax = 0;
00731     plot->Limits.ymin = 0;
00732     plot->Limits.ymax = 0;
00733   }
00734   
00735   if (symbol != NULL) {
00736     plot->Symbol.shape = symbol->shape;
00737     plot->Symbol.size  = symbol->size;
00738     plot->Symbol.color = symbol->color;
00739   }
00740   else {
00741     /* Default symbol */
00742     plot->Symbol.shape = 16;  /* 0=[], 3=+, 5=<>, 15=(), 16=[*], 21=(*) */
00743     plot->Symbol.size  = 15;  /* DISLIN default: 35 */
00744     plot->Symbol.color = 255; /* Foreground */
00745   }
00746 
00747   plot->mincut = mincut;
00748   plot->maxcut = maxcut;
00749 }
00750 
00751 /* 
00752 ;;; Local Variables: ***
00753 ;;; eval: (add-to-list 'c-font-lock-extra-types "PlotLimits") ***
00754 ;;; eval: (add-to-list 'c-font-lock-extra-types "PlotSymbol") ***
00755 ;;; eval: (add-to-list 'c-font-lock-extra-types "Plot") ***
00756 ;;; End: ***
00757 */

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