File: trafos.cc

package info (click to toggle)
healpix-cxx 3.50.0-3
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 3,008 kB
  • sloc: cpp: 13,882; ansic: 5,930; sh: 4,451; makefile: 216
file content (203 lines) | stat: -rw-r--r-- 6,384 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
/*
 *  This file is part of libcxxsupport.
 *
 *  libcxxsupport is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  libcxxsupport is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libcxxsupport; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 *  libcxxsupport is being developed at the Max-Planck-Institut fuer Astrophysik
 *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
 *  (DLR).
 */

/*! \file trafos.cc
 *  Celestial coordinate transformations.
 *
 *  Copyright (C) 2005-2011 Max-Planck-Society
 * \author Martin Reinecke
 */

#include "trafos.h"
#include "geom_utils.h"
#include "lsconstants.h"

using namespace std;

vec3 Trafo::xcc_dp_precess (const vec3 &iv, double iepoch,
  double oepoch)
  {
  // Z-axis rotation by OBL_LONG:
  double Tm = ((oepoch+iepoch)*0.5 - 1900.) *0.01;
  double gp_long  = degr2rad * ((oepoch-iepoch) * (50.2564+0.0222*Tm) / 3600.);
  double obl_long =
    degr2rad * (180. - (173. + (57.06+54.77*Tm) / 60.)) + gp_long*0.5;
  double dco = cos(obl_long), dso = sin(obl_long);
  vec3 ov (iv.x*dco-iv.y*dso, iv.x*dso+iv.y*dco, iv.z);

  // X-axis rotation by dE:
  double dE = degr2rad * ((oepoch-iepoch) * (0.4711-0.0007*Tm) / 3600.);
  double dce = cos(dE), dse = sin(dE);
  double temp = ov.y*dce - ov.z*dse;
  ov.z = ov.y*dse + ov.z*dce;
  ov.y = temp;

  // Z-axis rotation by GP_LONG - OBL_LONG:
  double dL = gp_long - obl_long;
  double dcl = cos(dL), dsl = sin(dL);
  temp = ov.x*dcl - ov.y*dsl;
  ov.y = ov.x*dsl + ov.y*dcl;
  ov.x = temp;

  return ov;
  }

double Trafo::get_epsilon (double epoch)
  {
  double T = (epoch - 1900.) * 0.01; 
  return
    degr2rad * (23.452294 - 0.0130125*T - 1.63889e-6*T*T + 5.02778e-7*T*T*T);
  }

/*! Routine to convert from ecliptic coordinates to equatorial (celestial)
    coordinates at the given epoch.  Adapted from the COBLIB routine by Dr.
    E. Wright. */
vec3 Trafo::xcc_dp_e_to_q (const vec3 &iv, double epoch)
  {
  double epsilon=get_epsilon(epoch);
  double dc = cos(epsilon), ds = sin(epsilon);
  return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
  }

vec3 Trafo::xcc_dp_q_to_e (const vec3 &iv, double epoch)
  {
  double epsilon=-get_epsilon(epoch);
  double dc = cos(epsilon), ds = sin(epsilon);
  return vec3 (iv.x, dc*iv.y-ds*iv.z, dc*iv.z+ds*iv.y);
  }

/*! Routine to convert galactic coordinates to ecliptic (celestial) 
    coordinates at the given epoch.  First the conversion to ecliptic 
    2000 is done, then if necessary the results are precessed. */
vec3 Trafo::xcc_dp_g_to_e (const vec3 &iv, double epoch)
  {
  static const rotmatrix T (-0.054882486,  0.494116468, -0.867661702,
                            -0.993821033, -0.110993846, -0.000346354,
                            -0.096476249,  0.862281440,  0.497154957);
  vec3 hv=T.Transform(iv);

  if (!approx(epoch,2000.))
    hv=xcc_dp_precess(hv,2000.,epoch);

  return hv;
  }

