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

extract/lib/proc_psf.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_sf_erf.h>
00020 
00021 /* === Doxygen Comment ======================================= */
00035 /* =========================================================== */
00036 
00037 int read_PSFgeo_param(void *anyfile, double *geo_param[]) 
00038 {
00039   int nggeo;
00040   
00041   if (RD_desc(anyfile,MPUP_NGGEO,INT,1,&nggeo) < 0) {
00042     print_error("Cannot read descriptor %s",MPUP_NGGEO);
00043     nggeo = -1;
00044   }
00045   else {
00046     *geo_param = (double *)malloc(3*nggeo*sizeof(double));
00047     if (RD_desc(anyfile,MPUP_GEOPAR,DOUBLE,3*nggeo,*geo_param) < 0) {
00048       print_error("Cannot read descriptor %s",MPUP_GEOPAR);
00049       free(*geo_param);
00050       nggeo = -1;
00051     }
00052     else
00053       print_msg("   Geometric PSF model: %d component(s)",nggeo);
00054   }
00055   return(nggeo);
00056 }
00057 
00058 /* === Doxygen Comment ======================================= */
00071 /* =========================================================== */
00072 
00073 int read_PSFglobal_param(void *anyfile, double *psf_param[]) 
00074 {
00075   int ngpsf;
00076   
00077   if (RD_desc(anyfile,MPUP_NGPSF,INT,1,&ngpsf) < 0) {
00078     print_error("Cannot read descriptor %s",MPUP_NGPSF);
00079     ngpsf = -1;
00080   }
00081   else {
00082     *psf_param = (double *)malloc((1+2*ngpsf)*sizeof(double));
00083     if (RD_desc(anyfile,MPUP_PSFPAR,DOUBLE,(1+2*ngpsf),*psf_param) < 0) {
00084       print_error("Cannot read descriptor %s",MPUP_PSFPAR);
00085       free(*psf_param);
00086       ngpsf = -1;
00087     }
00088     else {
00089       if (DEBUG) 
00090         print_msg("   Geometric pupil scale: %8.4f px",(*psf_param)[0]);
00091       print_msg("   Global PSF model: %d component(s)",ngpsf);
00092     }
00093   }
00094   return(ngpsf);
00095 }
00096 
00097 /* === Doxygen Comment ======================================= */
00106 /* =========================================================== */
00107 
00108 int read_PSFlocal_ng(void *anyfile) 
00109 {
00110   int ngloc;
00111   
00112   if (RD_desc(anyfile,MPUP_NGLOC,INT,1,&ngloc) < 0) {
00113     print_error("Cannot read descriptor %s",MPUP_NGLOC);
00114     ngloc = MAXINT;
00115   }
00116   else
00117     print_msg("   Local PSF model: %d component(s)",ABS(ngloc));
00118 
00119   return(ngloc);
00120 }
00121 
00122 /* === Doxygen Comment ======================================= */
00134 /* =========================================================== */
00135 
00136 int read_corePSF(void *anyfile, CorePSF *psf) 
00137 {
00138   int i,j;
00139   double *geo_param, *psf_param;
00140 
00141   /* ..... Geometric PSF .............................. */
00142 
00143   if ((psf->nggeo = read_PSFgeo_param(anyfile, &geo_param)) < 0)
00144     return(ERR_NODESC);
00145 
00146   psf->Igeo = (double *)malloc(psf->nggeo*sizeof(double));
00147   psf->Xgeo = (double *)malloc(psf->nggeo*sizeof(double));
00148   psf->Sgeo = (double *)malloc(psf->nggeo*sizeof(double));
00149   for (i=0; i<psf->nggeo; i++) {
00150     psf->Igeo[i] = geo_param[3*i  ];
00151     psf->Xgeo[i] = geo_param[3*i+1];
00152     psf->Sgeo[i] = geo_param[3*i+2];
00153     if (DEBUG) print_msg("\tGaussian %d: I = %f, X = %f, S = %f",i+1,
00154                          geo_param[3*i],geo_param[3*i+1],geo_param[3*i+2]);
00155   }
00156   free(geo_param);
00157 
00158   /* ..... Global PSF .............................. */
00159 
00160   if ((psf->ngpsf = read_PSFglobal_param(anyfile, &psf_param)) < 0)
00161     return(ERR_NODESC);
00162 
00163   psf->scale = psf_param[0];
00164 
00165   if (DEBUG) 
00166     for (i=0; i<psf->ngpsf; i++) 
00167       print_msg("\tUnconvolved gaussian %d: I = %f, S = %f", 
00168                 i+1,psf_param[2*i+1],psf_param[2*i+2]);
00169 
00170   /* ..... Local PSF .............................. */
00171 
00172   if ((psf->ngloc = read_PSFlocal_ng(anyfile)) == MAXINT)
00173     return(ERR_NODESC);
00174 
00175   /* ..... Core PSF .............................. */
00176 
00177   if (psf->ngloc > 0) {
00178     print_msg("   Core PSF = Geometric PSF (x) Global PSF");
00179   }
00180   else {
00181     print_msg("   Core PSF = Geometric PSF (global PSF *discarded*)");
00182     psf->ngloc = -psf->ngloc;
00183     psf->ngpsf = 0;
00184   }
00185 
00186   if (psf->ngloc > 1) {
00187     print_error("Multi-gaussian local PSF model not implemented yet");
00188     return(ERR_NOIMPL);
00189   }
00190 
00191   alloc2d(&(psf->Ipsf),psf->nggeo,MAX(psf->ngpsf,1),DOUBLE);
00192   alloc2d(&(psf->Spsf),psf->nggeo,MAX(psf->ngpsf,1),DOUBLE);
00193 
00194   psf->norm = 0;
00195   
00196   if (psf->ngpsf > 0) {
00197     /* Use global PSF derived previously: core PSF = geom PSF (x) glob PSF */
00198     for (i=0; i<psf->nggeo; i++) {   /* Loop on geom PSF components   */
00199       for (j=0; j<psf->ngpsf; j++) { /* Loop on global PSF components */
00200         psf->Spsf[i][j] = hypot(psf_param[2*j+2],psf->Sgeo[i]);
00201         psf->Ipsf[i][j] = M_SQRT2PI*psf_param[2*j+1]*psf->Igeo[i] * 
00202           psf_param[2*j+2]*psf->Sgeo[i] / psf->Spsf[i][j];
00203         psf->norm += psf->Ipsf[i][j]*M_SQRT2PI*psf->Spsf[i][j];
00204       }
00205     }
00206   }
00207   else {
00208     /* Do *not* use global PSF derived previously: core PSF = geom PSF */
00209     print_warning("Global PSF discarded: core PSF = geom PSF");
00210     for (i=0; i<psf->nggeo; i++) { /* Loop on geom PSF components */
00211       psf->Spsf[i][0] = psf->Sgeo[i];
00212       psf->Ipsf[i][0] = psf->Igeo[i];
00213       psf->norm += psf->Ipsf[i][0]*M_SQRT2PI*psf->Spsf[i][0];
00214     }
00215   }
00216   free(psf_param);
00217 
00218   return(OK);
00219 }
00220 
00221 /* === Doxygen Comment ======================================= */
00225 /* =========================================================== */
00226 
00227 void free_corePSF(CorePSF *psf) 
00228 {
00229   free(psf->Igeo); free(psf->Xgeo); free(psf->Sgeo);
00230   free2d(&(psf->Ipsf),DOUBLE); free2d(&(psf->Spsf),DOUBLE);
00231 
00232   psf->nggeo = psf->ngpsf = psf->ngloc = 0;
00233 }
00234 
00235 /* === Doxygen Comment ======================================= */
00246 /* =========================================================== */
00247 
00248 double xdisp_profile(CorePSF *psf, SpectrumCCD *spectrum, float x) 
00249 {
00250   double val, qe, qs, xs, xe;
00251   int l, p;
00252 
00253   xs = x - 0.5;
00254   xe = x + 0.5;
00255   for (val=0., l=0; l<psf->nggeo; l++) {     /* Loop on geom pupil */
00256     for (p=0; p<psf->ngpsf; p++) {           /* Loop on PSF profile */
00257       qs = (xs - psf->Xgeo[l])/spectrum->Scs[l][p]/M_SQRT2;
00258       qe = (xe - psf->Xgeo[l])/spectrum->Scs[l][p]/M_SQRT2;
00259       val += spectrum->Scs[l][p]*spectrum->Ics[l][p] * 
00260         ( gsl_sf_erf(qe) - gsl_sf_erf(qs) );
00261     }
00262   }
00263   return(val);
00264 }
00265 

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