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 #define 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;
00041   *emsig = sqrt( (SQ(parpeak[3]*eparpeak[3]) + SQ(parpeak[4]*eparpeak[4])) /
00042                  SQ(*msig) )/2;            /* Mean quadratic error */
00043 }
00044 
00045 /* === Doxygen Comment ======================================= */
00049 /* =========================================================== */
00050 
00051 void print_parpeak(double parpeak[6], double eparpeak[6], 
00052                    double msig, double emsig)
00053 {
00054   print_msg("   Relative offset: %.3g x %.3g px (+- %.3g x %.3g)", 
00055             parpeak[0],parpeak[1],eparpeak[0],eparpeak[1]);
00056   print_msg("   Intensity: %.3g (+- %.3g)", 
00057             parpeak[2],eparpeak[2]);
00058   print_msg("   Background: %.3g (+- %.3g)", parpeak[5],eparpeak[5]);
00059   print_msg("   Dispersion: %.2f x %.2f px (+- %.2f x %.2f)", 
00060             parpeak[3],parpeak[4],eparpeak[3],eparpeak[4]);
00061   print_msg("   -> Mean Dispersion: %.3g +- %.3g px",msig,emsig);
00062 }
00063 
00064 /* === Doxygen Comment ======================================= */
00084 /* =========================================================== */
00085 
00086 int comp_calibration_offsets(char *refname, float *ref_airmass, 
00087                              char *arcname, float *cal_airmass, float *cal_jdate,
00088                              char *arc2name, float *cal2_airmass, float *cal2_jdate,
00089                              double offset[2])
00090 {
00091   Channel channel, channel2;
00092   IMAGE2D refframe, arcframe;
00093 
00094   int status;
00095 
00096   double **ref_tab, **cal_tab;
00097   double parpeak[6], eparpeak[6], offset2[2];
00098   double msig, emsig;
00099 
00100   /* ----- Reference calibration frame ------------------------------ */
00101 
00102   /* Open frame & Check file class */
00103   if ((status = open_image(&refframe,refname,"I",PRE_CAL_FRAME,"reference arc")))
00104     return(status);
00105 
00106   /* Don't stop if input ref. arc has a missing CHANNEL */
00107   if (!(channel = read_channel(&refframe))) 
00108     print_warning("Cannot test channel coherence among calibration frames");
00109 
00110   *ref_airmass = read_airmass(&refframe); /* Read AIRMASS descriptor */
00111 
00112   /* Extract central part */
00113   alloc2d(&ref_tab,ARCWIDTH,ARCWIDTH,DOUBLE);
00114   if (extract_arc(&refframe, ARCWIDTH, ref_tab)) 
00115     return(UNKNOWN);
00116 
00117   close_frame(&refframe);
00118     
00119   /* ----- 1st calibration frame ------------------------------ */
00120 
00121   print_msg("Computing relative offset between reference (%s) and "
00122             "calibration (%s) frames...",refname,arcname);
00123 
00124   /* Open frame & Check file class */
00125   if ((status = open_image(&arcframe,arcname,"I",PRE_CAL_FRAME,"arc")))
00126     return(status);
00127 
00128   /* Check channel */
00129   if (channel && (channel2 = read_channel(&arcframe))!=channel) {
00130     print_error("reference channel (%s) and calibration channel (%s) are different",
00131                 channel_names[channel],channel_names[channel2]);
00132     return(ERR_BAD_PARAM);
00133   }
00134 
00135   *cal_airmass = read_airmass(&arcframe);
00136   *cal_jdate   = read_juliandate(&arcframe);
00137 
00138   /* Extract central part */
00139   alloc2d(&cal_tab,ARCWIDTH,ARCWIDTH,DOUBLE);
00140   if (extract_arc(&arcframe, ARCWIDTH, cal_tab)) 
00141     return(UNKNOWN);
00142 
00143   close_frame(&arcframe);
00144 
00145   /* Fit the X-corr. peak */
00146   if (Fit_Xpeak(ref_tab,cal_tab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,parpeak,eparpeak))
00147     return(UNKNOWN);
00148 
00149   if (DEBUG) {
00150     /* Convert correlation quantities to intrinsic quantities */
00151     comp_mean_disp(parpeak, eparpeak, &msig, &emsig);
00152     print_parpeak(parpeak, eparpeak, msig, emsig);
00153   }
00154     
00155   offset[0] = parpeak[0];
00156   offset[1] = parpeak[1];
00157   print_msg("   Relative offset: %.2f x %.2f px", offset[0],offset[1]);
00158 
00159   /* ----- 2nd calibration frame (if any) ------------------------------ */
00160 
00161   if (arc2name != NULL) {
00162 
00163     print_msg("Computing relative offset between reference (%s) and "
00164               "2nd calibration (%s) frames...",refname,arc2name);
00165 
00166     /* Open frame & Check file class */
00167     if ((status = open_image(&arcframe,arc2name,"I",PRE_CAL_FRAME,"arc (II)"))) 
00168       return(status);
00169 
00170     /* Check channel */
00171     if ((channel2 = read_channel(&arcframe)) != channel) {
00172       print_error("reference channel (%s) and calibration channel (%s) are different",
00173                   channel_names[channel],channel_names[channel2]);
00174       return(ERR_BAD_PARAM);
00175     }
00176 
00177     *cal2_airmass = read_airmass(&arcframe);
00178     *cal2_jdate   = read_juliandate(&arcframe);
00179     
00180     /* Extract central part */
00181     if (extract_arc(&arcframe, ARCWIDTH, cal_tab)) 
00182       return(UNKNOWN);
00183 
00184     close_frame(&arcframe);
00185 
00186     /* Fit the X-corr. peak */
00187     if (Fit_Xpeak(ref_tab,cal_tab,ARCWIDTH,ARCWIDTH,10,ARCWIDTH/4,parpeak,eparpeak))
00188       return(UNKNOWN);
00189 
00190     if (DEBUG) {
00191       /* Convert correlation quantities to intrinsic quantities */
00192       comp_mean_disp(parpeak, eparpeak, &msig, &emsig);
00193       print_parpeak(parpeak, eparpeak, msig, emsig);
00194     }
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(&cal_tab,DOUBLE); free2d(&ref_tab,DOUBLE);
00206 
00207   return(OK);
00208 }

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