trafos.cc

Go to the documentation of this file.
00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*! \file trafos.cc
00028  *  Celestial coordinate transformations.
00029  *
00030  *  Copyright (C) 2005 Max-Planck-Society
00031  * \author Martin Reinecke
00032  */
00033 
00034 #include "trafos.h"
00035 #include "lsconstants.h"
00036 
00037 using namespace std;
00038 
00039 vec3 Trafo::xcc_dp_precess (const vec3 &iv, double iepoch,
00040   double oepoch)
00041   {
00042   // Z-axis rotation by OBL_LONG:
00043   double Tm = ((oepoch+iepoch)*0.5 - 1900.) *0.01;
00044   double gp_long  = (oepoch-iepoch) * (50.2564+0.0222*Tm) / 3600.;
00045   double obl_long = 180. - (173. + (57.06+54.77*Tm) / 60.) + gp_long*0.5;
00046   double dco = cos(obl_long*degr2rad), dso = sin(obl_long*degr2rad);
00047   vec3 ov (iv.x*dco-iv.y*dso, iv.x*dso+iv.y*dco, iv.z);
00048 
00049   // X-axis rotation by dE:
00050   double dE = (oepoch-iepoch) * (0.4711-0.0007*Tm) / 3600.;
00051   double dce = cos(dE*degr2rad), dse = sin(dE*degr2rad);
00052   double temp = ov.y*dce - ov.z*dse;
00053   ov.z = ov.y*dse + ov.z*dce;
00054   ov.y = temp;
00055 
00056   // Z-axis rotation by GP_LONG - OBL_LONG:
00057   double dL = gp_long - obl_long;
00058   double dcl = cos(dL*degr2rad), dsl = sin(dL*degr2rad);
00059   temp = ov.x*dcl - ov.y*dsl;
00060   ov.y = ov.x*dsl + ov.y*dcl;
00061   ov.x = temp;
00062 
00063   return ov;
00064   }
00065 
00066 double Trafo::get_epsilon (double epoch)
00067   {
00068   double T = (epoch - 1900.) * 0.01; 
00069   double epsilon = 23.452294 - 0.0130125*T - 1.63889e-6*T*T + 5.02778e-7*T*T*T;
00070   return epsilon*degr2rad;
00071   }
00072 
00073 /*! Routine to convert from ecliptic coordinates to equatorial (celestial)
00074     coordinates at the given epoch.  Adapted from the COBLIB routine by Dr.
00075     E. Wright. */
00076 vec3 Trafo::xcc_dp_e_to_q (const vec3 &iv, double epoch)
00077   {
00078   double epsilon=get_epsilon(epoch);
00079   double dc = cos(epsilon), ds = sin(epsilon);
00080   return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
00081   }
00082 
00083 vec3 Trafo::xcc_dp_q_to_e (const vec3 &iv, double epoch)
00084   {
00085   double epsilon=-get_epsilon(epoch);
00086   double dc = cos(epsilon), ds = sin(epsilon);
00087   return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
00088   }
00089 
00090 /*! Routine to convert galactic coordinates to ecliptic (celestial) 
00091     coordinates at the given epoch.  First the conversion to ecliptic 
00092     2000 is done, then if necessary the results are precessed. */
00093 vec3 Trafo::xcc_dp_g_to_e (const vec3 &iv, double epoch)
00094   {
00095   static rotmatrix T (-0.054882486,  0.494116468, -0.867661702,
00096                       -0.993821033, -0.110993846, -0.000346354,
00097                       -0.096476249,  0.862281440,  0.497154957);
00098   vec3 hv=T.Transform(iv);
00099 
00100   if (!approx(epoch,2000.))
00101     hv=xcc_dp_precess(hv,2000.,epoch);
00102 
00103   return hv;
00104   }
00105 
00106 /*! Routine to convert ecliptic (celestial) coordinates at the given
00107     epoch to galactic coordinates.  The ecliptic coordinates are first
00108     precessed to 2000.0, then converted. */
00109 vec3 Trafo::xcc_dp_e_to_g (const vec3 &iv, double epoch)
00110   {
00111   static rotmatrix T (-0.054882486, -0.993821033, -0.096476249,
00112                        0.494116468, -0.110993846,  0.862281440,
00113                       -0.867661702, -0.000346354,  0.497154957);
00114   vec3 hv=iv;
00115   if (!approx(epoch,2000.))
00116     hv=xcc_dp_precess(hv,epoch,2000.);
00117 
00118   return T.Transform(hv);
00119   }
00120 
00121 /*! Function to convert between standard coordinate systems, including
00122     precession. */
00123 vec3 Trafo::xcc_v_convert(const vec3 &iv, double iepoch, double oepoch,
00124   coordsys isys,coordsys osys)
00125   {
00126   vec3 xv;
00127   if (isys == Ecliptic)
00128     xv=iv;
00129   else if (isys == Equatorial)
00130     xv = xcc_dp_q_to_e(iv,iepoch);
00131   else if (isys == Galactic)
00132     xv = xcc_dp_g_to_e(iv,iepoch);
00133   else
00134     throw Message_error("Unknown input coordinate system");
00135 
00136   vec3 yv = approx(iepoch,oepoch) ? xv : xcc_dp_precess(xv,iepoch,oepoch);
00137 
00138   vec3 ov;
00139   if (osys == Ecliptic)
00140     ov = yv;
00141   else if (osys == Equatorial)
00142     ov = xcc_dp_e_to_q(yv,oepoch);
00143   else if (osys == Galactic)
00144     ov = xcc_dp_e_to_g(yv,oepoch);
00145   else
00146     throw Message_error("Unknown output coordinate system");
00147 
00148   return ov;
00149   }
00150 
00151 void Trafo::coordsys2matrix (double iepoch, double oepoch,
00152   coordsys isys, coordsys osys, rotmatrix &matrix)
00153   {
00154   vec3 v1p = xcc_v_convert(vec3(1,0,0),iepoch,oepoch,isys,osys),
00155        v2p = xcc_v_convert(vec3(0,1,0),iepoch,oepoch,isys,osys),
00156        v3p = xcc_v_convert(vec3(0,0,1),iepoch,oepoch,isys,osys);
00157   v1p.Normalize(); v2p.Normalize(); v3p.Normalize();
00158   matrix=rotmatrix(v1p,v2p,v3p);
00159   }