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

extract/lib/optics_signature.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_poly.h>
00020 
00021 /* === Doxygen Comment ======================================= */
00048 /* =========================================================== */
00049 
00050 double blaze_function(const Grism *grism, int order, double lambda)
00051 {
00052   double ngl,nep, rholbda,nglsinA, i,r, Theta, BF;
00053   
00054   ngl = sellmeier_index(lambda,grism->prism.sellcoeffs);   /* Prism index */
00055   nep = sellmeier_index(lambda,grism->grating.sellcoeffs); /* Grating index */
00056 
00057   rholbda = grism->g_per_mm*(lambda*1e-7);   /* mm-1 * mm = no dimension */
00058   nglsinA = ngl*sin(grism->A);
00059 
00060   i = asin( nglsinA/nep             ) - grism->blaze; /* Eq. 7 & 8 */
00061   r = asin( nglsinA - order*rholbda ) - grism->blaze; /* Eq. 7 & 9 */
00062 
00063   Theta = M_PI/rholbda*cos(grism->blaze) * (nep*sin(i) - sin(r)); /* Eq. 3 */
00064   
00065   BF = SQ(sin(Theta)/Theta);                 /* Blaze function */
00066   
00067   return(BF);
00068 }
00069 
00070 /* === Doxygen Comment ======================================= */
00078 /* =========================================================== */
00079 
00080 double lateral_color_ampl(const CameraCollimator *element, double lambda) 
00081 {
00082   int i;
00083   double a, dl, pow;
00084 
00085   dl = lambda - element->LateralColor.lbdaref;
00086   for (a=0., pow=dl, i=0; i<3; i++) {
00087     a += element->LateralColor.coef[i]*pow;
00088     pow *= dl;
00089   }
00090   return(a);
00091 }
00092 
00093 /* === Doxygen Comment ======================================= */
00108 /* =========================================================== */
00109 
00110 double collimator_forward(const CameraCollimator *collimator, double gamma,
00111                           double lambda, double r)
00112 {
00113   double tantheta, rr;
00114 
00115   rr = r/collimator->focal;
00116   tantheta  = lateral_color_ampl(collimator, lambda)/gamma;
00117   tantheta += 1. + SQ(rr)*collimator->distor;
00118   tantheta *= rr;
00119 
00120   return(tantheta);
00121 }
00122 
00123 /* === Doxygen Comment ======================================= */
00136 /* =========================================================== */
00137 
00138 double camera_forward(const CameraCollimator *camera, double lambda, double tantheta)
00139 {
00140   double r;
00141 
00142   r  = lateral_color_ampl(camera, lambda);
00143   r += 1. + SQ(tantheta)*camera->distor;
00144   r *= tantheta*camera->focal;
00145 
00146   return(r);
00147 }
00148 
00149 /* === Doxygen Comment ======================================= */
00158 /* =========================================================== */
00159 
00160 double invert_camcoll(double y, double e, double b)
00161 {
00162   int nsol;
00163   double x0,x1,x2, x;
00164 
00165   /* Trivial cases */
00166   if (y==0) return(0.);
00167   else if (b==0 && e==0) {
00168     print_error("ARGH! No solution in invert_camcoll");
00169     exit_session(UNKNOWN);
00170   }
00171   else if (b==0) return(cbrt(y/e)); /* cbrt returns the cube-root */
00172   else if (e==0) return(y/b);
00173 
00174   nsol = gsl_poly_solve_cubic(0.,b/e,-y/e, &x0,&x1,&x2);
00175   
00176   if (nsol==3) {
00177     /*
00178     print_warning("ARGH! 3 real solutions in invert_camcoll:\n"
00179                   "      %f=x*(%f*x^2+%f) => x0=%f, x1=%f, x2=%f",
00180                   y,e,b,x0,x1,x2);
00181     x = MAXDOUBLE;
00182     if (x0>0. && x0<x) x = x0;
00183     if (x1>0. && x1<x) x = x1;
00184     if (x2>0. && x2<x) x = x2;
00185     if (x == MAXDOUBLE) {
00186       print_error("ARGH! No positive solution in invert_camcoll");
00187       exit_session(UNKNOWN);
00188     }
00189     else
00190       print_warning("ARGH! Returning smallest positive root %f...",x);
00191     */
00192     /* Returning central one */
00193     x = x1;
00194   }
00195   else x = x0;
00196 
00197   return(x);
00198 }
00199 
00200 /* === Doxygen Comment ======================================= */
00217 /* =========================================================== */
00218 
00219 #define camera_backward(camera,lambda,r) \
00220         invert_camcoll((r)/(camera)->focal,(camera)->distor, \
00221                        1.+lateral_color_ampl((camera),lambda))
00222 
00223 /* === Doxygen Comment ======================================= */
00242 /* =========================================================== */
00243 
00244 #define collimator_backward(collimator,gamma,lambda,tantheta) \
00245         (collimator)->focal * \
00246         invert_camcoll((tantheta),(collimator)->distor, \
00247                1.+lateral_color_ampl((collimator),(lambda))/(gamma))
00248 
00249 /* === Doxygen Comment ======================================= */
00257 /* =========================================================== */
00258 
00259 void rotation_2D(double *x, double *y, double theta) 
00260 {
00261   double costheta, sintheta, x1, y1;
00262 
00263   costheta = cos(theta);
00264   sintheta = sin(theta);
00265 
00266   x1 =  (*x)*costheta + (*y)*sintheta;
00267   y1 = -(*x)*sintheta + (*y)*costheta;
00268 
00269   *x = x1;
00270   *y = y1;
00271 }
00272 
00273 /* === Doxygen Comment ======================================= */
00277 /* =========================================================== */
00278 
00279 #define x_rotation(x,y,z,t) (rotation_2D((y),(z),(t)))
00280 
00281 /* === Doxygen Comment ======================================= */
00285 /* =========================================================== */
00286 
00287 #define y_rotation(x,y,z,t) (rotation_2D((x),(z),(t)))
00288 
00289 /* === Doxygen Comment ======================================= */
00293 /* =========================================================== */
00294 
00295 #define z_rotation(x,y,z,t) (rotation_2D((x),(y),(t)))
00296 
00297 /* === Doxygen Comment ======================================= */
00307 /* =========================================================== */
00308 
00309 void refraction(double *x, double *y, double *z, double n1, double n2) 
00310 {
00311   *x *= n1/n2;
00312   *y *= n1/n2;
00313   *z  = sqrt(1. - ( SQ(*x) + SQ(*y) ));
00314 }
00315 
00316 /* === Doxygen Comment ======================================= */
00331 /* =========================================================== */
00332 
00333 void grating_forward(double *x, double *y, double *z, 
00334                      double lambda, double n, double g_per_mm, int order) 
00335 {
00336   *x = (*x)*n;
00337   *y = (*y)*n - order * (lambda*1e-7) * g_per_mm;
00338   *z = sqrt(1. - (SQ(*x)+SQ(*y)));
00339 }
00340 
00341 /* === Doxygen Comment ======================================= */
00355 /* =========================================================== */
00356 
00357 void grating_backward(double *x, double *y, double *z, 
00358                       double lambda, double n, double g_per_mm, int order) 
00359 {
00360   double x1, y1, z1;
00361 
00362   x1 = (*x)/n;
00363   y1 = (*y + order * (lambda*1e-7) * g_per_mm)/n;
00364   z1 = sqrt(1. - (SQ(x1)+SQ(y1)));
00365 
00366   *x = x1;
00367   *y = y1;
00368   *z = z1;
00369 }
00370 
00371 /* === Doxygen Comment ======================================= */
00385 /* =========================================================== */
00386 
00387 double sellmeier_index(double lambda, const double sellcoeffs[])
00388 {
00389   int i;
00390   double n, l2;
00391 
00392   l2 = SQ(lambda*1.e-4);
00393   for (n=1., i=0; i<3; i++) n += sellcoeffs[i]*l2/(l2 - sellcoeffs[i+3]);
00394   n = sqrt(n);
00395 
00396   return(n);
00397 }
00398 
00399 /* === Doxygen Comment ======================================= */
00410 /* =========================================================== */
00411 
00412 void grism_forward(const Grism *grism, double lambda, 
00413                    double *x, double *y, double *z, int order)
00414 {
00415   double n;
00416   
00417   n = sellmeier_index(lambda, grism->prism.sellcoeffs);
00418 
00419   /* Get aligned w/ grism */
00420   x_rotation(x,y,z, grism->tiltx);
00421   y_rotation(x,y,z, grism->tilty);
00422   z_rotation(x,y,z, M_PI+grism->rotz);
00423 
00424   /* Air/glass interface */
00425   refraction(x,y,z, 1.,n);
00426 
00427   /* Prism */
00428   x_rotation(x,y,z, grism->A);
00429 
00430   /* Grating */
00431   grating_forward(x,y,z, lambda,n,grism->g_per_mm,order);
00432 
00433   /* Back to original orientation */
00434   x_rotation(x,y,z, -grism->A);
00435   z_rotation(x,y,z, -(M_PI+grism->rotz));
00436   y_rotation(x,y,z, -grism->tilty);
00437   x_rotation(x,y,z, -grism->tiltx);
00438 }
00439 
00440 /* === Doxygen Comment ======================================= */
00451 /* =========================================================== */
00452 
00453 void grism_backward(const Grism *grism, double lambda, 
00454                    double *x, double *y, double *z, int order)
00455 {
00456   double n;
00457   
00458   n = sellmeier_index(lambda, grism->prism.sellcoeffs);
00459 
00460   /* Get aligned with (A-tilted) grating */
00461   x_rotation(x,y,z, grism->tiltx);
00462   y_rotation(x,y,z, grism->tilty);
00463   z_rotation(x,y,z, M_PI+grism->rotz);
00464   x_rotation(x,y,z, grism->A);
00465 
00466   /* Grating */
00467   grating_backward(x,y,z, lambda,n,grism->g_per_mm,order);
00468 
00469   /* Prism */
00470   x_rotation(x,y,z, -grism->A);
00471 
00472   /* Glass/air interface */
00473   refraction(x,y,z, n,1.);
00474 
00475   /* Back to original orientation */
00476   z_rotation(x,y,z, -(M_PI+grism->rotz));
00477   y_rotation(x,y,z, -grism->tilty);
00478   x_rotation(x,y,z, -grism->tiltx);
00479 }
00480 
00481 /* === Doxygen Comment ======================================= */
00501 /* =========================================================== */
00502 
00503 void snifs_optics_forward(const SnifsOptics *optics, 
00504                           double xmla,double ymla, double *xccd,double *yccd, 
00505                           double lambda, int order)
00506 {
00507   double rmla, rccd, x,y,z, gamma;
00508   double phimla, phi, phiccd, tantheta, costheta, sintheta, tmp;
00509   
00510   /* Spectrograph magnification (fcam/fcoll) */
00511   gamma = optics->camera.focal / optics->collimator.focal;
00512 
00513   /* Initial position: (xmla,ymla) -> (rmla,phimla) */
00514   rmla = hypot(xmla,ymla);
00515   phimla = atan2(ymla,xmla); /* polar angle wrt x-axis in [-Pi,Pi] */
00516 
00517   /* Collimator: r>0 -> tan(theta)>0 */
00518   tantheta = collimator_forward(&(optics->collimator), gamma, lambda, rmla);
00519   phi = phimla + M_PI;
00520 
00521   /* Direction (normed) vector: (0<theta<Pi/2,phi) -> (x,y,z) */
00522   tmp = SQ(tantheta);
00523   costheta = sqrt(  1./(1. + tmp) );
00524   sintheta = sqrt( tmp/(1. + tmp) );
00525   x = cos(phi)*sintheta;
00526   y = sin(phi)*sintheta;
00527   z = costheta;
00528 
00529   /* Grism: (x,y,z) -> (x,y,z) */
00530   grism_forward(&(optics->grism), lambda, &x,&y,&z, order);
00531 
00532   /* Inclination: (x,y,z) -> (0<theta<Pi/2,phi) */
00533   if (z<=0.) {
00534     print_error("ARGH! cannot get out of the grism");
00535     exit_session(UNKNOWN);
00536   }
00537   tantheta = hypot(x,y)/z;
00538   phi = atan2(y,x);
00539 
00540   /* Camera: tan(theta)>0 -> r>0 */
00541   rccd = camera_forward(&(optics->camera), lambda, tantheta);
00542   phiccd = phi + M_PI;
00543 
00544   /* Final position: (rccd,phiccd) -> (xccd,yccd) */
00545   *xccd = rccd*cos(phiccd);
00546   *yccd = rccd*sin(phiccd);
00547 }
00548 
00549 /* === Doxygen Comment ======================================= */
00568 /* =========================================================== */
00569 
00570 void snifs_optics_backward(const SnifsOptics *optics, 
00571                            double xccd,double yccd, double *xmla,double *ymla, 
00572                            double lambda, int order)
00573 {
00574   double rmla, rccd, x,y,z, gamma;
00575   double phimla, phi, phiccd, tantheta, costheta, sintheta, tmp;
00576   
00577   /* Spectrograph magnification (fcam/fcoll) */
00578   gamma = optics->camera.focal / optics->collimator.focal;
00579 
00580   /* Final position: (xccd,yccd) -> (rccd,phiccd) */
00581   rccd = hypot(xccd,yccd);
00582   phiccd = atan2(yccd,xccd); /* polar angle wrt x-axis in [-Pi,Pi] */
00583 
00584   /* Camera backward: tan(theta)>0 <- r>0 */
00585   tantheta = camera_backward(&(optics->camera), lambda, rccd);
00586   phi = phiccd - M_PI;
00587 
00588   /* Direction (normed) vector: (0<theta<Pi/2,phi) -> (x,y,z) */
00589   tmp = SQ(tantheta);
00590   costheta = sqrt(  1./(1. + tmp) );
00591   sintheta = sqrt( tmp/(1. + tmp) );
00592   x = cos(phi)*sintheta;
00593   y = sin(phi)*sintheta;
00594   z = costheta;
00595 
00596   /* Grism backward: (x,y,z) -> (x,y,z) */
00597   grism_backward(&(optics->grism), lambda, &x,&y,&z, order);
00598 
00599   /* Inclination: (x,y,z) -> (0<theta<Pi/2,phi) */
00600   if (z<=0.) {
00601     print_error("ARGH! cannot get in the grism");
00602     exit_session(UNKNOWN);
00603   }
00604   tantheta = hypot(x,y)/z;
00605   phi = atan2(y,x);
00606 
00607   /* Collimator backward: r>0 <- tan(theta)>0 */
00608   rmla = collimator_backward(&(optics->collimator), gamma, lambda, tantheta);
00609   phimla = phi - M_PI;
00610 
00611   /* Initial position: (rmla,phimla) -> (xmla,ymla) */
00612   *xmla = rmla*cos(phimla);
00613   *ymla = rmla*sin(phimla);
00614 }
00615 
00616 /* === Doxygen Comment ======================================= */
00620 /* =========================================================== */
00621 
00622 void snifs_optics_test_forward(const SnifsOptics *optics, 
00623                                double xmla,double ymla, 
00624                                double *xccd,double *yccd, 
00625                                double lambda, int order)
00626 {
00627   double rmla, rccd, x,y,z, xtmp,ytmp,ztmp, gamma;
00628   double phimla, phi, phiccd, tantheta, costheta, sintheta, tmp;
00629 
00630   print_msg("** Testing SNIFS optics (O%d) forward for lambda=%.1f AA",
00631             order,lambda);
00632   
00633   /* Spectrograph magnification (fcam/fcoll) */
00634   gamma = optics->camera.focal / optics->collimator.focal;
00635 
00636   print_msg("   Initial position: xmla,ymla = %g mm, %g mm",xmla,ymla);
00637 
00638   /* Initial position: (xmla,ymla) -> (rmla,phimla) */
00639   rmla = hypot(xmla,ymla);
00640   phimla = atan2(ymla,xmla); /* polar angle wrt x-axis in [-Pi,Pi] */
00641 
00642   print_msg("   Initial position: rmla,phimla = %g mm, %g deg",
00643             rmla,phimla/M_PI*180.);
00644 
00645   /* Collimator: r>0 -> tan(theta)>0 */
00646   tantheta = collimator_forward(&(optics->collimator), gamma, lambda, rmla);
00647   phi = phimla + M_PI;
00648 
00649   print_msg("   Collimator: tan(theta),phi = %g, %g deg",
00650             tantheta,phi/M_PI*180.);
00651 
00652   /* Collimator backward: r>0 <- tan(theta)>0 */
00653   rmla = collimator_backward(&(optics->collimator), gamma, lambda, tantheta);
00654   phimla = phi - M_PI;
00655 
00656   print_msg("   Collimator backward: rmla,phimla = %g mm, %g deg",
00657             rmla,phimla/M_PI*180.);
00658 
00659   /* Direction (normed) vector: (0<theta<Pi/2,phi) -> (x,y,z) */
00660   tmp = SQ(tantheta);
00661   costheta = sqrt(  1./(1. + tmp) );
00662   sintheta = sqrt( tmp/(1. + tmp) );
00663   x = cos(phi)*sintheta;
00664   y = sin(phi)*sintheta;
00665   z = costheta;
00666 
00667   print_msg("   Direction: x,y,z = %g, %g, %g",x,y,z);
00668 
00669   /* Grism: (x,y,z) -> (x,y,z) */
00670 #ifdef MEGA_VERBOSE
00671   if (0) {
00672     /* detail the grism stage */
00673     Grism *grism=&(optics->grism);
00674     double n;
00675     
00676     n = sellmeier_index(lambda, grism->prism.sellcoeffs);
00677     print_msg("   Refractive index: n(%.1f AA) = %g",lambda,n);
00678 
00679     /* Get aligned w/ grism */
00680     x_rotation(&x,&y,&z, grism->tiltx);
00681     y_rotation(&x,&y,&z, grism->tilty);
00682     z_rotation(&x,&y,&z, M_PI+grism->rotz);
00683     print_msg("   Grism rotations (%g,%g,%g deg): "
00684               "x,y,z = %g, %g, %g",
00685               grism->tiltx/M_PI*180.,grism->tilty/M_PI*180.,
00686               grism->rotz/M_PI*180.,x,y,z);
00687 
00688     /* Air/glass interface */
00689     refraction(&x,&y,&z, 1.,n);
00690     print_msg("   Grism refraction: x,y,z = %g, %g, %g",x,y,z);
00691 
00692     /* Prism */
00693     x_rotation(&x,&y,&z, grism->A);
00694     print_msg("   Prism (%g deg): x,y,z = %g, %g, %g",
00695               grism->A/M_PI*180.,x,y,z);
00696 
00697     /* Grating */
00698     grating_forward(&x,&y,&z, lambda,n,grism->g_per_mm,order);
00699     print_msg("   Grating (%.1f g/mm, O%d): x,y,z = %g, %g, %g",
00700               grism->g_per_mm,order,x,y,z);
00701 
00702     /* Back to original orientation */
00703     x_rotation(&x,&y,&z, -grism->A);
00704     print_msg("   Prism back-rot: x,y,z = %g, %g, %g",x,y,z);
00705     z_rotation(&x,&y,&z, -(M_PI+grism->rotz));
00706     y_rotation(&x,&y,&z, -grism->tilty);
00707     x_rotation(&x,&y,&z, -grism->tiltx);
00708     print_msg("   Grism back-rot: x,y,z = %g, %g, %g",x,y,z);
00709   }
00710 #endif
00711 
00712   grism_forward(&(optics->grism), lambda, &x,&y,&z, order);
00713   
00714   print_msg("   Grism: x,y,z = %g, %g, %g",x,y,z);
00715 
00716   /* Grism backward: (x,y,z) -> (x,y,z) */
00717   xtmp = x;
00718   ytmp = y;
00719   ztmp = z;
00720   grism_backward(&(optics->grism), lambda, &xtmp,&ytmp,&ztmp, order);
00721 
00722   print_msg("   Grism backward: x,y,z = %g, %g, %g",xtmp,ytmp,ztmp);
00723 
00724  /* Inclination: (x,y,z) -> (0<theta<Pi/2,phi) */
00725   if (z<=0.) {
00726     print_error("ARGH! cannot get out of the grism");
00727     exit_session(UNKNOWN);
00728   }
00729   tantheta = hypot(x,y)/z;
00730   phi = atan2(y,x);
00731 
00732   print_msg("   Direction: tan(theta),phi = %g, %g deg",
00733             tantheta,phi/M_PI*180.);
00734 
00735   /* Camera: tan(theta)>0 -> r>0 */
00736   rccd = camera_forward(&(optics->camera), lambda, tantheta);
00737   phiccd = phi + M_PI;
00738 
00739   print_msg("   Camera: rccd,phiccd = %g mm, %g deg",
00740             rccd,phiccd/M_PI*180.);
00741 
00742   /* Camera backward: tan(theta)>0 <- r>0 */
00743   tantheta = camera_backward(&(optics->camera), lambda, rccd);
00744   phi = phiccd - M_PI;
00745 
00746   print_msg("   Camera backward: tan(theta),phi = %g, %g deg",
00747             tantheta,phi/M_PI*180.);
00748 
00749   /* Final position: (rccd,phiccd) -> (xccd,yccd) */
00750   *xccd = rccd*cos(phiccd);
00751   *yccd = rccd*sin(phiccd);
00752 
00753   print_msg("   Final position: xccd,yccd = %g mm, %g mm",*xccd,*yccd);
00754 }
00755 
00756 /* === Doxygen Comment ======================================= */
00768 /* =========================================================== */
00769 
00770 void snifs_optics_CCD2MLA(const SnifsOptics *optics, 
00771                           double xccd,double yccd, double *xmla,double *ymla, 
00772                           double pixsize,double lbdaref)
00773 { 
00774   double xpup,ypup;
00775   
00776   /* [CCD/px] -> [CCD/mm] */
00777   xpup = (xccd - optics->center.xc)*pixsize;
00778   ypup = (yccd - optics->center.yc)*pixsize;
00779 
00780   /* [CCD/mm] -> [MLA/mm] */
00781   snifs_optics_backward(optics, xpup,ypup, xmla,ymla, lbdaref,1);
00782 }
00783 
00784 /* === Doxygen Comment ======================================= */
00797 /* =========================================================== */
00798 
00799 void snifs_optics_MLA2CCD(const SnifsOptics *optics, 
00800                           double xmla,double ymla, double *xccd,double *yccd, 
00801                           double pixsize,double lbda, int order)
00802 { 
00803   /* [MLA/mm] -> [CCD/mm] */
00804   snifs_optics_forward(optics, xmla,ymla, xccd,yccd,lbda,order);
00805         
00806   /* [CCD/mm] -> [CCD/px] */
00807   *xccd = (*xccd)/pixsize + optics->center.xc;
00808   *yccd = (*yccd)/pixsize + optics->center.yc;
00809 }
00810 
00811 /* 
00812 ;;; Local Variables: ***
00813 ;;; eval: (add-to-list 'c-font-lock-extra-types "Grism") ***
00814 ;;; eval: (add-to-list 'c-font-lock-extra-types "CameraCollimator") ***
00815 ;;; eval: (add-to-list 'c-font-lock-extra-types "SnifsOptics") ***
00816 ;;; End: ***
00817 */

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