LevelS C++ support library  3.82
trafos.cc
Go to the documentation of this file.
1 /*
2  * This file is part of libcxxsupport.
3  *
4  * libcxxsupport is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * libcxxsupport is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with libcxxsupport; if not, write to the Free Software
16  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17  */
18 
19 /*
20  * libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
21  * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22  * (DLR).
23  */
24 
25 /*! \file trafos.cc
26  * Celestial coordinate transformations.
27  *
28  * Copyright (C) 2005-2011 Max-Planck-Society
29  * \author Martin Reinecke
30  */
31 
32 #include "trafos.h"
33 #include "geom_utils.h"
34 #include "lsconstants.h"
35 
36 using namespace std;
37 
38 vec3 Trafo::xcc_dp_precess (const vec3 &iv, double iepoch,
39  double oepoch)
40  {
41  // Z-axis rotation by OBL_LONG:
42  double Tm = ((oepoch+iepoch)*0.5 - 1900.) *0.01;
43  double gp_long = degr2rad * ((oepoch-iepoch) * (50.2564+0.0222*Tm) / 3600.);
44  double obl_long =
45  degr2rad * (180. - (173. + (57.06+54.77*Tm) / 60.)) + gp_long*0.5;
46  double dco = cos(obl_long), dso = sin(obl_long);
47  vec3 ov (iv.x*dco-iv.y*dso, iv.x*dso+iv.y*dco, iv.z);
48 
49  // X-axis rotation by dE:
50  double dE = degr2rad * ((oepoch-iepoch) * (0.4711-0.0007*Tm) / 3600.);
51  double dce = cos(dE), dse = sin(dE);
52  double temp = ov.y*dce - ov.z*dse;
53  ov.z = ov.y*dse + ov.z*dce;
54  ov.y = temp;
55 
56  // Z-axis rotation by GP_LONG - OBL_LONG:
57  double dL = gp_long - obl_long;
58  double dcl = cos(dL), dsl = sin(dL);
59  temp = ov.x*dcl - ov.y*dsl;
60  ov.y = ov.x*dsl + ov.y*dcl;
61  ov.x = temp;
62 
63  return ov;
64  }
65 
66 double Trafo::get_epsilon (double epoch)
67  {
68  double T = (epoch - 1900.) * 0.01;
69  return
70  degr2rad * (23.452294 - 0.0130125*T - 1.63889e-6*T*T + 5.02778e-7*T*T*T);
71  }
72 
73 /*! Routine to convert from ecliptic coordinates to equatorial (celestial)
74  coordinates at the given epoch. Adapted from the COBLIB routine by Dr.
75  E. Wright. */
76 vec3 Trafo::xcc_dp_e_to_q (const vec3 &iv, double epoch)
77  {
78  double epsilon=get_epsilon(epoch);
79  double dc = cos(epsilon), ds = sin(epsilon);
80  return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
81  }
82 
83 vec3 Trafo::xcc_dp_q_to_e (const vec3 &iv, double epoch)
84  {
85  double epsilon=-get_epsilon(epoch);
86  double dc = cos(epsilon), ds = sin(epsilon);
87  return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
88  }
89 
90 /*! Routine to convert galactic coordinates to ecliptic (celestial)
91  coordinates at the given epoch. First the conversion to ecliptic
92  2000 is done, then if necessary the results are precessed. */
93 vec3 Trafo::xcc_dp_g_to_e (const vec3 &iv, double epoch)
94  {
95  static const rotmatrix T (-0.054882486, 0.494116468, -0.867661702,
96  -0.993821033, -0.110993846, -0.000346354,
97  -0.096476249, 0.862281440, 0.497154957);
98  vec3 hv=T.Transform(iv);
99 
100  if (!approx(epoch,2000.))
101  hv=xcc_dp_precess(hv,2000.,epoch);
102 
103  return hv;
104  }
105 
106 /*! Routine to convert ecliptic (celestial) coordinates at the given
107  epoch to galactic coordinates. The ecliptic coordinates are first
108  precessed to 2000.0, then converted. */
109 vec3 Trafo::xcc_dp_e_to_g (const vec3 &iv, double epoch)
110  {
111  static const rotmatrix T (-0.054882486, -0.993821033, -0.096476249,
112  0.494116468, -0.110993846, 0.862281440,
113  -0.867661702, -0.000346354, 0.497154957);
114  vec3 hv=iv;
115  if (!approx(epoch,2000.))
116  hv=xcc_dp_precess(hv,epoch,2000.);
117 
118  return T.Transform(hv);
119  }
120 
121 /*! Function to convert between standard coordinate systems, including
122  precession. */
123 vec3 Trafo::xcc_v_convert(const vec3 &iv, double iepoch, double oepoch,
124  coordsys isys,coordsys osys)
125  {
126  vec3 xv;
127  if (isys == Ecliptic)
128  xv=iv;
129  else if (isys == Equatorial)
130  xv = xcc_dp_q_to_e(iv,iepoch);
131  else if (isys == Galactic)
132  xv = xcc_dp_g_to_e(iv,iepoch);
133  else
134  planck_fail("Unsupported input coordinate system");
135 
136  vec3 yv = approx(iepoch,oepoch) ? xv : xcc_dp_precess(xv,iepoch,oepoch);
137 
138  vec3 ov;
139  if (osys == Ecliptic)
140  ov = yv;
141  else if (osys == Equatorial)
142  ov = xcc_dp_e_to_q(yv,oepoch);
143  else if (osys == Galactic)
144  ov = xcc_dp_e_to_g(yv,oepoch);
145  else
146  planck_fail("Unsupported output coordinate system");
147 
148  return ov;
149  }
150 
151 void Trafo::coordsys2matrix (double iepoch, double oepoch,
152  coordsys isys, coordsys osys, rotmatrix &matrix)
153  {
154  vec3 v1p = xcc_v_convert(vec3(1,0,0),iepoch,oepoch,isys,osys),
155  v2p = xcc_v_convert(vec3(0,1,0),iepoch,oepoch,isys,osys),
156  v3p = xcc_v_convert(vec3(0,0,1),iepoch,oepoch,isys,osys);
157  v1p.Normalize(); v2p.Normalize(); v3p.Normalize();
158  matrix=rotmatrix(v1p,v2p,v3p);
159  }
160 
161 Trafo::Trafo (double iepoch, double oepoch, coordsys isys, coordsys osys)
162  { coordsys2matrix (iepoch, oepoch, isys, osys, mat); }
163 
164 pointing Trafo::operator() (const pointing &ptg) const
165  { return pointing(operator()(vec3(ptg))); }
166 
167 void Trafo::rotatefull (const pointing &ptg, pointing &newptg,
168  double &delta_psi) const
169  {
170  vec3 vec (ptg);
171  vec3 east (-vec.y,vec.x,0.);
172  vec3 newvec = operator()(vec);
173  vec3 neweast = operator()(east);
174  delta_psi = orientation(newvec,neweast)+halfpi;
175  newptg = newvec;
176  }
177 
178 void Trafo::rotatefull (pointing &ptg, double &psi) const
179  {
180  vec3 vec (ptg);
181  vec3 east (-vec.y,vec.x,0.);
182  vec3 newvec = operator()(vec);
183  vec3 neweast = operator()(east);
184  psi += orientation(newvec,neweast)+halfpi;
185  ptg = newvec;
186  }
187 
188 void Trafo::rotatefull (const vec3 &vec, vec3 &newvec, double &delta_psi) const
189  {
190  vec3 east (-vec.y,vec.x,0.);
191  newvec = operator()(vec);
192  vec3 neweast = operator()(east);
193  delta_psi = orientation(newvec,neweast)+halfpi;
194  }
195 
196 void Trafo::rotatefull (vec3 &vec, double &psi) const
197  {
198  vec3 east (-vec.y,vec.x,0.);
199  vec3 newvec = operator()(vec);
200  vec3 neweast = operator()(east);
201  psi += orientation(newvec,neweast)+halfpi;
202  vec = newvec;
203  }
vec3_t::y
T y
Definition: vec3.h:46
geom_utils.h
vec3_t::z
T z
Definition: vec3.h:46
Trafo::Trafo
Trafo(double iepoch, double oepoch, coordsys isys, coordsys osys)
Definition: trafos.cc:161
Trafo::operator()
vec3 operator()(const vec3 &vec) const
Definition: trafos.h:63
vec3_t::x
T x
Definition: vec3.h:46
vec3
vec3_t< float64 > vec3
Definition: vec3.h:143
orientation
double orientation(const vec3 &loc, const vec3 &dir)
Definition: geom_utils.h:47
approx
bool approx(F a, F b, F epsilon=1e-5)
Definition: math_utils.h:44
vec3_t
Definition: vec3.h:43
planck_fail
#define planck_fail(msg)
Definition: error_handling.h:69
vec3_t::Normalize
void Normalize()
Definition: vec3.h:72
Trafo::rotatefull
void rotatefull(const pointing &ptg, pointing &newptg, double &delta_psi) const
Definition: trafos.cc:167
Generated on Thu Jul 28 2022 17:32:06 for LevelS C++ support library