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

extract/source/fit_mpupgeo.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 /* === Doxygen Comment ======================================== */
00025 /* =========================================================== */
00026 
00027 long nllsqfit_mpupgeo(long *mode, long *npts, long *npar, long *ldfj, 
00028                       double par[], double f[], double fjac[], 
00029                       long *nstate, long *iuser, double x_n_sig[])
00030 {
00031   int i,n, ng1s;
00032   double *x, *sig, *xx, *ex;
00033 
00034   x  = x_n_sig;           /* x   in the 1st half of x_n_sig */
00035   sig= x_n_sig + (*npts); /* sig in the 2nd half of x_n_sig */
00036 
00037   /* Number of gaussians on *one* side (*not* including central one) */
00038   ng1s = (*npar - 2)/3;
00039 
00040   xx = (double*)malloc((2*ng1s+1)*sizeof(double));
00041   ex = (double*)malloc((2*ng1s+1)*sizeof(double));
00042 
00043   for (n=0; n<*npts; n++) {
00044     /* Function & Derivatives */
00045     xx[0] = x[n]/par[1];
00046     ex[0] = exp(-SQ(xx[0])/2.);
00047     f[n]  = par[0]*ex[0];
00048     fjac[n]           = ex[0];                         /* df/dI0 */
00049     fjac[n+  (*npts)] = par[0]*ex[0]*SQ(xx[0])/par[1]; /* df/dS0 */
00050     for (i=1; i<=ng1s; i++) {
00051       xx[2*i-1] = (x[n]-par[3*i])/par[3*i+1]; /* - */
00052       ex[2*i-1] = exp(-SQ(xx[2*i-1])/2.);
00053       xx[2*i  ] = (x[n]+par[3*i])/par[3*i+1]; /* + */
00054       ex[2*i  ] = exp(-SQ(xx[2*i  ])/2.);
00055       f[n] += par[3*i-1]*(ex[2*i-1] + ex[2*i]);
00056 
00057       fjac[n+(3*i-1)*(*npts)] = ex[2*i-1] + ex[2*i];     /* df/dIi */
00058       fjac[n+(3*i  )*(*npts)] = par[3*i-1]/par[3*i+1]*
00059         (ex[2*i-1]*xx[2*i-1] - ex[2*i]*xx[2*i]);         /* df/dXi */
00060       fjac[n+(3*i+1)*(*npts)] = par[3*i-1]/par[3*i+1]*
00061         (ex[2*i-1]*SQ(xx[2*i-1]) + ex[2*i]*SQ(xx[2*i])); /* df/dSi */
00062     }
00063   }
00064 
00065   free(xx); free(ex);
00066 
00067   return(OK);
00068 }
00069 
00070 /* === Doxygen Comment ======================================= */
00074 /* =========================================================== */
00075 
00076 void print_mpupgeo_par(double par[], double dpar[], int ng1s)
00077 {
00078   int i,n,ng;
00079 
00080   if (ng1s < 0) {
00081     n = -ng1s;
00082     ng = 2*n+1;
00083     print_msg("#define GEOPUP%dG_I0 %.3f",ng,1.);
00084     print_msg("#define GEOPUP%dG_S0 %.3f",ng,par[1]);
00085     for (i=1; i<=n; i++) {
00086       print_msg("#define GEOPUP%dG_I%d %.3f",ng,i,par[3*i-1]/par[0]);
00087       print_msg("#define GEOPUP%dG_X%d %.3f",ng,i,par[3*i  ]);
00088       print_msg("#define GEOPUP%dG_S%d %.3f",ng,i,par[3*i+1]);
00089     }    
00090   }
00091   else if (dpar == NULL) { /* No errors: initial guess */
00092     print_msg("Gaussian #0: I=%.2f, S=%.2f",par[0],par[1]);
00093     for (i=1; i<=ng1s; i++) {
00094       print_msg("Gaussian #%d: I=%.2f, X=%.2f, S=%.2f",
00095                 i,par[3*i-1],par[3*i],par[3*i+1]);
00096     }
00097   }
00098   else {
00099     print_msg("Gaussian #0: I=%.2f +-%.2f, S=%.2f +-%.2f",
00100               par[0],dpar[0],par[1],dpar[1]);
00101     for (i=1; i<=ng1s; i++) {
00102       print_msg("Gaussian #%d: I=%.2f +-%.2f, X=%.2f +-%.2f, S=%.2f +-%.2f",
00103                 i,par[3*i-1],dpar[3*i-1],par[3*i],dpar[3*i],
00104                 par[3*i+1],dpar[3*i+1]);
00105     }
00106   }
00107 }
00108 
00109 /* === Doxygen Comment ======================================= */
00114 /* =========================================================== */
00115 
00116 void save_mpupgeo_par(double par[],int ng1s, char *name,int npts,
00117                       float start,float step)
00118 {
00119   SPECTRUM tmp;
00120   char tmpname[lg_name+1];
00121   int i,n;
00122   float x,y;
00123 
00124   print_msg("Writing global fit to %s",name);
00125   sprintf(tmpname,"%s",name);
00126   create_spec(&tmp,tmpname,npts,start,step,FLOAT,"Mpup profile FIT",NULL);
00127   for (n=0; n<npts; n++) {
00128     x = coord_spec(&tmp,n);
00129     y = par[0]*exp(-SQ(x/par[1])/2.);
00130     for (i=1; i<=ng1s; i++) {
00131       y += par[3*i-1]*(exp(-SQ((x-par[3*i])/par[3*i+1])/2.) + 
00132                        exp(-SQ((x+par[3*i])/par[3*i+1])/2.));
00133     }
00134     WR_spec(&tmp, n, y);
00135   }
00136   close_spec(&tmp);
00137 
00138   print_msg("Writing individual components [0-%d] to %s_i",ng1s,name);
00139   sprintf(tmpname,"%s_0",name);
00140   create_spec(&tmp,tmpname,npts,start,step,FLOAT,"Mpup profile FIT",NULL);
00141   for (n=0; n<npts; n++) {
00142     x = coord_spec(&tmp,n);
00143     y = par[0]*exp(-SQ(x/par[1])/2.);
00144     WR_spec(&tmp, n, y);
00145   }
00146   close_spec(&tmp);
00147 
00148   for (i=1; i<=ng1s; i++) {
00149     sprintf(tmpname,"%s_%d",name,i);
00150     create_spec(&tmp,tmpname,npts,start,step,FLOAT,"Mpup profile FIT",NULL);
00151     for (n=0; n<npts; n++) {
00152       x = coord_spec(&tmp,n);
00153       y = par[3*i-1]*(exp(-SQ((x-par[3*i])/par[3*i+1])/2.) + 
00154                       exp(-SQ((x+par[3*i])/par[3*i+1])/2.));
00155       WR_spec(&tmp, n, y);
00156     }
00157     close_spec(&tmp);
00158   }
00159 }
00160 
00161 /* === Doxygen Comment ======================================= */
00176 /* =========================================================== */
00177 
00178 int main(int argc, char **argv)
00179 {
00180   char **argval, **arglabel;
00181 
00182   SPECTRUM raw,fit;
00183   char rawname[lg_name+1], fitname[lg_name+1];
00184   int i,ng,npts,ng1s,npar,status,niter,nitermax;
00185   float d1,d2,fac,N,start,step;  
00186   double *x,*y,*f,*par,*dpar,*lpar;
00187   double y1,y2,chi2,gof;
00188 
00189   set_purpose("Multi-gaussian 1D-fit of the geometric pupil profile");
00190   set_arglist("-d1 2239 -d2 613 -ng 5 -n 200 -fac 1.6 "
00191               "-raw mpup_profile -fit mpup_profile_5gauss");
00192 
00193   init_snifs("$Name:  $");
00194   init_session(argv,argc,&arglabel,&argval);
00195 
00196   if (DEBUG) {
00197     print_msg("$Id: fit_mpupgeo.c,v 1.8 2004/10/06 16:56:32 ycopin Exp $");
00198     print_msg("Compilation: %s, %s",__DATE__,__TIME__);
00199   }
00200 
00201   /* Parameters ============================== */
00202 
00203   get_argval(0,"%f",&d1);   /* primary mirror   */
00204   get_argval(1,"%f",&d2);   /* secondary mirror */
00205   get_argval(2,"%d",&ng);   /* nb of gaussians  */
00206   get_argval(3,"%d",&npts); /* nb of points     */
00207   get_argval(4,"%f",&fac);  /* expansion factor */
00208   get_argval(5,"%s",rawname); /* raw profile    */
00209   get_argval(6,"%s",fitname); /* adjusted prof. */
00210 
00211   if (!(ng%2)) {
00212     print_error("Nb of gaussians (%d) should be odd",ng);
00213     exit_session(ERR_BAD_PARAM);
00214   }
00215 
00216   print_msg("o Primary mirror (mm): %.1f",d1);
00217   print_msg("o Central obsc. (mm): %.1f",d2);
00218   print_msg("o Nb of gaussians: %d",ng);
00219   print_msg("o Nb of points: %d",npts);
00220   print_msg("o Expansion factor: %.1f",fac);
00221   print_msg("o Raw mpup profile: %s",rawname);
00222   print_msg("o Adjusted mpup profile: %s",fitname);
00223 
00224   /* Raw profile ============================== */
00225 
00226   N = d2/d1;
00227   start = -0.5*fac;
00228   step  = -2.*start/(npts-1);
00229 
00230   x = (double *)malloc(npts*sizeof(double));
00231   y = (double *)malloc(npts*sizeof(double));
00232   f = (double *)malloc(npts*sizeof(double));
00233 
00234   create_spec(&raw,rawname,npts,start,step,FLOAT,"Mpup profile",NULL);
00235   for (i=0; i<npts; i++) {
00236     x[i] = start + i*step;
00237     y1   = (ABS(2*x[i])  <1) ? sqrt(1-SQ(2*x[i])  ) : 0.; /* Primary       */
00238     y2   = (ABS(2*x[i]/N)<1) ? sqrt(1-SQ(2*x[i]/N)) : 0.; /* Central obsc. */
00239     y[i] = y1 - N*y2;
00240     WR_spec(&raw, i, y[i]);
00241   }
00242   close_spec(&raw);
00243 
00244   /* Adjustment ============================== */
00245 
00246   /* Parameters: I0,S0,I1,X1,S1,I2,X2,S2,... = 2 + (ng-1)/2*3 */
00247 
00248   ng1s = (ng-1)/2;    /* nb of gaussians on *one* side (central one not incl.)*/
00249   npar = 2 + ng1s*3;
00250   par  = (double *)malloc(npar*sizeof(double));
00251   lpar = (double *)calloc(npar,sizeof(double));
00252   dpar = (double *)malloc(npar*sizeof(double));
00253 
00254   par[0] = 2./ng;              /* Central gaussian intensity */
00255   par[1] = 0.5/(ng1s+2);        /*                  sigma     */
00256   for (i=1; i<=ng1s; i++) {
00257     par[3*i-1] = 2.5/ng;        /* i-th gaussian intensity */
00258     par[3*i  ] = (i+1)*par[1];  /*               position  */
00259     par[3*i+1] = par[1];        /*               sigma     */
00260   }
00261 
00262   print_msg("Initial estimates:");
00263   print_mpupgeo_par(par,NULL,ng1s);
00264 
00265   if (DEBUG) save_mpupgeo_par(par,ng1s,"dbg_fit",npts,start,step);
00266 
00267   /* Simple least-sq. fit */
00268   status = nllsqfit_bnd(x, y, NULL, npts, par, lpar, NULL, npar,
00269                         *nllsqfit_mpupgeo, f, NULL, dpar, &chi2,&gof);
00270 
00271   if (status == 4) { /* No enough iterations to converge, try again */
00272     nitermax = 3;
00273     niter = 1;
00274     print_msg("Fit #%d/%d (status=%d, chi2=%f)",niter,nitermax,status,chi2);
00275     while (status == 4 && niter < nitermax) {
00276       status = nllsqfit_bnd(x, y, NULL, npts, par, lpar, NULL, npar,
00277                             *nllsqfit_mpupgeo, f, NULL, dpar, &chi2,&gof);
00278       niter ++;
00279       print_msg("Fit #%d/%d (status=%d, chi2=%f)",niter,nitermax,status,chi2);
00280     }
00281   }
00282 
00283   gof = sqrt(chi2/npts);
00284   print_msg("Fit (status=%d, chi2=%f, RMS=%f):",status,chi2,gof);
00285   print_mpupgeo_par(par,dpar, ng1s);
00286   print_msg("Normalized parameters:");
00287   print_mpupgeo_par(par,dpar,-ng1s);
00288 
00289   save_mpupgeo_par(par,ng1s,fitname,npts,start,step);
00290   open_spec(&fit,fitname,"IO");
00291   WR_desc(&fit,"NGAUSS",INT,1,   &ng);
00292   WR_desc(&fit,"CHI2",  DOUBLE,1,&chi2);
00293   WR_desc(&fit,"RMS",   DOUBLE,1,&gof);
00294   close_spec(&fit);
00295 
00296   free(par); free(dpar);
00297   free(x); free(y); free(f);
00298 
00299   exit_session(OK);
00300   return(OK);
00301 }

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