/*! Routine to convert ecliptic (celestial) coordinates at the given
    epoch to galactic coordinates.  The ecliptic coordinates are first
    precessed to 2000.0, then converted. */
vec3 Trafo::xcc_dp_e_to_g (const vec3 &iv, double epoch)
  {
  static const rotmatrix T (-0.054882486, -0.993821033, -0.096476249,
                             0.494116468, -0.110993846,  0.862281440,
                            -0.867661702, -0.000346354,  0.497154957);
  vec3 hv=iv;
  if (!approx(epoch,2000.))
    hv=xcc_dp_precess(hv,epoch,2000.);

  return T.Transform(hv);
  }

/*! Function to convert between standard coordinate systems, including
    precession. */
vec3 Trafo::xcc_v_convert(const vec3 &iv, double iepoch, double oepoch,
  coordsys isys,coordsys osys)
  {
  vec3 xv;
  if (isys == Ecliptic)
    xv=iv;
  else if (isys == Equatorial)
    xv = xcc_dp_q_to_e(iv,iepoch);
  else if (isys == Galactic)
    xv = xcc_dp_g_to_e(iv,iepoch);
  else
    planck_fail("Unsupported input coordinate system");

  vec3 yv = approx(iepoch,oepoch) ? xv : xcc_dp_precess(xv,iepoch,oepoch);

  vec3 ov;
  if (osys == Ecliptic)
    ov = yv;
  else if (osys == Equatorial)
    ov = xcc_dp_e_to_q(yv,oepoch);
  else if (osys == Galactic)
    ov = xcc_dp_e_to_g(yv,oepoch);
  else
    planck_fail("Unsupported output coordinate system");

  return ov;
  }

void Trafo::coordsys2matrix (double iepoch, double oepoch,
  coordsys isys, coordsys osys, rotmatrix &matrix)
  {
  vec3 v1p = xcc_v_convert(vec3(1,0,0),iepoch,oepoch,isys,osys),
       v2p = xcc_v_convert(vec3(0,1,0),iepoch,oepoch,isys,osys),
       v3p = xcc_v_convert(vec3(0,0,1),iepoch,oepoch,isys,osys);
  v1p.Normalize(); v2p.Normalize(); v3p.Normalize();
  matrix=rotmatrix(v1p,v2p,v3p);
  }

Trafo::Trafo (double iepoch, double oepoch, coordsys isys, coordsys osys)
  { coordsys2matrix (iepoch, oepoch, isys, osys, mat); }

pointing Trafo::operator() (const pointing &ptg) const
  { return pointing(operator()(vec3(ptg))); }

void Trafo::rotatefull (const pointing &ptg, pointing &newptg,
  double &delta_psi) const
  {
  vec3 vec (ptg);
  vec3 east (-vec.y,vec.x,0.);
  vec3 newvec = operator()(vec);
  vec3 neweast = operator()(east);
  delta_psi = orientation(newvec,neweast)+halfpi;
  newptg = newvec;
  }

void Trafo::rotatefull (pointing &ptg, double &psi) const
  {
  vec3 vec (ptg);
  vec3 east (-vec.y,vec.x,0.);
  vec3 newvec = operator()(vec);
  vec3 neweast = operator()(east);
  psi += orientation(newvec,neweast)+halfpi;
  ptg = newvec;
  }

void Trafo::rotatefull (const vec3 &vec, vec3 &newvec, double &delta_psi) const
  {
  vec3 east (-vec.y,vec.x,0.);
  newvec = operator()(vec);
  vec3 neweast = operator()(east);
  delta_psi = orientation(newvec,neweast)+halfpi;
  }

void Trafo::rotatefull (vec3 &vec, double &psi) const
  {
  vec3 east (-vec.y,vec.x,0.);
  vec3 newvec = operator()(vec);
  vec3 neweast = operator()(east);
  psi += orientation(newvec,neweast)+halfpi;
  vec = newvec;
  }