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

extract/lib/proc_extract.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 
00020 #include <gsl/gsl_sf_erf.h>
00021 #include <gsl/gsl_errno.h>
00022 #include <gsl/gsl_spline.h>
00023 
00024 /* === Doxygen Comment ======================================= */
00034 /* =========================================================== */
00035 
00036 int good_neighbor(float *ylo, float *yup, int *neighbors, float y)
00037 {
00038   int up, lo;
00039   float trans;
00040 
00041   /* which neighbor is higher? */
00042   if ( (ylo[neighbors[0]]+yup[neighbors[0]]) > 
00043        (ylo[neighbors[1]]+yup[neighbors[1]]) ) { /* 0 is higher */
00044     up = neighbors[0];
00045     lo = neighbors[1];
00046   } 
00047   else {                                     /* 1 is higher */
00048     up = neighbors[1];
00049     lo = neighbors[0];
00050   }
00051   /* compute the transition in the dark zone between the two neighbors */
00052   trans = (yup[lo] + ylo[up])/2.;
00053   /* so, which one is the good one? */
00054   if (y < trans) return(lo);
00055   else           return(up);
00056 }
00057 
00058 /* === Doxygen Comment ======================================= */
00069 /* =========================================================== */
00070 
00071 int good_neighbor_spectrum(SpectrumCCD spectra[], int neighbors[], float y)
00072 {
00073   int up, lo;
00074   float trans;
00075 
00076   /* which neighbor is higher? */
00077   if ( (spectra[neighbors[0]].total.ylo + spectra[neighbors[0]].total.yup) > 
00078        (spectra[neighbors[1]].total.ylo + spectra[neighbors[1]].total.yup) ) { 
00079     /* 0 is higher */
00080     up = neighbors[0];
00081     lo = neighbors[1];
00082   } 
00083   else {
00084     /* 1 is higher */
00085     up = neighbors[1];
00086     lo = neighbors[0];
00087   }
00088   /* compute the transition in the dark zone between the two neighbors */
00089   trans = (spectra[lo].total.yup + spectra[up].total.ylo)/2.;
00090   /* so, which one is the good one? */
00091   if (y < trans) return(lo);
00092   else           return(up);
00093 }
00094 
00095 /* === Doxygen Comment ======================================= */
00119 /* =========================================================== */
00120 
00121 void interpolate_y2lbda(int n, double y[], double signal[], double variance[],
00122                         int yofl_deg, float yofl_coeffs[],
00123                         SPECTRUM *spectrum, SPECTRUM *noise)
00124 {
00125   gsl_interp_accel *acc, *acc2;
00126   gsl_spline *spline, *spline2;
00127   int i, NOISE;
00128   double yl;
00129   double *ycopy, *scopy, *vcopy;
00130   
00131   /* NOISE if both variance and output variance spectrum are set */
00132   NOISE = (variance != NULL && noise != NULL);
00133   
00134   /* Copy input arrays so that they're ordered increasingly according to y */
00135   ycopy = (double *)malloc(n*sizeof(double));
00136   scopy = (double *)malloc(n*sizeof(double));
00137   if (NOISE) vcopy = (double *)malloc(n*sizeof(double));
00138   for (i=0; i<n; i++) {
00139     ycopy[i] = y[n-1-i];
00140     scopy[i] = signal[n-1-i];
00141     if (NOISE) vcopy[i] = variance[n-1-i];
00142   }
00143   
00144   acc = gsl_interp_accel_alloc();
00145   /* Use linear interpolation (could be 'cspline' or 'akima') */
00146   spline = gsl_spline_alloc(gsl_interp_linear, n);
00147   gsl_spline_init(spline, ycopy, scopy, n);
00148 
00149   if (NOISE) {
00150     acc2 = gsl_interp_accel_alloc();
00151     /* Use linear interpolation (could be 'cspline' or 'akima') */
00152     spline2 = gsl_spline_alloc(gsl_interp_linear, n);
00153     gsl_spline_init(spline2, ycopy, vcopy, n);
00154   }
00155 
00156   /* Interpolate flux on wavelength ramp */
00157   for (i=0; i<spectrum->npts; i++) {
00158     yl = Polynom(yofl_deg,yofl_coeffs,       /* y(lbda) */
00159                  coord_spec(spectrum,i));    /* lbda */
00160     WR_spec(spectrum,i, gsl_spline_eval(spline, yl, acc));
00161     if (NOISE)
00162       WR_spec(noise,i, gsl_spline_eval(spline2, yl, acc2));
00163   }
00164 
00165   gsl_spline_free(spline);
00166   gsl_interp_accel_free(acc);
00167   free(ycopy);
00168   free(scopy);
00169   if (NOISE) {
00170     gsl_spline_free(spline2);
00171     gsl_interp_accel_free(acc2);
00172     free(vcopy);
00173   }
00174 }

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