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

extract/lib/deprecated.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 ======================================= */
00033 /* =========================================================== */
00034 
00035 long nllsqfit_localPSF(long *mode, long *npts, long *npar, long *ldfj, 
00036                        double par[], double fvec[], double fjac[], 
00037                        long *nstate, long *iuser, double x_n_sig[])
00038 {
00039   long j,g,p, npsf, ngeo;
00040   double qe, qs, Ge, Gs, xe, xs, profil;
00041   double *x, *sig, *Xs, **Is, **Ss;
00042 
00043 #define Xpup (par[0]*FMP_XPUPNORM + (*npts)/2.)
00044 #define Ipup (par[1]*FMP_IPUPNORM)
00045 #define Spup (par[2]*FMP_SPUPNORM)
00046 
00047   ngeo = glnggeo;
00048   npsf = MAX(1,glngpsf); /* Handle the case: glngpsf = 0 */
00049 
00050   x    = x_n_sig;           /* x   in the 1st half of x_n_sig */
00051   sig  = x_n_sig + (*npts); /* sig in the 2nd half of x_n_sig */
00052 
00053   Xs = (double *)malloc(ngeo*sizeof(double));
00054   alloc2d(&Is,ngeo,npsf,DOUBLE);
00055   alloc2d(&Ss,ngeo,npsf,DOUBLE);
00056 
00057   /* Convolution of the core PSF (ngeo*npsf gaussians) with the
00058      up-to-date local PSF model (1 gaussian) */
00059   for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00060     Xs[g] = Xpup - glXgeo[g];
00061     for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00062       Ss[g][p] = hypot(glSpsf[g][p],Spup);
00063       Is[g][p] = M_SQRT2PI*Ipup*glIpsf[g][p]*Spup*glSpsf[g][p]/Ss[g][p];
00064     }
00065   }
00066 
00067   xs = -0.5; /* pixel boundaries */
00068   xe =  0.5;
00069 
00070   /* Loop on mpup profile pixels */
00071   for (j=0; j<*npts; j++) {
00072     /* Initialization */
00073     if (*mode != 1) fvec[j] = 0.;
00074     if (*mode != 0) for (p=0; p<*npar; p++) fjac[j+p*(*ldfj)] = 0.;
00075 
00076     for (g=0; g<ngeo; g++) {       /* Loop on geom PSF components   */
00077       for (p=0; p<npsf; p++) {     /* Loop on global PSF components */
00078         qs = (xs - Xs[g])/Ss[g][p]/M_SQRT2;
00079         qe = (xe - Xs[g])/Ss[g][p]/M_SQRT2;
00080         Gs = Is[g][p]*exp(-SQ(qs));
00081         Ge = Is[g][p]*exp(-SQ(qe));
00082         profil = sqrt(M_PI_2)*Is[g][p]*Ss[g][p]*
00083           (gsl_sf_erf(qe) - gsl_sf_erf(qs));
00084         if (*mode != 1) fvec[j] += profil;
00085         if (*mode == 0) continue;
00086         /* df/dXpup */
00087         fjac[j          ] -= Ge - Gs;
00088         /* df/dIpup */
00089         fjac[j + (*ldfj)] += profil/Ipup;
00090         /* df/dSpup */
00091         fjac[j+2*(*ldfj)] += 
00092           ( (SQ(glSpsf[g][p]/Spup) + 1.)*profil -
00093             ((xe-Xs[g])*Ge - (xs-Xs[g])*Gs) ) * Spup/SQ(Ss[g][p]);
00094       }
00095     }
00096     /* Normalization */
00097     fvec[j] /= sig[j];
00098     fjac[j]           *= FMP_XPUPNORM/sig[j];
00099     fjac[j + (*ldfj)] *= FMP_IPUPNORM/sig[j];
00100     fjac[j+2*(*ldfj)] *= FMP_SPUPNORM/sig[j];
00101 
00102     xs += 1.; /* next pixel */
00103     xe += 1.;
00104   }
00105 
00106   free(Xs); free2d(&Is,DOUBLE); free2d(&Ss,DOUBLE);
00107 
00108   return(OK);
00109 
00110 #undef Xpup
00111 #undef Ipup
00112 #undef Spup
00113 }
00114 
00115 /* === Doxygen Comment ======================================= */
00126 /* =========================================================== */
00127 
00128 Channel snifs_channel(const float lbdaref) 
00129 {
00130   Channel channel;
00131   
00132   if (lbdaref>SNIFSB_FILT_L && lbdaref<SNIFSB_FILT_H) 
00133     channel = BLUE_CHANNEL;
00134   else if (lbdaref>SNIFSR_FILT_L && lbdaref<SNIFSR_FILT_H) 
00135     channel = RED_CHANNEL;
00136   else {
00137     print_error("Ref. wavelength %.1f AA does not match any channel",lbdaref);
00138     return(UNKNOWN_CHANNEL);
00139   }
00140 
00141   print_msg("Ref. wavelength %.1f AA corresponds to %s channel",lbdaref,
00142             channel_names[channel]);
00143 
00144   return(channel);
00145 }

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