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

extract/lib/proc_shift.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 static int const ARCWIDTH=128; 
00021 
00022 /* === Doxygen Comment ======================================= */
00032 /* =========================================================== */
00033 
00034 void comp_mean_disp(double parpeak[6], double eparpeak[6], 
00035                     double *msig, double *emsig) 
00036 {
00037   parpeak[3] /= M_SQRT2; eparpeak[3] /= M_SQRT2;
00038   parpeak[4] /= M_SQRT2; eparpeak[4] /= M_SQRT2;
00039 
00040   *msig  = hypot(parpeak[3],parpeak[4])/M_SQRT2; /* Quadratic mean */
00041   *emsig = sqrt( (SQ(parpeak[3]*eparpeak[3]) + 
00042                   SQ(parpeak[4]*eparpeak[4]))/2 ) / *msig;
00043 }
00044 
00045 /* === Doxygen Comment ======================================= */
00051 /* =========================================================== */
00052 
00053 void print_parpeak(double parpeak[6], double eparpeak[6])
00054 {
00055   double msig, emsig;
00056   
00057   print_msg("   Relative offset: %.3g x %.3g px (+- %.3g x %.3g)", 
00058             parpeak[0],parpeak[1],eparpeak[0],eparpeak[1]);
00059   print_msg("   Intensity: %.3g (+- %.3g)", 
00060             parpeak[2],eparpeak[2]);
00061   print_msg("   Background: %.3g (+- %.3g)", parpeak[5],eparpeak[5]);
00062   print_msg("   Dispersion: %.2f x %.2f px (+- %.2f x %.2f)", 
00063             parpeak[3],parpeak[4],eparpeak[3],eparpeak[4]);
00064   comp_mean_disp(parpeak, eparpeak, &msig, &emsig);
00065   print_msg("   -> Mean Dispersion: %.3g +- %.3g px",msig,emsig);
00066 }
00067 
00068 /* === Doxygen Comment ======================================= */
00087 /* =========================================================== */
00088 
00089 int comp_calibration_offsets(char *refName, 
00090                              char *arcName, float *arcAirmass, float *arcJDate,
00091                              char *arc2Name,float *arc2Airmass,float *arc2JDate,
00092                              float offset[2])
00093 {
00094   Channel channel, channel2;
00095   IMAGE2D refFrame, arcFrame;
00096 
00097   int status;
00098 
00099   float offset2[2];
00100 
00101   double **refTab, **calTab;
00102   double parpeak[6], eparpeak[6];
00103 
00104   /* ----- Reference calibration frame ------------------------------ */
00105 
00106   /* Open frame & Check file class */
00107   if ((status = open_image(&refFrame,refName,"I",
00108                            PRE_CAL_FRAME,"reference arc")))
00109     return(status);
00110 
00111   /* Don't stop if input ref. arc has a missing CHANNEL */
00112   if (!(channel = read_channel(&refFrame))) 
00113     print_warning("Cannot test channel coherence among calibration frames");
00114 
00115   /* Extract central part */
00116   alloc2d(&refTab,ARCWIDTH,ARCWIDTH,DOUBLE);
00117   if (extract_arc(&refFrame, ARCWIDTH, refTab)) 
00118     return(UNKNOWN);
00119 
00120   close_frame(&refFrame);
00121     
00122   /* ----- 1st calibration frame ------------------------------ */
00123 
00124   print_msg("Computing relative offset between reference (%s) and "
00125             "calibration (%s) frames...",refName,arcName);
00126 
00127   /* Open frame & Check file class */
00128   if ((status = open_image(&arcFrame,arcName,"I",PRE_CAL_FRAME,"arc")))
00129     return(status);
00130 
00131   /* Check channel */
00132   if (channel && (channel2 = read_channel(&arcFrame))!=channel) {
00133     print_error("reference channel (%s) and calibration channel (%s) "
00134                 "are different",channel_names[channel],channel_names[channel2]);
00135     return(ERR_BAD_PARAM);
00136   }
00137 
00138   *arcAirmass = read_airmass(&arcFrame);
00139   *arcJDate   = read_jdate(&arcFrame);
00140 
00141   /* Extract central part */
00142   alloc2d(&calTab,ARCWIDTH,ARCWIDTH,DOUBLE);
00143   if (extract_arc(&arcFrame, ARCWIDTH, calTab)) 
00144     return(UNKNOWN);
00145 
00146   close_frame(&arcFrame);
00147 
00148   /* Fit the X-corr. peak */
00149   if (Fit_Xpeak(refTab,calTab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,
00150                 parpeak,eparpeak))
00151     return(UNKNOWN);
00152 
00153   if (DEBUG) /* Convert correlation quantities to intrinsic quantities */
00154     print_parpeak(parpeak, eparpeak);
00155 
00156   offset[0] = parpeak[0];
00157   offset[1] = parpeak[1];
00158   print_msg("   Relative offset: %.2f x %.2f px", offset[0],offset[1]);
00159 
00160   /* ----- 2nd calibration frame (if any) ------------------------------ */
00161 
00162   if (arc2Name != NULL) {
00163 
00164     print_msg("Computing relative offset between reference (%s) and "
00165               "2nd calibration (%s) frames...",refName,arc2Name);
00166 
00167     /* Open frame & Check file class */
00168     if ((status = open_image(&arcFrame,arc2Name,"I",PRE_CAL_FRAME,"arc (II)"))) 
00169       return(status);
00170 
00171     /* Check channel */
00172     if (channel && (channel2 = read_channel(&arcFrame))!=channel) {
00173       print_error("reference channel (%s) and calibration channel (%s) "
00174                   "are different",
00175                   channel_names[channel],channel_names[channel2]);
00176       return(ERR_BAD_PARAM);
00177     }
00178 
00179     *arc2Airmass = read_airmass(&arcFrame);
00180     *arc2JDate   = read_jdate(&arcFrame);
00181     
00182     /* Extract central part */
00183     if (extract_arc(&arcFrame, ARCWIDTH, calTab)) 
00184       return(UNKNOWN);
00185 
00186     close_frame(&arcFrame);
00187 
00188     /* Fit the X-corr. peak */
00189     if (Fit_Xpeak(refTab,calTab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,
00190                   parpeak,eparpeak))
00191       return(UNKNOWN);
00192 
00193     if (DEBUG) /* Convert correlation quantities to intrinsic quantities */
00194       print_parpeak(parpeak, eparpeak);
00195     
00196     offset2[0] = parpeak[0];
00197     offset2[1] = parpeak[1];
00198     print_msg("   Relative offset: %.2f x %.2f px", offset2[0],offset2[1]);
00199 
00200     /* Average offset */
00201     offset[0] = (offset[0]+offset2[0])/2.;
00202     offset[1] = (offset[1]+offset2[1])/2.;
00203   } 
00204 
00205   free2d(&calTab,DOUBLE); free2d(&refTab,DOUBLE);
00206 
00207   return(OK);
00208 }
00209 
00210 int comp_specres(char *arcName, float *specres)
00211 {
00212   IMAGE2D arcFrame;
00213 
00214   int status;
00215 
00216   double **calTab;
00217   double parpeak[6], eparpeak[6];
00218 
00219   /* Open frame & Check file class */
00220   if ((status = open_image(&arcFrame,arcName,"I",PRE_CAL_FRAME,"arc")))
00221     return(status);
00222 
00223   /* Extract central part */
00224   alloc2d(&calTab,ARCWIDTH,ARCWIDTH,DOUBLE);
00225   if (extract_arc(&arcFrame, ARCWIDTH, calTab)) 
00226     return(UNKNOWN);
00227 
00228   close_frame(&arcFrame);
00229 
00230   /* Fit the X-corr. peak */
00231   if (Fit_Xpeak(calTab,calTab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,
00232                 parpeak,eparpeak))
00233     return(UNKNOWN);
00234 
00235   if (DEBUG) /* Convert correlation quantities to intrinsic quantities */
00236     print_parpeak(parpeak, eparpeak);
00237     
00238   *specres = parpeak[4]/M_SQRT2; /* Spectral resolution [px] */
00239 
00240   if (DEBUG)
00241     print_msg("   Spectral resolution: %.3f px", specres);
00242 
00243   free2d(&calTab,DOUBLE);
00244 
00245   return(OK);
00246 }

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