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

extract/lib/fit_mpup.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 #ifdef HAVE_LIBDISLIN
00022 #include <dislin.h>
00023 #endif
00024 
00026 #define DISLIN 1
00027 
00028 /* Global variables needed by \c nllsqfit_distMLA and \c model_distMLA */
00029 int *glmlaml, *glmlanl;                      /* Mpup position indices */      
00030 double *glmlaxpup, *glmlaypup;               /* Mpup positions */
00031 double glmlaxlc, glmlaylc;                   /* Central lens position */
00032 
00033 /* Global variables describing the geometric PSF */
00034 extern int glnggeo;                          /* Nb of gaussian components */
00035 extern double *glIgeo, *glXgeo, *glSgeo;
00036 
00037 /* Global variables describing the global PSF */
00038 extern int glngpsf;                          /* Nb of gaussian components */
00039 
00040 /* Global variables describing the core PSF = geom. PSF (x) global PSF */
00041 extern double **glIpsf, **glSpsf; /* Global var. needed by nllsqfit_globalPSF */
00042 
00043 /* Global variables describing the local PSF */
00044 extern int glngloc;                          /* Nb of gaussian components */
00045 
00046 /* Global variables needed by nllsqfit_globalPSF */
00047 int *glpcen;               /* Length of the central mpup profiles */      
00048 
00049 /* === Doxygen Comment ======================================= */
00063 /* =========================================================== */
00064 
00065 double gaussian_mpup(double arr[], int npts, double start, double step, 
00066                      double par[], double dpar[], double fit[])
00067 {
00068   int i,status;
00069   double chi2, gof, rms;
00070   
00071   /* Initial guess */
00072   par[1] = npts/2.;            /* The mpup is already roughly centered */
00073   par[2] = 1;                          /* Anything sensible will do it */
00074   par[0] = arr[((int)par[1])]*par[2]*M_SQRT2PI; /* Normalized gaussian */
00075 
00076   /* Simple least-sq. fit */
00077   status = nllsqfit_bnd(NULL, arr, NULL, npts, par, NULL, NULL, 3,
00078                         *nllsqfit_gaussian, fit, NULL, dpar, &chi2,&gof);
00079 
00080   /* Convert parameters from pixel to world coordinates */
00081   if (!status) {
00082     par[1]   = start + par[1]*step; dpar[1] *= step; /* mean */
00083     par[2]  *= step;                dpar[2] *= step; /* disp */
00084     rms = sqrt(chi2/npts);
00085   }
00086   else {
00087     for (i=0; i<3; i++) par[i] = 0.;
00088     rms = -status;
00089     for (i=0; i<npts; i++) fit[i] = 0.;      /* Zero-ing the fit profile */
00090   }
00091 
00092   return(rms);
00093 }
00094 
00095 /* ##### MLA distortion ############################## */
00096 
00097 /* === Doxygen Comment ======================================= */
00110 /* =========================================================== */
00111 
00112 long nllsqfit_distMLA(long *mode, long *npup, long *npar, long *ldfj, 
00113                       double *par, double f[], double fjac[], long *nstate, 
00114                       long *iuser, double user[])
00115 {
00116   int i;
00117   double r, r2, K, K1, K2, dxd, dyd;
00118   double xd, yd, xnd, ynd, dx, dy, xrm, yrm;
00119   double coeff, tmp;
00120 
00121 #define dist  par[0]*FMP_LENSSNORM /* MLA lenssize */
00122 #define angle par[1]*FMP_ANGLENORM /* MLA angle */
00123 #define xc    par[2]*FMP_POSITNORM /* Center of distortion */
00124 #define yc    par[3]*FMP_POSITNORM
00125 #define alpha par[4]*FMP_ALPHANORM /* Distortion: alpha amplitude */
00126 #ifdef USE_R4_DIST
00127 #define beta  par[5]*SQ(FMP_ALPHANORM) /* Distortion: beta amplitude */
00128 #endif
00129 
00130   for (i=0; i<*npup ;i++) { /* Loop over mpups */
00131     /* Rotation */
00132     K1 = glmlaml[i]*cos(angle) - glmlanl[i]*sin(angle);
00133     K2 = glmlaml[i]*sin(angle) + glmlanl[i]*cos(angle);
00134     /* Undistorted position from central lens */
00135     xnd = glmlaxlc + K1*dist;
00136     ynd = glmlaylc + K2*dist;
00137     /* Radial distortion: rd = r(1 + alpha*r2 [+ beta*r4]) */
00138     dx = xnd - xc;
00139     dy = ynd - yc;
00140     r = hypot(dx,dy); /* Distance to distortion center */
00141     r2 = SQ(r);
00142     K = 1 + alpha*r2; /* Expansion factor */
00143     coeff = alpha;
00144 #ifdef USE_R4_DIST
00145     K     += beta*SQ(r2);
00146     coeff += 2*beta*r2;
00147 #endif
00148     /* Distorted position from central lens */
00149     xd = xc + K*dx;
00150     yd = yc + K*dy;
00151     /* Distance between real and theoretical distorted positions */
00152     xrm = glmlaxpup[i] - xd;
00153     yrm = glmlaypup[i] - yd;
00154 
00155     /* Subfunction */
00156     f[i] = SQ(xrm) + SQ(yrm);
00157     /* f[i] = sqrt(f[i]); */ /* shall we give f or f^2 ? */
00158 
00159     if (*mode) {
00160       /* Partial derivatives with respect to... */
00161 
00162       /* dist */
00163       tmp = 2*coeff*(K1*dx + K2*dy);
00164       dxd = K*K1 + dx*tmp;
00165       dyd = K*K2 + dy*tmp;
00166       fjac[i + *npup * 0] = -2*(xrm*dxd + yrm*dyd)*FMP_LENSSNORM;
00167       /* fjac[i + *npup * 0] /= 2*f[i]; */
00168 
00169       /* angle */
00170       tmp = 2*coeff*(K1*dy - K2*dx);
00171       dxd = -K*K2 + dx*tmp;
00172       dyd =  K*K1 + dy*tmp;
00173       fjac[i + *npup * 1] = -2*dist*(xrm*dxd + yrm*dyd)*FMP_ANGLENORM;
00174       /* fjac[i + *npup * 1] /= 2*f[i]; */
00175 
00176       /* xc */
00177       dxd = -2*coeff*SQ(dx) - K + 1;
00178       dyd = -2*coeff*dx*dy;
00179       fjac[i + *npup * 2] = -2*(xrm*dxd + yrm*dyd)*FMP_POSITNORM;
00180       /* fjac[i + *npup * 2] /= 2*f[i]; */
00181 
00182       /* yc */
00183       dxd = -2*coeff*dx*dy;
00184       dyd = -2*coeff*SQ(dy) - K + 1;
00185       fjac[i + *npup * 3] = -2*(xrm*dxd + yrm*dyd)*FMP_POSITNORM;
00186       /* fjac[i + *npup * 3] /= 2*f[i]; */
00187 
00188       /* alpha */
00189       fjac[i + *npup * 4] = -2*r2*(xrm*dx + yrm*dy)*FMP_ALPHANORM;
00190       /* fjac[i + *npup * 4] /= 2*f[i]; */
00191 
00192 #ifdef USE_R4_DIST
00193       /* beta */
00194       fjac[i + *npup * 5] = fjac[i + *npup * 4]*r2*FMP_ALPHANORM;
00195 #endif
00196     } /* Derivatives */
00197   } /* Loop over mpups */
00198   return(OK);
00199 
00200 #undef dist
00201 #undef angle
00202 #undef xc
00203 #undef yc
00204 #undef alpha
00205 #ifdef USE_R4_DIST
00206 #undef beta
00207 #endif
00208 }
00209 
00210 /* === Doxygen Comment ======================================= */
00217 /* =========================================================== */
00218 
00219 void model_distMLA(int npup, double dpar[], int npar, 
00220                    float xd[],float yd[],float rd[])
00221 {
00222   int i;
00223   double r2, K, K1, K2, xnd, ynd, dx, dy;
00224   
00225 #define dist  dpar[0] /* MLA lenssize */
00226 #define angle dpar[1] /* MLA angle */
00227 #define xc    dpar[2] /* Center of distortion */
00228 #define yc    dpar[3]
00229 #define alpha dpar[4] /* Distortion: alpha */
00230 #ifdef USE_R4_DIST
00231 #define beta  dpar[5] /* Distortion: beta */
00232 #endif
00233 
00234   for (i=0; i<npup ;i++) { /* Loop over mpups */
00235     /* Rotation */
00236     K1 = glmlaml[i]*cos(angle) - glmlanl[i]*sin(angle);
00237     K2 = glmlaml[i]*sin(angle) + glmlanl[i]*cos(angle);
00238     /* Undistorted position from central lens */
00239     xnd = glmlaxlc + K1*dist;
00240     ynd = glmlaylc + K2*dist;
00241     /* Radial distortion: rd = r(1 + alpha*r2) */
00242     dx = xnd - xc;
00243     dy = ynd - yc;
00244     r2 = SQ(dx) + SQ(dy); /* Distance to distortion center */
00245     K = 1 + alpha*r2;
00246 #ifdef USE_R4_DIST
00247     K += beta*SQ(r2);
00248 #endif
00249     xd[i] = xc + K*dx;
00250     yd[i] = yc + K*dy;
00251     /* Distance between real and theoretical positions */
00252     rd[i] = hypot(glmlaxpup[i]-xd[i],glmlaypup[i]-yd[i]);
00253   }
00254 
00255 #ifdef HAVE_LIBDISLIN
00256   if (DEBUG && !DISLIN) 
00257     print_warning("Distortion graphical output desactivated (%s:%d)",
00258                   __FILE__,__LINE__);
00259   if (DEBUG && DISLIN) {
00260     char title[4*lg_ident+1];
00261     int nlev = 10, ntri,nmax, *i1ray,*i2ray,*i3ray;
00262     float *x0,*y0, *x1,*y1, min,max, levels[nlev];
00263 
00264     Plot plot;
00265 
00266     /* Graphic output initialization */
00267     sprintf(title,"dislin_distortion.eps");
00268     print_msg("DISLIN: Output graphics in %s",title);
00269     setfil(title);                           /* output filename */
00270 
00271     plot_initialize(&plot, 1, "EPS", NULL, &(plot.Symbol), 0,0);
00272     plot_set_limits(xd,yd,npup,&(plot.Limits));
00273     plot_increase_limits(&(plot.Limits));
00274     plot_start(&plot,0);
00275     plot_axes(&plot,"x [px]","y [px]");
00276 
00277     /* Mpup index */
00278     x0 = (float *)malloc((npup+3)*sizeof(float)); /* needed by triang */
00279     y0 = (float *)malloc((npup+3)*sizeof(float));
00280     x1 = (float *)malloc(    npup*sizeof(float));
00281     y1 = (float *)malloc(    npup*sizeof(float));
00282     for (min=MAXFLOAT,max=-MAXFLOAT,i=0; i<npup ;i++) {
00283       x0[i] = xd[i];
00284       y0[i] = yd[i];
00285       x1[i] = xd[i] + (glmlaxpup[i]-xd[i])*50;
00286       y1[i] = yd[i] + (glmlaypup[i]-yd[i])*50;
00287       min = MIN(min,rd[i]);
00288       max = MAX(max,rd[i]);
00289     }
00290     /* Levels */
00291     for (i=0; i<nlev ;i++) levels[i] = min + (i+1)*(max-min)/(nlev+1);
00292 
00293     sprintf(title,"Distortion residuals (rd -> robs)\n"
00294             "%d mpups, scale=50, %d levels [%.2f-%.2f]\n"
00295             "d=%.2f, angle=%.2f, alpha=%.2fe-8\n",
00296             npup,nlev,levels[0],levels[nlev-1],
00297             dist,angle/M_PI*180,alpha*1e8);
00298 #ifdef USE_R4_DIST
00299     sprintf(title,"%sbeta=%.2fe-16",title,beta*1e16);
00300 #endif
00301     plot_title(title);
00302 
00303     /* Triangulation: use x0 and y0 arrays (npup+3 points) */
00304     nmax = 2*npup + 1;
00305     i1ray = (int *)malloc(nmax*sizeof(int));
00306     i2ray = (int *)malloc(nmax*sizeof(int));
00307     i3ray = (int *)malloc(nmax*sizeof(int));
00308     ntri = triang(x0,y0,npup, i1ray,i2ray,i3ray, nmax); 
00309 
00310     /* Contours (contri cannot label contours) */
00311     for (i=0; i<nlev ;i++) {
00312       setclr((i+1) * 25);
00313       contri(xd,yd,rd,npup, i1ray,i2ray,i3ray,ntri, levels[i]);
00314     }
00315     //    confll(xd,yd,rd,npup, i1ray,i2ray,i3ray,ntri, levels,nlev); /* Shaded contours */
00316     free(i1ray); free(i2ray); free(i3ray);
00317     /* Mpups */
00318     color("fore");
00319     for (i=0; i<npup ;i++) rlsymb(3,glmlaxpup[i],glmlaypup[i]);
00320     field(xd,yd, x1,y1, npup, 1301);         /* Vector field */
00321     free(x0); free(y0); free(x1); free(y1);
00322     /* Center of distortion */
00323     color("red"); rlsymb(1,xc,yc);
00324     /* Graphic output conclusion */
00325     disfin();
00326   }
00327 #endif /* HAVE_LIBDISLIN */
00328 
00329 #undef dist
00330 #undef angle
00331 #undef xc
00332 #undef yc
00333 #undef alpha
00334 #ifdef USE_R4_DIST
00335 #undef beta
00336 #endif
00337 }
00338 
00339 /* ##### Global PSF ############################## */
00340 
00341 /* === Doxygen Comment ======================================= */
00351 /* =========================================================== */
00352 
00353 double set_geoPSF(float scale)
00354 {
00355   int i;
00356   double norm;
00357   
00358   switch (glnggeo) {
00359   case 3: 
00360     glIgeo[0] = GEOPUP3G_I0;       /* Central gaussian */
00361     glXgeo[0] = 0.; 
00362     glSgeo[0] = GEOPUP3G_S0*scale;
00363     glIgeo[1] = GEOPUP3G_I1;       /* Left gaussian    */
00364     glXgeo[1] =-GEOPUP3G_X1*scale; 
00365     glSgeo[1] = GEOPUP3G_S1*scale;
00366     glIgeo[2] = GEOPUP3G_I1;       /* Right gaussian   */
00367     glXgeo[2] = GEOPUP3G_X1*scale; 
00368     glSgeo[2] = GEOPUP3G_S1*scale;
00369     break;
00370   case 5:
00371     glIgeo[0] = GEOPUP5G_I0;       /* Central gaussian   */
00372     glXgeo[0] = 0.; 
00373     glSgeo[0] = GEOPUP5G_S0*scale;
00374     glIgeo[1] = GEOPUP5G_I1;       /* 1st-left gaussian  */
00375     glXgeo[1] =-GEOPUP5G_X1*scale; 
00376     glSgeo[1] = GEOPUP5G_S1*scale;
00377     glIgeo[2] = GEOPUP5G_I1;       /* 1st-right gaussian */
00378     glXgeo[2] = GEOPUP5G_X1*scale; 
00379     glSgeo[2] = GEOPUP5G_S1*scale;
00380     glIgeo[3] = GEOPUP5G_I2;       /* 2nd-left gaussian  */
00381     glXgeo[3] =-GEOPUP5G_X2*scale; 
00382     glSgeo[3] = GEOPUP5G_S2*scale;
00383     glIgeo[4] = GEOPUP5G_I2;       /* 2nd-right gaussian */
00384     glXgeo[4] = GEOPUP5G_X2*scale; 
00385     glSgeo[4] = GEOPUP5G_S2*scale;
00386     break;
00387   case 7:
00388     glIgeo[0] = GEOPUP7G_I0;       /* Central gaussian   */
00389     glXgeo[0] = 0.; 
00390     glSgeo[0] = GEOPUP7G_S0*scale;
00391     glIgeo[1] = GEOPUP7G_I1;       /* 1st-left gaussian  */
00392     glXgeo[1] =-GEOPUP7G_X1*scale; 
00393     glSgeo[1] = GEOPUP7G_S1*scale;
00394     glIgeo[2] = GEOPUP7G_I1;       /* 1st-right gaussian */
00395     glXgeo[2] = GEOPUP7G_X1*scale; 
00396     glSgeo[2] = GEOPUP7G_S1*scale;
00397     glIgeo[3] = GEOPUP7G_I2;       /* 2nd-left gaussian  */
00398     glXgeo[3] =-GEOPUP7G_X2*scale; 
00399     glSgeo[3] = GEOPUP7G_S2*scale;
00400     glIgeo[4] = GEOPUP7G_I2;       /* 2nd-right gaussian */
00401     glXgeo[4] = GEOPUP7G_X2*scale; 
00402     glSgeo[4] = GEOPUP7G_S2*scale;
00403     glIgeo[5] = GEOPUP7G_I3;       /* 3rd-left gaussian  */
00404     glXgeo[5] =-GEOPUP7G_X3*scale; 
00405     glSgeo[5] = GEOPUP7G_S3*scale;
00406     glIgeo[6] = GEOPUP7G_I3;       /* 3rd-right gaussian */
00407     glXgeo[6] = GEOPUP7G_X3*scale; 
00408     glSgeo[6] = GEOPUP7G_S3*scale;
00409     break;
00410   default:
00411     print_error("Geometric PSF: 3|5|7 gaussians (%d)",glnggeo);
00412     exit_session(ERR_NOIMPL);
00413   }
00414 
00415   for (norm=0., i=0; i<glnggeo; i++) norm += glIgeo[i]*M_SQRT2PI*glSgeo[i];
00416 
00417   return(norm);
00418 }
00419 
00420 /* === Doxygen Comment ======================================= */
00424 /* =========================================================== */
00425 
00426 double geoPSF(double x)
00427 {
00428   int i;
00429   double y=0;
00430   
00431   for (i=0; i<glnggeo; i++) 
00432     y += glIgeo[i]*exp(-0.5*SQ( (x-glXgeo[i])/glSgeo[i] ));
00433 
00434   return(y);
00435 }
00436 
00437 /* === Doxygen Comment ======================================= */
00451 /* =========================================================== */
00452 
00453 long nllsqfit_globalPSF(long *mode, long *npts, long *npar, long *ldfj, 
00454                         double *par, double fvec[], double fjac[], 
00455                         long *nstate, long *iuser, double x_n_sig[])
00456 {
00457   long i,j,k,g,p, npup, npsf, ngeo;
00458   double qe, qs, Ge, Gs, xe, xs, profil;
00459   double *sig, **Xs, **Is, **Ss;
00460 
00461   /* Global PSF components (i=0..npsf-1): intensity (for i>0) and sigma */
00462 #define Ipsf(i) (i==0 ? 1. : par[2*i-1])
00463 #define Spsf(i)              par[2*i]
00464   /* Central mpups (i=0..npup-1): intensity and position */
00465 #define Ipup(i) (par[2*npsf+2*i-1]*FMP_IPSFNORM)
00466 #define Xpup(i) (par[2*npsf+2*i  ]*FMP_XPSFNORM + glpcen[i]/2.)
00467 
00468   ngeo = glnggeo;
00469   npsf = glngpsf;
00470   npup = ((*npar) - 2*npsf + 1)/2; /* Number of central mpups */
00471 
00472   sig  = x_n_sig + (*npts); /* sig in the 2nd half of x_n_sig */
00473 
00474   if (npsf == 0) { /* That should never be the case here */
00475     print_error("Global PSF has 0 gaussian component!");
00476     *mode = -1;
00477     return(UNKNOWN);
00478   }
00479 
00480   alloc2d(&Is,ngeo,npsf,DOUBLE);
00481   alloc2d(&Ss,ngeo,npsf,DOUBLE);
00482   alloc2d(&Xs,npup,ngeo,DOUBLE);
00483 
00484   /* Convolution of the geometric PSF (ngeo gaussians) with the
00485      up-to-date global PSF model (npsf gaussians) */
00486   for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00487     for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00488       Ss[g][p] = hypot(Spsf(p),glSgeo[g]);
00489       Is[g][p] = M_SQRT2PI*Ipsf(p)*glIgeo[g]*Spsf(p)*glSgeo[g]/Ss[g][p];
00490     }
00491   }
00492 
00493   /* Loop on central mpups */
00494   for (k=i=0; i<npup; i++) { /* k=0...*npts-1 */
00495     xs = -0.5; /* pixel boundaries */
00496     xe =  0.5;
00497 
00498     /* Loop on geom PSF components */
00499     for (g=0; g<ngeo; g++) Xs[i][g] = Xpup(i) - glXgeo[g];
00500 
00501     for (j=0; j<glpcen[i]; j++, k++) { /* Loop on mpup profile pixels */
00502 
00503       /* Initialization */
00504       if (*mode != 1) fvec[k] = 0.;
00505       if (*mode != 0) for (p=0; p<*npar; p++) fjac[k+p*(*ldfj)] = 0.;
00506 
00507       for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00508         for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00509           qs = (xs - Xs[i][g])/Ss[g][p]/M_SQRT2;
00510           qe = (xe - Xs[i][g])/Ss[g][p]/M_SQRT2;
00511           Gs = Ipup(i)*Is[g][p]*exp(-SQ(qs));
00512           Ge = Ipup(i)*Is[g][p]*exp(-SQ(qe));
00513           profil = sqrt(M_PI_2)*Ipup(i)*Is[g][p]*Ss[g][p] *
00514             (gsl_sf_erf(qe) - gsl_sf_erf(qs));
00515           if (*mode != 1) fvec[k] += profil;
00516           if (*mode == 0) continue;
00517           /* df/dIpsf */
00518           if (p>0) fjac[k + (2*p-1)*(*ldfj)] += profil/Ipsf(p);
00519           /* df/dSpsf */
00520           fjac[k + (2*p)*(*ldfj)] += 
00521             ( (SQ(glSgeo[g]/Spsf(p))+1.)*profil -
00522               ((xe-Xs[i][g])*Ge - (xs-Xs[i][g])*Gs) ) * Spsf(p)/SQ(Ss[g][p]);
00523           /* df/dIpup */
00524           fjac[k + (2*npsf+2*i-1)*(*ldfj)] += profil/Ipup(i)*FMP_IPSFNORM;
00525           /* df/dXpup */
00526           fjac[k + (2*npsf+2*i  )*(*ldfj)] -= (Ge - Gs)     *FMP_XPSFNORM;
00527         }
00528       }
00529       /* Error normalization */
00530       fvec[k] /= sig[k];
00531       for (p=0; p<*npar; p++) fjac[k+p*(*ldfj)] /= sig[k];
00532 
00533       xs += 1.; /* next pixel */
00534       xe += 1.;
00535     } /* Loop on mpup profile pixels */
00536   }   /* Loop on central mpups */
00537 
00538   free2d(&Is,DOUBLE); free2d(&Ss,DOUBLE); free2d(&Xs,DOUBLE);
00539 
00540   return(OK);
00541 
00542 #undef Ipsf
00543 #undef Spsf
00544 #undef Ipup
00545 #undef Xpup
00546 }
00547 
00548 /* === Doxygen Comment ======================================= */
00552 /* =========================================================== */
00553 
00554 double globalPSF(double x, double psfpar[])
00555 {
00556   int i;
00557   double y=0;
00558   
00559   for (i=0; i<glngpsf; i++) 
00560     y += psfpar[2*i+1]*exp(-0.5*SQ( x/psfpar[2*i+2] ));
00561 
00562   return(y);
00563 }
00564 
00565 /* === Doxygen Comment ======================================= */
00574 /* =========================================================== */
00575 
00576 double deconv_globalPSF(int ngpsf, double psfpar[], double decfac)
00577 {
00578   int i,i0;
00579   double sigmin,Sdec,Idec,Sint, norm;
00580   
00581   print_msg("   Partial deconvolution factor: %.2f",decfac);
00582 
00583   /* Find the smallest sigma, used as ref. for the partial deconvolution */
00584   for (sigmin=MAXDOUBLE, i0=i=0; i<ngpsf; i++) { /* Loop on global PSF comp. */
00585     if (psfpar[2*i+2] < sigmin) {
00586       i0 = i;
00587       sigmin = psfpar[2*i+2];
00588     }
00589   }
00590 
00591   Sdec = sigmin*sqrt(decfac);
00592   Sint = sqrt(SQ(sigmin) - SQ(Sdec));
00593   Idec = sigmin*psfpar[2*i0+1]/(M_SQRT2PI*Sdec*Sint);
00594   print_msg("   Partial deconvolution gaussian (from C%d): "
00595             "Idec = %.3f, Sdec = %.3f",i0, Idec, Sdec);
00596 
00597   for (norm=0., i=0; i<ngpsf; i++) {  /* Loop on global PSF components */
00598     if (DEBUG) 
00599       print_msg("   C%d Before: Ipsf = %.3f, Spsf = %.3f",
00600                 i,psfpar[2*i+1],psfpar[2*i+2]);
00601     Sint          = sqrt(SQ(psfpar[2*i+2]) - SQ(Sdec));
00602     psfpar[2*i+1] = psfpar[2*i+1]*psfpar[2*i+2]/(M_SQRT2PI*Sdec*Sint*Idec);
00603     psfpar[2*i+2] = Sint;
00604     norm += psfpar[2*i+1]*M_SQRT2PI*psfpar[2*i+2];
00605     if (DEBUG) 
00606       print_msg("      After:  Ipsf = %.3f, Spsf = %.3f",
00607                 psfpar[2*i+1],psfpar[2*i+2]);
00608   }
00609 
00610   return(norm);
00611 }
00612 
00613 /* ##### Local PSF ############################## */
00614 
00615 /* === Doxygen Comment ======================================= */
00624 /* =========================================================== */
00625 
00626 double set_corePSF(double psfpar[]) 
00627 {
00628   int i,j;
00629   double norm;
00630   
00631   /* Use the global PSF derived previously: core PSF = geom PSF (x) glob PSF */
00632   if (glngpsf > 0) {
00633     for (norm=0,i=0; i<glnggeo; i++) { /* Loop on geom PSF components   */
00634       for (j=0; j<glngpsf; j++) {      /* Loop on global PSF components */
00635         glSpsf[i][j] = hypot(psfpar[2*j+2],glSgeo[i]);
00636         glIpsf[i][j] = M_SQRT2PI*psfpar[2*j+1]*glIgeo[i]*
00637           psfpar[2*j+2]*glSgeo[i] / glSpsf[i][j];
00638         norm += glIpsf[i][j]*M_SQRT2PI*glSpsf[i][j];
00639       }
00640     }
00641   }
00642   /* Do *not* use the global PSF derived previously: core PSF = geom PSF */
00643   else {
00644     print_warning("Global PSF discarded: core PSF = geom PSF");
00645     for (norm=0,i=0; i<glnggeo; i++) { /* Loop on geom PSF components */
00646       glSpsf[i][0] = glSgeo[i];
00647       glIpsf[i][0] = glIgeo[i];
00648       norm += glIpsf[i][0]*M_SQRT2PI*glSpsf[i][0];
00649     }
00650   }
00651 
00652   return(norm);
00653 }
00654 
00655 /* === Doxygen Comment ======================================= */
00665 /* =========================================================== */
00666 
00667 long nllsqfit_localPSFn(long *mode, long *npts, long *npar, long *ldfj, 
00668                         double *par, double fvec[], double fjac[], 
00669                         long *nstate, long *iuser, double x_n_sig[])
00670 {
00671   long j,g,p,l, npsf, ngeo, nloc;
00672   double qe, qs, Ge, Gs, xe, xs, profil;
00673   double *x, *sig, *Xs, *Is, *Ss;
00674 
00675 #define Xpup    (par[0]    *FMP_XPUPNORM + (*npts)/2.)
00676 #define Ipup(i) (par[2*i+1]*FMP_IPUPNORM)
00677 #define Spup(i) (par[2*i+2]*FMP_SPUPNORM)
00678 
00679   ngeo = glnggeo;
00680   npsf = MAX(1,glngpsf); /* Handle the case: glngpsf = 0 */
00681   nloc = glngloc;
00682 
00683   x    = x_n_sig;           /* x   in the 1st half of x_n_sig */
00684   sig  = x_n_sig + (*npts); /* sig in the 2nd half of x_n_sig */
00685 
00686   Xs = (double *)malloc(ngeo*sizeof(double));
00687   Is = (double *)malloc(ngeo*npsf*nloc*sizeof(double));
00688   Ss = (double *)malloc(ngeo*npsf*nloc*sizeof(double));
00689 
00690 #define Is(g,p,l) (Is[g*npsf*nloc + p*nloc + l])
00691 #define Ss(g,p,l) (Ss[g*npsf*nloc + p*nloc + l])
00692 
00693   /* Convolution of the core PSF (ngeo*npsf gaussians) with the
00694      up-to-date local PSF model (nloc gaussian) */
00695   for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00696     Xs[g] = Xpup - glXgeo[g];
00697     for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00698       for (l=0; l<nloc; l++) {   /* Loop on local PSF components  */
00699         Ss(g,p,l) = hypot(glSpsf[g][p],Spup(l));
00700         Is(g,p,l) = 
00701           M_SQRT2PI*Ipup(l)*glIpsf[g][p]*Spup(l)*glSpsf[g][p]/Ss(g,p,l);
00702       }
00703     }
00704   }
00705 
00706   xs = -0.5; /* pixel boundaries */
00707   xe =  0.5;
00708 
00709   /* Loop on mpup profile pixels */
00710   for (j=0; j<*npts; j++) {
00711     /* Initialization */
00712     if (*mode != 1) fvec[j] = 0.;
00713     if (*mode != 0) for (p=0; p<*npar; p++) fjac[j+p*(*ldfj)] = 0.;
00714 
00715     for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00716       for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00717         for (l=0; l<nloc; l++) {   /* Loop on local PSF components  */
00718           qs = (xs - Xs[g])/Ss(g,p,l)/M_SQRT2;
00719           qe = (xe - Xs[g])/Ss(g,p,l)/M_SQRT2;
00720           Gs = Is(g,p,l)*exp(-SQ(qs));
00721           Ge = Is(g,p,l)*exp(-SQ(qe));
00722           profil = sqrt(M_PI_2)*Is(g,p,l)*Ss(g,p,l)*
00723             (gsl_sf_erf(qe) - gsl_sf_erf(qs));
00724           if (*mode != 1) fvec[j] += profil;
00725           if (*mode == 0) continue;
00726           /* df/dXpup */
00727           fjac[j] -= (Ge - Gs)             *FMP_XPUPNORM;
00728           /* df/dIpup(l) */
00729           fjac[j + (2*l+1)*(*ldfj)] += profil/Ipup(l)*FMP_IPUPNORM;
00730           /* df/dSpup(l) */
00731           fjac[j + (2*l+2)*(*ldfj)] +=                FMP_SPUPNORM*
00732             ( (SQ(glSpsf[g][p]/Spup(l)) + 1.)*profil -
00733               ((xe-Xs[g])*Ge - (xs-Xs[g])*Gs) ) * Spup(l)/SQ(Ss(g,p,l));
00734         }
00735       }
00736     }
00737     /* Normalization */
00738     fvec[j] /= sig[j];
00739     for (p=0; p<*npar; p++) fjac[j+p*(*ldfj)] /= sig[j];
00740 
00741     xs += 1.; /* next pixel */
00742     xe += 1.;
00743   }
00744 
00745   free(Xs); free(Is); free(Ss);
00746 
00747   return(OK);
00748 
00749 #undef Is
00750 #undef Ss
00751 
00752 #undef Ipup
00753 #undef Xpup
00754 #undef Spup
00755 }
00756 
00757 /* === Doxygen Comment ======================================= */
00773 /* =========================================================== */
00774 
00775 double final_mpup(double xarr[], int npts, double start, double step, 
00776                   double xloc[], double dxloc[], double xfit[])
00777 {
00778   int i,status,ngloc;
00779   double chi2, gof, rms, *lxloc;
00780 
00781   ngloc = glngloc;
00782 
00783   /* Boudaries on parameters */
00784   lxloc = (double *)calloc((1+2*ngloc),sizeof(double));
00785   lxloc[0] = -MAXDOUBLE; /* All local PSF param >0 except X */
00786 
00787   /* Initial guess normalization */
00788   xloc[0] = (xloc[0]-npts/2.)/FMP_XPUPNORM;
00789   xloc[1] /= FMP_IPUPNORM;
00790   xloc[2] /= FMP_SPUPNORM;
00791   for (i=1; i<ngloc; i++) {
00792     xloc[2*i+1] /= FMP_IPUPNORM;
00793     xloc[2*i+2] /= FMP_SPUPNORM;
00794   }
00795 
00796   /* Least-squares */
00797   /* dxloc[0] = 0; *//* Eastern egg: Print-out: 1, Derivatives: 2 */
00798 
00799   status = nllsqfit_bnd(NULL,xarr,NULL,npts,xloc,lxloc,NULL,1+2*ngloc,
00800                         *nllsqfit_localPSFn, xfit, NULL, dxloc, &chi2,&gof);
00801 
00802   /* Renormalization */
00803   if (!status) {
00804     xloc[0] = xloc[0]*FMP_XPUPNORM+npts/2. + start; dxloc[0] *= FMP_XPUPNORM;
00805     xloc[1] *= FMP_IPUPNORM;                        dxloc[1] *= FMP_IPUPNORM;
00806     xloc[2] *= FMP_SPUPNORM;                        dxloc[2] *= FMP_SPUPNORM;
00807     for (i=1; i<ngloc; i++) {
00808       xloc[2*i+1] *= FMP_IPUPNORM; dxloc[2*i+1] *= FMP_IPUPNORM;
00809       xloc[2*i+2] *= FMP_SPUPNORM; dxloc[2*i+2] *= FMP_SPUPNORM;
00810     }
00811     rms = sqrt(chi2/npts);
00812   }
00813   else {
00814     for (i=0; i<3; i++)     xloc[i]     = 0.;
00815     for (i=1; i<ngloc; i++) xloc[2*i+1] = 0.;
00816     rms = -status;
00817     for (i=0; i<npts; i++) xfit[i] = 0.;      /* Zero-ing the fit profile */
00818   }
00819 
00820   free(lxloc);
00821 
00822   return(rms);
00823 }
00824 

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