LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Coord.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
34 #include <cmath>
35 #include <limits>
36 #include <cstdio>
37 #include <iomanip>
38 
39 #include "Eigen/Core"
40 #include "Eigen/Geometry"
41 
42 #include "lsst/pex/exceptions.h"
43 #include "boost/algorithm/string.hpp"
44 #include "boost/tuple/tuple.hpp"
45 #include "boost/format.hpp"
46 
47 #include "lsst/afw/coord/Coord.h"
48 #include "lsst/daf/base/DateTime.h"
49 
50 namespace afwCoord = lsst::afw::coord;
51 namespace ex = lsst::pex::exceptions;
52 namespace afwGeom = lsst::afw::geom;
53 namespace dafBase = lsst::daf::base;
54 
55 namespace {
56 
57 typedef std::map<std::string, afwCoord::CoordSystem> CoordSystemMap;
58 
59 CoordSystemMap const getCoordSystemMap() {
60  CoordSystemMap idMap;
61  idMap["FK5"] = afwCoord::FK5;
62  idMap["ICRS"] = afwCoord::ICRS;
63  idMap["ECLIPTIC"] = afwCoord::ECLIPTIC;
64  idMap["GALACTIC"] = afwCoord::GALACTIC;
65  idMap["ELON"] = afwCoord::ECLIPTIC;
66  idMap["GLON"] = afwCoord::GALACTIC;
67  idMap["TOPOCENTRIC"] = afwCoord::TOPOCENTRIC;
68  return idMap;
69 }
70 
75 afwGeom::Angle haversine(afwGeom::Angle const dAlpha,
76  afwGeom::Angle const dDelta,
77  double const cosDelta1,
78  double const cosDelta2
79  )
80 {
81  double const havDDelta = std::sin(dDelta/2.0) * std::sin(dDelta/2.0);
82  double const havDAlpha = std::sin(dAlpha/2.0) * std::sin(dAlpha/2.0);
83  double const havD = havDDelta + cosDelta1 * cosDelta2 * havDAlpha;
84  double const sinDHalf = std::sqrt(havD);
85  afwGeom::Angle dist = (2.0 * std::asin(sinDHalf)) * afwGeom::radians;
86  return dist;
87 }
88 
89 double const epochTolerance = 1.0e-12;
90 
91 // Put a pair of coordinates in a common FK5 system
92 std::pair<afwCoord::Fk5Coord, afwCoord::Fk5Coord> commonFk5(afwCoord::Coord const& c1,
93  afwCoord::Coord const& c2
94  )
95 {
96  // make sure they're fk5
97  afwCoord::Fk5Coord fk51 = c1.toFk5();
98  afwCoord::Fk5Coord fk5tmp = c2.toFk5();
99 
100  // make sure they have the same epoch
101  afwCoord::Fk5Coord fk52;
102  if (fabs(fk51.getEpoch() - fk5tmp.getEpoch()) > epochTolerance) {
103  fk52 = fk5tmp.precess(fk51.getEpoch());
104  } else {
105  fk52 = fk5tmp;
106  }
107 
108  return std::make_pair(fk51, fk52);
109 }
110 
111 } // end anonymous namespace
112 
113 
118  static CoordSystemMap idmap = getCoordSystemMap();
119  if (idmap.find(system) != idmap.end()) {
120  return idmap[system];
121  } else {
122  throw LSST_EXCEPT(ex::InvalidParameterError, "System " + system + " not defined.");
123  }
124 }
125 
126 
127 namespace {
128 
129 double const NaN = std::numeric_limits<double>::quiet_NaN();
130 double const JD2000 = 2451544.50;
131 
132 
133 /*
134  * A local class to handle dd:mm:ss coordinates
135  *
136  * This class allows a decimal or dd:mm:ss coordinate to be
137  * disassembed into d, m, and s components.
138  * It's in an anonymous namespace, but there are public functions
139  * which perform the transformations directly:
140  *
141  * --> std::string dmsStr = degreesToDmsString(double deg);
142  * --> double deg = dmsStringToDegrees(std::string dms);
143  */
144 class Dms {
145 public:
146  Dms() {};
147 
148  // note that isSouth is needed to specify coords between dec = 0, and dec = -1
149  // otherwise, d = -0 gets carried as d = 0 ... need a way to specify it explicitly
150  Dms(int const d, int const m, double const s, bool const isSouth=false) {
151  sign = (d < 0 || isSouth) ? -1 : 1;
152  deg = std::abs(d);
153  min = m;
154  sec = s;
155  };
156  // unit could be "degrees" or "hours"
157  Dms(afwGeom::Angle const deg00, afwGeom::AngleUnit const unit = afwGeom::degrees) {
158  double deg0 = deg00.asAngularUnits(unit);
159  double const absVal = std::fabs(deg0);
160  sign = (deg0 >= 0) ? 1 : -1;
161  deg = static_cast<int>(std::floor(absVal));
162  min = static_cast<int>(std::floor((absVal - deg)*60.0));
163  sec = ((absVal - deg)*60.0 - min)*60.0;
164  }
165 
166  int deg;
167  int min;
168  double sec;
169  int sign;
170 };
171 
172 
173 
178 afwCoord::Coord const& GalacticPoleInFk5()
179 {
180  static afwCoord::Coord pole(192.85950*afwGeom::degrees, 27.12825*afwGeom::degrees, 2000.0); // C&O
181  return pole;
182 }
183 
184 afwCoord::Coord const& Fk5PoleInGalactic()
185 {
186  static afwCoord::Coord pole(122.93200 * afwGeom::degrees, 27.12825 * afwGeom::degrees, 2000.0); // C&O
187  return pole;
188 }
189 
194 afwGeom::Angle meanSiderealTimeGreenwich(
195  double const jd
196  ) {
197  double const T = (jd - 2451545.0)/36525.0;
198  return (280.46061837 + 360.98564736629*(jd - 2451545.0) + 0.000387933*T*T - (T*T*T/38710000.0)) * afwGeom::degrees;
199 }
200 
201 
202 /*
203  * A pair of utility functions to go from cartesian to spherical
204  */
205 double const atPoleEpsilon = 0.0; //std::numeric_limits<double>::epsilon();
206 afwGeom::Angle pointToLongitude(lsst::afw::geom::Point3D const &p3d, double const defaultLongitude=0.0) {
207  afwGeom::Angle lon;
208  if (fabs(p3d.getX()) <= atPoleEpsilon && fabs(p3d.getY()) <= atPoleEpsilon) {
209  lon = afwGeom::Angle(0.0);
210  } else {
211  lon = std::atan2(p3d.getY(), p3d.getX()) * afwGeom::radians;
212  lon.wrap();
213  }
214  return lon;
215 }
216 afwGeom::Angle pointToLatitude(lsst::afw::geom::Point3D const &p3d) {
217  afwGeom::Angle lat;
218  if ( fabs(p3d.getX()) <= atPoleEpsilon && fabs(p3d.getY()) <= atPoleEpsilon) {
219  lat = (p3d.getZ() >= 0) ? afwGeom::Angle(afwGeom::HALFPI) : afwGeom::Angle(-afwGeom::HALFPI);
220  } else {
221  lat = std::asin(p3d.getZ()) * afwGeom::radians;
222  }
223  return lat;
224 }
225 std::pair<afwGeom::Angle, afwGeom::Angle> pointToLonLat(lsst::afw::geom::Point3D const &p3d, double const defaultLongitude=0.0, bool normalize=true) {
226  std::pair<afwGeom::Angle, afwGeom::Angle> lonLat;
227 
228  if (normalize) {
229  double const inorm = 1.0/p3d.asEigen().norm();
230  double const x = inorm*p3d.getX();
231  double const y = inorm*p3d.getY();
232  double const z = inorm*p3d.getZ();
233  if (fabs(x) <= atPoleEpsilon && fabs(y) <= atPoleEpsilon) {
234  lonLat.first = 0.0 * afwGeom::radians;
235  lonLat.second = ((z >= 0) ? 1.0 : -1.0) * afwGeom::HALFPI * afwGeom::radians;
236  } else {
237  lonLat.first = atan2(y, x) * afwGeom::radians;
238  lonLat.first.wrap();
239  lonLat.second = asin(z) * afwGeom::radians;
240  }
241  } else {
242  if (fabs(p3d.getX()) <= atPoleEpsilon && fabs(p3d.getY()) <= atPoleEpsilon) {
243  lonLat.first = 0.0 * afwGeom::radians;
244  lonLat.second = ((p3d.getZ() >= 0) ? 1.0 : -1.0) * afwGeom::HALFPI * afwGeom::radians;
245  } else {
246  lonLat.first = atan2(p3d.getY(), p3d.getX()) * afwGeom::radians;
247  lonLat.first.wrap();
248  lonLat.second = asin(p3d.getZ()) * afwGeom::radians;
249  }
250  }
251  return lonLat;
252 }
253 
254 
255 } // end anonymous namespace
256 
257 
258 
259 /******************* Public functions ********************/
260 
261 static std::string angleToXmsString(afwGeom::Angle const a, afwGeom::AngleUnit const unit) {
262 
263  Dms dms(a, unit);
264 
265  // make sure rounding won't give 60.00 for sec or min
266  if ( (60.00 - dms.sec) < 0.005 ) {
267  dms.sec = 0.0;
268  dms.min += 1;
269  if (dms.min == 60) {
270  dms.min = 0;
271  dms.deg += 1;
272  if (dms.deg == 360) {
273  dms.deg = 0;
274  }
275  }
276  }
277 
278  std::string fmt("%02d:%02d:%05.2f");
279  std::string s = (boost::format(fmt) % dms.deg % dms.min % dms.sec).str();
280  if (dms.sign < 0) {
281  s = "-" + s;
282  }
283  return s;
284 
285 }
286 
287 
294  return angleToXmsString(a, afwGeom::degrees);
295 }
296 
301  return angleToXmsString(a, afwGeom::hours);
302 }
303 
307 static afwGeom::Angle xmsStringToAngle(
308  std::string const dms,
309  afwGeom::AngleUnit unit
310  ) {
311  if (dms.find(":") == std::string::npos) {
312  throw LSST_EXCEPT(ex::InvalidParameterError,
313  (boost::format("String is not in xx:mm:ss format: %s") % dms).str());
314  }
315  std::vector<std::string> elements;
316  boost::split(elements, dms, boost::is_any_of(":"));
317  if (elements.size() != 3) {
318  throw LSST_EXCEPT(ex::InvalidParameterError,
319  (boost::format("Could not parse string as xx:mm:ss format: %s") % dms).str());
320  }
321  int const deg = abs(atoi(elements[0].c_str()));
322  int const min = atoi(elements[1].c_str());
323  double const sec = atof(elements[2].c_str());
324 
325  afwGeom::Angle ang = (deg + min/60.0 + sec/3600.0) * unit;
326  if ( (elements[0].c_str())[0] == '-' ) {
327  ang *= -1.0;
328  }
329  return ang;
330 }
331 
336  std::string const hms
337  ){
338  return xmsStringToAngle(hms, afwGeom::hours);
339 }
340 
345  std::string const dms
346  ) {
347  return xmsStringToAngle(dms, afwGeom::degrees);
348 }
349 
355  double const epoch
356  ) {
357  double const T = (epoch - 2000.0)/100.0;
358  return (23.0 + 26.0/60.0 + (21.448 - 46.82*T - 0.0006*T*T - 0.0018*T*T*T)/3600.0) * afwGeom::degrees;
359 }
360 
361 
362 /* ============================================================
363  *
364  * class Coord
365  *
366  * ============================================================*/
367 
368 
374  lsst::afw::geom::Point2D const &p2d,
375  afwGeom::AngleUnit unit,
376  double const epoch
377  ) :
378  _longitude(NaN), _latitude(NaN), _epoch(epoch) {
379  _longitude = afwGeom::Angle(p2d.getX(), unit);
380  _longitude.wrap();
381  _latitude = afwGeom::Angle(p2d.getY(), unit);
382  _verifyValues();
383 }
384 
389  afwGeom::Point3D const &p3d,
390  double const epoch,
391  bool normalize,
392  afwGeom::Angle const defaultLongitude
393  ) :
394  _longitude(0. * afwGeom::radians),
395  _latitude(0. * afwGeom::radians),
396  _epoch(epoch) {
397 
398  std::pair<afwGeom::Angle, afwGeom::Angle> lonLat = pointToLonLat(p3d, defaultLongitude, normalize);
399  _longitude = lonLat.first;
400  _latitude = lonLat.second;
401  _epoch = epoch;
402 }
403 
409  afwGeom::Angle const ra,
410  afwGeom::Angle const dec,
411  double const epoch
412  ) :
413  _longitude(ra), _latitude(dec), _epoch(epoch) {
414  _longitude.wrap();
415  _verifyValues();
416 }
417 
423  std::string const ra,
424  std::string const dec,
425  double const epoch
426  ) :
427  _longitude(hmsStringToAngle(ra)),
428  _latitude(dmsStringToAngle(dec)),
429  _epoch(epoch) {
430  _longitude.wrap();
431  _verifyValues();
432 }
433 
441 afwCoord::Coord::Coord() : _longitude(afwGeom::Angle(NaN)), _latitude(afwGeom::Angle(NaN)), _epoch(NaN) {}
442 
443 
448  if (_latitude.asRadians() < -afwGeom::HALFPI || _latitude.asRadians() > afwGeom::HALFPI) {
449  throw LSST_EXCEPT(ex::InvalidParameterError,
450  (boost::format("Latitude coord must be: -PI/2 <= lat <= PI/2 (%f).") %
451  _latitude).str());
452  }
453 }
454 
461  afwGeom::Angle const longitude,
462  afwGeom::Angle const latitude,
463  double const epoch
464  ) {
465  _longitude = longitude;
466  _latitude = latitude;
467  _longitude.wrap();
468  _epoch = epoch;
469  _verifyValues();
470 }
471 
472 
478  // treat HOURS specially, they must mean hours for RA, degrees for Dec
479  if (unit == afwGeom::hours) {
480  return afwGeom::Point2D(getLongitude().asHours(), getLatitude().asDegrees());
481  } else {
482  return afwGeom::Point2D(getLongitude().asAngularUnits(unit), getLatitude().asAngularUnits(unit));
483  }
484 }
485 
486 
492  double lng = getLongitude();
493  double lat = getLatitude();
494  double const x = std::cos(lng) * std::cos(lat);
495  double const y = std::sin(lng) * std::cos(lat);
496  double const z = std::sin(lat);
497  return afwGeom::Point3D(x, y, z);
498 }
499 
500 
508  Coord const &poleTo,
509  Coord const &poleFrom
510  ) const {
511  double const alphaGP = poleFrom[0];
512  double const deltaGP = poleFrom[1];
513  double const lCP = poleTo[0];
514 
515  double const alpha = getLongitude();
516  double const delta = getLatitude();
517 
518  afwGeom::Angle const l = (lCP - std::atan2(std::sin(alpha - alphaGP),
519  std::tan(delta)*std::cos(deltaGP) - std::cos(alpha - alphaGP)*std::sin(deltaGP))) * afwGeom::radians;
520  afwGeom::Angle const b = std::asin( (std::sin(deltaGP)*std::sin(delta) + std::cos(deltaGP)*std::cos(delta)*std::cos(alpha - alphaGP))) * afwGeom::radians;
521  return Coord(l, b);
522 }
523 
524 
530  Coord const &axis,
531  afwGeom::Angle const theta
532  ) {
533 
534  double const c = std::cos(theta);
535  double const mc = 1.0 - c;
536  double const s = std::sin(theta);
537 
538  // convert to cartesian
539  afwGeom::Point3D const x = getVector();
540  afwGeom::Point3D const u = axis.getVector();
541  double const ux = u[0];
542  double const uy = u[1];
543  double const uz = u[2];
544 
545  // rotate
546  afwGeom::Point3D xprime;
547  xprime[0] = (ux*ux + (1.0 - ux*ux)*c)*x[0] + (ux*uy*mc - uz*s)*x[1] + (ux*uz*mc + uy*s)*x[2];
548  xprime[1] = (uy*uy + (1.0 - uy*uy)*c)*x[1] + (uy*uz*mc - ux*s)*x[2] + (ux*uy*mc + uz*s)*x[0];
549  xprime[2] = (uz*uz + (1.0 - uz*uz)*c)*x[2] + (uz*ux*mc - uy*s)*x[0] + (uy*uz*mc + ux*s)*x[1];
550 
551  // in-situ
552  std::pair<afwGeom::Angle, afwGeom::Angle> lonLat = pointToLonLat(xprime);
553  _longitude = lonLat.first;
554  _latitude = lonLat.second;
555 }
556 
557 
567  afwGeom::Angle const phi,
568  afwGeom::Angle const arcLen
569  ) {
570 
571  // let v = vector in the direction arcLen points (tangent to surface of sphere)
572  // thus: |v| = arcLen
573  // angle phi = orientation of v in a tangent plane, measured wrt to a parallel of declination
574 
575  // To do the rotation, use rotate() method.
576  // - must provide an axis of rotation: take the cross product r x v to get that axis (pole)
577 
578  // get the vector r
579  Eigen::Vector3d r = getVector().asEigen();
580 
581  // Get the vector v:
582  // let u = unit vector lying on a parallel of declination
583  // let w = unit vector along line of longitude = r x u
584  // the vector v must satisfy the following:
585  // |v| = arcLen
586  // r . v = 0
587  // u . v = |v| cos(phi) = arcLen*cos(phi)
588  // w . v = |v| sin(phi) = arcLen*sin(phi)
589 
590  // v is a linear combination of u and w
591  // v = arcLen*cos(phi)*u + arcLen*sin(phi)*w
592 
593  // Thus, we must:
594  // - create u vector
595  // - solve w vector (r cross u)
596  // - compute v
597  Eigen::Vector3d u;
598  u << -std::sin(getLongitude()), std::cos(getLongitude()), 0.0;
599  Eigen::Vector3d w = r.cross(u);
600  Eigen::Vector3d v = arcLen * std::cos(phi)*u + arcLen * std::sin(phi)*w;
601 
602  // take r x v to get the axis
603  Eigen::Vector3d axisVector = r.cross(v);
604  axisVector.normalize();
605  Coord axisCoord = Coord(afwGeom::Point3D(axisVector), getEpoch());
606 
607  rotate(axisCoord, arcLen);
608 
609  // now get the position angle at our destination
610  // u2 . v2 = arcLen*cos(phi2)
611  // w2 . v2 = arcLen*sin(phi2)
612  // if we normalize v2:
613  // phi2 = atan2(w2.v2, u2.v2)
614  //
615  // we need to compute u2, and then rotate v (exactly as we rotated r) to get v2
616  Eigen::Vector3d r2 = getVector().asEigen();
617  Eigen::Vector3d u2;
618  u2 << -std::sin(getLongitude()), std::cos(getLongitude()), 0.0;
619  Eigen::Vector3d w2 = r2.cross(u2);
620 
621  // make v a unit vector and rotate v exactly as we rotated r
622  v.normalize();
623  Coord v2Coord = Coord(afwGeom::Point3D(v), getEpoch());
624  v2Coord.rotate(axisCoord, arcLen);
625  Eigen::Vector3d v2 = v2Coord.getVector().asEigen();
626 
627  afwGeom::Angle phi2 = std::atan2(w2.dot(v2), u2.dot(v2)) * afwGeom::radians;
628 
629  return phi2;
630 }
631 
632 
633 
638 
639  switch (system) {
640  case FK5:
641  {
642  Fk5Coord c1 = this->toFk5();
643  return boost::shared_ptr<Fk5Coord>(new Fk5Coord(c1.getLongitude(),
644  c1.getLatitude(),
645  c1.getEpoch()));
646  }
647  break;
648  case ICRS:
649  {
650  IcrsCoord c2 = this->toIcrs();
651  return boost::shared_ptr<IcrsCoord>(new IcrsCoord(c2.getLongitude(),
652  c2.getLatitude()));
653  }
654  break;
655  case GALACTIC:
656  {
657  GalacticCoord c4 = this->toGalactic();
658  return boost::shared_ptr<GalacticCoord>(new GalacticCoord(c4.getLongitude(),
659  c4.getLatitude()));
660  }
661  break;
662  case ECLIPTIC:
663  {
664  EclipticCoord c5 = this->toEcliptic();
665  return boost::shared_ptr<EclipticCoord>(new EclipticCoord(c5.getLongitude(),
666  c5.getLatitude(),
667  c5.getEpoch()));
668  }
669  break;
670  case TOPOCENTRIC:
671  throw LSST_EXCEPT(ex::InvalidParameterError,
672  "Cannot make Topocentric with convert() (must also specify Observatory).\n"
673  "Instantiate TopocentricCoord() directly.");
674  break;
675  default:
676  throw LSST_EXCEPT(ex::InvalidParameterError,
677  "Undefined CoordSystem: only FK5, ICRS, GALACTIC, ECLIPTIC allowed.");
678  break;
679 
680  }
681 
682 }
683 
684 
690  Coord const &c
691  ) const {
692 
693  // work in Fk5, no matter what two derived classes we're given (eg Fk5 and Galactic)
694  // we'll put them in the same system.
695  std::pair<afwCoord::Fk5Coord, afwCoord::Fk5Coord> const& fk5 = commonFk5(*this, c);
696  afwGeom::Angle const alpha1 = fk5.first.getRa();
697  afwGeom::Angle const delta1 = fk5.first.getDec();
698  afwGeom::Angle const alpha2 = fk5.second.getRa();
699  afwGeom::Angle const delta2 = fk5.second.getDec();
700 
701  return haversine(alpha1 - alpha2, delta1 - delta2, std::cos(delta1), std::cos(delta2));
702 }
703 
704 
712 std::pair<afwGeom::Angle, afwGeom::Angle>
714  ) const
715 {
716  // work in Fk5, no matter what two derived classes we're given (eg Fk5 and Galactic)
717  // we'll put them in the same system.
718  std::pair<afwCoord::Fk5Coord, afwCoord::Fk5Coord> const& fk5 = commonFk5(*this, c);
719  afwGeom::Angle const alpha1 = fk5.first.getRa();
720  afwGeom::Angle const delta1 = fk5.first.getDec();
721  afwGeom::Angle const alpha2 = fk5.second.getRa();
722  afwGeom::Angle const delta2 = fk5.second.getDec();
723 
724  afwGeom::Angle const dAlpha = alpha1-alpha2;
725  afwGeom::Angle const dDelta = delta1-delta2;
726 
727  double const cosDelta1 = std::cos(delta1);
728  double const cosDelta2 = std::cos(delta2);
729 
730  afwGeom::Angle separation = haversine(dAlpha, dDelta, cosDelta1, cosDelta2);
731 
732  // Formula from http://www.movable-type.co.uk/scripts/latlong.html
733  double const y = std::sin(dAlpha)*cosDelta2;
734  double const x = cosDelta1*std::sin(delta2) - std::sin(delta1)*cosDelta2*std::cos(dAlpha);
735  afwGeom::Angle bearing = std::atan2(y, x)*afwGeom::radians - 90.0*afwGeom::degrees;
736 
737  return std::make_pair(bearing, separation);
738 }
739 
747 std::pair<afwGeom::Angle, afwGeom::Angle>
749  ) const
750 {
751  // work in Fk5, no matter what two derived classes we're given (eg Fk5 and Galactic)
752  // we'll put them in the same system.
753  std::pair<afwCoord::Fk5Coord, afwCoord::Fk5Coord> const& fk5 = commonFk5(*this, c);
754  afwGeom::Angle const alpha1 = fk5.first.getRa();
755  afwGeom::Angle const delta1 = fk5.first.getDec();
756  afwGeom::Angle const alpha2 = fk5.second.getRa();
757  afwGeom::Angle const delta2 = fk5.second.getDec();
758 
759  // This is a projection of coord2 to the tangent plane at coord1
760  double const sinDelta1 = std::sin(delta1);
761  double const cosDelta1 = std::cos(delta1);
762  double const sinDelta2 = std::sin(delta2);
763  double const cosDelta2 = std::cos(delta2);
764  double const cosAlphaDiff = std::cos(alpha2 - alpha1);
765  double const sinAlphaDiff = std::sin(alpha2 - alpha1);
766 
767  double const div = cosDelta1 * cosAlphaDiff * cosDelta2 + sinDelta1 * sinDelta2;
768  double const xi = cosDelta1 * sinAlphaDiff / div;
769  double const eta = (cosDelta1 * cosAlphaDiff * sinDelta2 - sinDelta1 * cosDelta2) / div;
770 
771  return std::make_pair(xi*afwGeom::radians, eta*afwGeom::radians);
772 }
773 
774 
778 afwCoord::Fk5Coord afwCoord::Coord::toFk5(double const epoch) const {
779  return Fk5Coord(getLongitude(), getLatitude(), getEpoch()).precess(epoch);
780 }
785  return Fk5Coord(getLongitude(), getLatitude(), getEpoch());
786 }
787 
793  return this->toFk5().toIcrs();
794 }
795 
800  return this->toFk5().toGalactic();
801 }
802 
807  return this->toFk5(epoch).toEcliptic();
808 }
813  return this->toFk5().toEcliptic();
814 }
815 
820  Observatory const &obs,
821  lsst::daf::base::DateTime const &obsDate
822  ) const {
823  return this->toFk5().toTopocentric(obs, obsDate);
824 }
825 
826 
827 
828 
829 /* ============================================================
830  *
831  * class Fk5Coord
832  *
833  * ============================================================*/
834 
835 
840  return Fk5Coord(getLongitude(), getLatitude(), getEpoch()).precess(epoch);
841 }
846  return Fk5Coord(getLongitude(), getLatitude(), getEpoch());
847 }
848 
854 
855  // only do the precession to 2000 if we're not already there.
856  if ( fabs(getEpoch() - 2000.0) > epochTolerance ) {
857  afwCoord::Fk5Coord c = precess(2000.0);
858  return IcrsCoord(c.getLongitude(), c.getLatitude());
859  } else {
860  return IcrsCoord(getLongitude(), getLatitude());
861  }
862 }
863 
864 
869 
870  // if we're epoch==2000, we can transform, otherwise we need to precess first
871  Fk5Coord c;
872  if ( fabs(getEpoch() - 2000.0) > epochTolerance ) {
873  c = precess(2000.0);
874  } else {
875  c = *this;
876  }
877 
878  Coord ct = c.transform(Fk5PoleInGalactic(), GalacticPoleInFk5());
879  return GalacticCoord(ct.getLongitude(), ct.getLatitude());
880 
881 }
882 
887  afwGeom::Angle const eclPoleIncl = eclipticPoleInclination(epoch);
888  Coord const eclPoleInEquatorial(270.0 * afwGeom::degrees, (90.0 * afwGeom::degrees) - eclPoleIncl, epoch);
889  Coord const equPoleInEcliptic(90.0 * afwGeom::degrees, (90.0 * afwGeom::degrees) - eclPoleIncl, epoch);
890  Coord c = transform(equPoleInEcliptic, eclPoleInEquatorial);
891  return EclipticCoord(c.getLongitude(), c.getLatitude(), epoch);
892 }
897  return this->toEcliptic(getEpoch());
898 }
899 
904  Observatory const &obs,
905  lsst::daf::base::DateTime const &obsDate
906  ) const {
907 
908  // make sure we precess to the epoch
909  Fk5Coord fk5 = precess(obsDate.get(dafBase::DateTime::EPOCH));
910 
911  // greenwich sidereal time
912  afwGeom::Angle theta0 = meanSiderealTimeGreenwich(obsDate.get(dafBase::DateTime::JD));
913  theta0.wrap();
914 
915  // lat/long of the observatory
916  afwGeom::Angle const phi = obs.getLatitude();
917  afwGeom::Angle const L = obs.getLongitude();
918 
919  // ra/dec of the target
920  afwGeom::Angle const alpha = fk5.getRa();
921  afwGeom::Angle const delta = fk5.getDec();
922 
923  afwGeom::Angle const H = theta0 + L - alpha;
924 
925  // compute the altitude, h
926  double const sinh = std::sin(phi)* std::sin(delta) + std::cos(phi) * std::cos(delta) * std::cos(H);
927  afwGeom::Angle const h = std::asin(sinh) * afwGeom::radians;
928 
929  // compute the azimuth, A
930  double const tanAnumerator = std::sin(H);
931  double const tanAdenominator = (std::cos(H) * std::sin(phi) - std::tan(delta) * std::cos(phi));
932 
933  // Equations used here assume azimuth is with respect to South
934  // but we use the North as our origin ... must add 180 deg
935  afwGeom::Angle A = (180.0*afwGeom::degrees) + atan2(tanAnumerator, tanAdenominator)* afwGeom::radians;
936  A.wrap();
937 
938  return TopocentricCoord(A, h, obsDate.get(dafBase::DateTime::EPOCH), obs);
939 }
940 
941 
942 
948  double const epochTo
949  ) const {
950 
951  // return a copy if the epochs are the same
952  if ( fabs(getEpoch() - epochTo) < epochTolerance) {
953  return Fk5Coord(getLongitude(), getLatitude(), getEpoch());
954  }
955 
958  double const jd0 = dateFrom.get(dafBase::DateTime::JD);
959  double const jd = dateTo.get(dafBase::DateTime::JD);
960 
961  double const T = (jd0 - JD2000)/36525.0;
962  double const t = (jd - jd0)/36525.0;
963  double const tt = t*t;
964  double const ttt = tt*t;
965 
966  afwGeom::Angle const xi = ((2306.2181 + 1.39656*T - 0.000139*T*T)*t +
967  (0.30188 - 0.000344*T)*tt + 0.017998*ttt) * afwGeom::arcseconds;
968  afwGeom::Angle const z = ((2306.2181 + 1.39656*T - 0.000139*T*T)*t +
969  (1.09468 + 0.000066*T)*tt + 0.018203*ttt) * afwGeom::arcseconds;
970  afwGeom::Angle const theta = ((2004.3109 - 0.85330*T - 0.000217*T*T)*t -
971  (0.42665 + 0.000217*T)*tt - 0.041833*ttt) * afwGeom::arcseconds;
972 
973  Fk5Coord fk5 = this->toFk5();
974  afwGeom::Angle const alpha0 = fk5.getRa();
975  afwGeom::Angle const delta0 = fk5.getDec();
976 
977  double const a = std::cos(delta0) * std::sin((alpha0 + xi));
978  double const b = std::cos(theta) * std::cos(delta0) * std::cos((alpha0 + xi)) - std::sin(theta) * std::sin(delta0);
979  double const c = std::sin(theta) * std::cos(delta0) * std::cos((alpha0 + xi)) + std::cos(theta) * std::sin(delta0);
980 
981  afwGeom::Angle const alpha = (std::atan2(a,b) + z) * afwGeom::radians;
982  afwGeom::Angle const delta = std::asin(c) * afwGeom::radians;
983 
984  return Fk5Coord(alpha, delta, epochTo);
985 }
986 
987 
988 /* ============================================================
989  *
990  * class IcrsCoord
991  *
992  * ============================================================*/
993 
997 void afwCoord::IcrsCoord::reset(afwGeom::Angle const longitude, afwGeom::Angle const latitude) {
998  Coord::reset(longitude, latitude, 2000.0);
999 }
1000 
1005  return Fk5Coord(getLongitude(), getLatitude(), 2000.0).precess(epoch);
1006 }
1011  return Fk5Coord(getLongitude(), getLatitude(), 2000.0);
1012 }
1013 
1018  return IcrsCoord(getLongitude(), getLatitude());
1019 }
1020 
1021 
1022 
1023 
1024 /* ============================================================
1025  *
1026  * class GalacticCoord
1027  *
1028  * ============================================================*/
1029 
1033 void afwCoord::GalacticCoord::reset(afwGeom::Angle const longitudeDeg, afwGeom::Angle const latitudeDeg) {
1034  Coord::reset(longitudeDeg, latitudeDeg, 2000.0);
1035 }
1036 
1041  // transform to fk5
1042  // galactic coords are ~constant, and the poles used are for epoch=2000, so we get J2000
1043  Coord c = transform(GalacticPoleInFk5(), Fk5PoleInGalactic());
1044  return Fk5Coord(c.getLongitude(), c.getLatitude(), 2000.0).precess(epoch);
1045 }
1050  return this->toFk5(2000.0);
1051 }
1052 
1057  return GalacticCoord(getLongitude(), getLatitude());
1058 }
1059 
1060 
1061 
1062 
1063 /* ============================================================
1064  *
1065  * class EclipticCoord
1066  *
1067  * ============================================================*/
1068 
1073  return EclipticCoord(getLongitude(), getLatitude(), getEpoch()).precess(epoch);
1074 }
1079  return EclipticCoord(getLongitude(), getLatitude(), getEpoch());
1080 }
1081 
1082 
1087  afwGeom::Angle const eclPoleIncl = eclipticPoleInclination(epoch);
1088  afwGeom::Angle ninety = 90. * afwGeom::degrees;
1089  Coord const eclipticPoleInFk5(270.0 * afwGeom::degrees, ninety - eclPoleIncl, epoch);
1090  Coord const fk5PoleInEcliptic(ninety, ninety - eclPoleIncl, epoch);
1091  Coord c = transform(eclipticPoleInFk5, fk5PoleInEcliptic);
1092  return Fk5Coord(c.getLongitude(), c.getLatitude(), epoch);
1093 }
1098  return this->toFk5(getEpoch());
1099 }
1100 
1101 
1108  double const epochTo
1109  ) const {
1110  return this->toFk5().precess(epochTo).toEcliptic();
1111 }
1112 
1113 
1114 
1115 /* ============================================================
1116  *
1117  * class TopocentricCoord
1118  *
1119  * ============================================================*/
1120 
1125 
1126  afwGeom::Angle const phi = _obs.getLatitude();
1127  afwGeom::Angle const L = _obs.getLongitude();
1128 
1129  // Equations used here assume azimuth is with respect to South
1130  // but we use the North as our origin.
1131  afwGeom::Angle A = getAzimuth() + 180.0*afwGeom::degrees;
1132  A.wrap();
1133  afwGeom::Angle const h = getAltitude();
1134 
1135 
1136  double const jd = dafBase::DateTime(epoch,
1139  afwGeom::Angle theta0 = meanSiderealTimeGreenwich(jd);
1140  theta0.wrap();
1141 
1142  double const tanHnum = std::sin(A);
1143  double const tanHdenom = std::cos(A)*std::sin(phi) + std::tan(h)*std::cos(phi);
1144  afwGeom::Angle H = std::atan2(tanHnum, tanHdenom) * afwGeom::radians;
1145 
1146  afwGeom::Angle const alpha = theta0 + L - H;
1147  double const sinDelta = std::sin(phi)*std::sin(h) - std::cos(phi)*std::cos(h)*std::cos(A);
1148  afwGeom::Angle const delta = (std::asin(sinDelta)) * afwGeom::radians;
1149 
1150  return Fk5Coord(alpha, delta, epoch);
1151 }
1156  return this->toFk5(getEpoch());
1157 }
1158 
1159 
1164  lsst::afw::coord::Observatory const &obs,
1165  lsst::daf::base::DateTime const &date
1166  ) const
1167 {
1168  if (obs != _obs) {
1169  throw LSST_EXCEPT(ex::InvalidParameterError,
1170  (boost::format("Expected observatory %s, saw %s") % _obs % obs).str());
1171  }
1172  if (fabs(date.get() - getEpoch()) > std::numeric_limits<double>::epsilon()) {
1173  throw LSST_EXCEPT(ex::InvalidParameterError,
1174  (boost::format("Expected date %g, saw %g") % getEpoch() % date.get()).str());
1175  }
1176 
1177  return TopocentricCoord(getLongitude(), getLatitude(), getEpoch(), _obs);
1178 }
1179 
1186  return TopocentricCoord(getLongitude(), getLatitude(), getEpoch(), _obs);
1187 }
1188 
1189 
1190 
1191 /* ===============================================================================
1192  *
1193  * Factory function definitions:
1194  *
1195  * Each makeCoord() function has overloaded variants with and without epoch specified.
1196  * It was tempting to specify default values, but some coordinate systems have
1197  * no epoch, and this would be confusing to the user. In its current form,
1198  * any epochless coordinate system requested with an epoch specified will throw an exception.
1199  * The epochless factories will always work, but assume default epoch = 2000.
1200  *
1201  * ===============================================================================
1202  */
1203 
1212  CoordSystem const system,
1213  afwGeom::Angle const ra,
1214  afwGeom::Angle const dec,
1215  double const epoch
1216 ) {
1217 
1218  switch (system) {
1219  case FK5:
1220  return boost::shared_ptr<Fk5Coord>(new Fk5Coord(ra, dec, epoch));
1221  break;
1222  case ICRS:
1223  throw LSST_EXCEPT(ex::InvalidParameterError,
1224  "ICRS has no epoch, use overloaded makeCoord with args (system, ra, dec).");
1225  break;
1226  case GALACTIC:
1227  throw LSST_EXCEPT(ex::InvalidParameterError,
1228  "Galactic has no epoch, use overloaded makeCoord with (system, ra, dec).");
1229  break;
1230  case ECLIPTIC:
1231  return boost::shared_ptr<EclipticCoord>(new EclipticCoord(ra, dec, epoch));
1232  break;
1233  case TOPOCENTRIC:
1234  throw LSST_EXCEPT(ex::InvalidParameterError,
1235  "Cannot make Topocentric with makeCoord() (must also specify Observatory).\n"
1236  "Instantiate TopocentricCoord() directly.");
1237  break;
1238  default:
1239  throw LSST_EXCEPT(ex::InvalidParameterError,
1240  "Undefined CoordSystem: only FK5, ICRS, GALACTIC, ECLIPTIC, and TOPOCENTRIC allowed.");
1241  break;
1242 
1243  }
1244 
1245 }
1246 
1247 
1248 
1256  CoordSystem const system,
1257  afwGeom::Angle const ra,
1258  afwGeom::Angle const dec
1259 ) {
1260 
1261  switch (system) {
1262  case FK5:
1263  return boost::shared_ptr<Fk5Coord>(new Fk5Coord(ra, dec, 2000.0));
1264  break;
1265  case ICRS:
1266  return boost::shared_ptr<IcrsCoord>(new IcrsCoord(ra, dec));
1267  break;
1268  case GALACTIC:
1269  return boost::shared_ptr<GalacticCoord>(new GalacticCoord(ra, dec));
1270  break;
1271  case ECLIPTIC:
1272  return boost::shared_ptr<EclipticCoord>(new EclipticCoord(ra, dec, 2000.0));
1273  break;
1274  case TOPOCENTRIC:
1275  throw LSST_EXCEPT(ex::InvalidParameterError,
1276  "Cannot make Topocentric with makeCoord() (must also specify Observatory).\n"
1277  "Instantiate TopocentricCoord() directly.");
1278  break;
1279  default:
1280  throw LSST_EXCEPT(ex::InvalidParameterError,
1281  "Undefined CoordSystem: only FK5, ICRS, GALACTIC, ECLIPTIC, and TOPOCENTRIC allowed.");
1282  break;
1283 
1284  }
1285 
1286 }
1287 
1288 
1289 
1297  CoordSystem const system,
1298  lsst::afw::geom::Point3D const &p3d,
1299  double const epoch,
1300  bool normalize,
1301  afwGeom::Angle const defaultLongitude
1302 ) {
1303  Coord c(p3d, 2000.0, normalize, defaultLongitude);
1304  return makeCoord(system, c.getLongitude(), c.getLatitude(), epoch);
1305 }
1313  CoordSystem const system,
1314  lsst::afw::geom::Point3D const &p3d,
1315  bool normalize,
1316  afwGeom::Angle const defaultLongitude
1317 ) {
1318  Coord c(p3d, 2000.0, normalize, defaultLongitude);
1319  return makeCoord(system, c.getLongitude(), c.getLatitude());
1320 }
1321 
1328  CoordSystem const system,
1329  lsst::afw::geom::Point2D const &p2d,
1330  afwGeom::AngleUnit unit,
1331  double const epoch
1332 ) {
1333  if (unit == afwGeom::hours) {
1334  return makeCoord(system, afwGeom::Angle(p2d.getX(), afwGeom::hours), afwGeom::Angle(p2d.getY(), afwGeom::degrees), epoch);
1335  } else {
1336  return makeCoord(system, afwGeom::Angle(p2d.getX(), unit), afwGeom::Angle(p2d.getY(), unit), epoch);
1337  }
1338 }
1339 
1347  CoordSystem const system,
1348  lsst::afw::geom::Point2D const &p2d,
1349  afwGeom::AngleUnit unit
1350 ) {
1351  if (unit == afwGeom::hours) {
1352  return makeCoord(system, afwGeom::Angle(p2d.getX(), afwGeom::hours), afwGeom::Angle(p2d.getY(), afwGeom::degrees));
1353  } else {
1354  return makeCoord(system, afwGeom::Angle(p2d.getX(), unit), afwGeom::Angle(p2d.getY(), unit));
1355  }
1356 }
1357 
1364  CoordSystem const system,
1365  std::string const ra,
1366  std::string const dec,
1367  double const epoch
1368  ) {
1369  return makeCoord(system, dmsStringToAngle(ra), dmsStringToAngle(dec), epoch);
1370 }
1377  CoordSystem const system,
1378  std::string const ra,
1379  std::string const dec
1380  ) {
1381  return makeCoord(system, dmsStringToAngle(ra), dmsStringToAngle(dec));
1382 }
1383 
1384 
1385 
1390  CoordSystem const system
1391  ) {
1392  switch (system) {
1393  case FK5:
1394  return boost::shared_ptr<Fk5Coord>(new Fk5Coord());
1395  break;
1396  case ICRS:
1397  return boost::shared_ptr<IcrsCoord>(new IcrsCoord());
1398  break;
1399  case GALACTIC:
1400  return boost::shared_ptr<GalacticCoord>(new GalacticCoord());
1401  break;
1402  case ECLIPTIC:
1403  return boost::shared_ptr<EclipticCoord>(new EclipticCoord());
1404  break;
1405  case TOPOCENTRIC:
1406  throw LSST_EXCEPT(ex::InvalidParameterError,
1407  "Cannot make Topocentric with makeCoord() (must also specify Observatory).\n"
1408  "Instantiate TopocentricCoord() directly.");
1409  break;
1410  default:
1411  throw LSST_EXCEPT(ex::InvalidParameterError,
1412  "Undefined CoordSystem: only FK5, ICRS, GALACTIC, ECLIPTIC, allowed.");
1413  break;
1414 
1415  }
1416 }
1417 
1418 std::ostream & afwCoord::operator<<(std::ostream & os, afwCoord::Coord const & coord) {
1419  auto const className = coord.getClassName();
1420  os << (boost::format("%s(%.7f, %.7f")
1421  % className
1422  % coord[0].asDegrees()
1423  % coord[1].asDegrees()).str();
1424  if (className == "TopocentricCoord") {
1425  os << (boost::format(", %.12f, (%s)")
1426  % coord.getEpoch()
1427  % dynamic_cast<afwCoord::TopocentricCoord const &>(coord).getObservatory()).str();
1428  } else if (className != "IcrsCoord" && className != "GalacticCoord") {
1429  os << (boost::format(", %.2f") % coord.getEpoch()).str();
1430  }
1431  os << ")";
1432  return os;
1433 }
int y
virtual Fk5Coord toFk5() const
Convert ourself from Ecliptic to Fk5 (use current epoch)
Definition: Coord.cc:1097
virtual IcrsCoord toIcrs() const
Icrs converter for IcrsCoord. (ie. a no-op)
Definition: Coord.cc:1017
boost::shared_ptr< Coord > Ptr
Definition: Coord.h:72
virtual EclipticCoord toEcliptic() const
Convert ourself to Ecliptic: lambda, beta (use existing epoch)
Definition: Coord.cc:812
Class for handling dates/times, including MJD, UTC, and TAI.
Definition: DateTime.h:58
lsst::afw::coord::Coord Coord
Definition: misc.h:36
lsst::afw::geom::Angle _longitude
Definition: Coord.h:146
virtual GalacticCoord toGalactic() const
Convert ourself to Galactic: l, b.
Definition: Coord.cc:799
virtual void reset(lsst::afw::geom::Angle const longitude, lsst::afw::geom::Angle const latitude)
special reset() overload to make sure no epoch can be set
Definition: Coord.cc:997
A class to handle Galactic coordinates (inherits from Coord)
Definition: Coord.h:237
lsst::afw::geom::Angle hmsStringToAngle(std::string const hms)
Convert a hh:mm:ss string to Angle.
Definition: Coord.cc:335
double asAngularUnits(AngleUnit const &units) const
Definition: Angle.h:118
AngleUnit const radians
constant with units of radians
Definition: Angle.h:91
virtual Fk5Coord toFk5() const
Convert ourself to Fk5 (ie. a no-op): RA, Dec (keep current epoch)
Definition: Coord.cc:845
Store information about an observatory ... lat/long, elevation.
Definition: Observatory.h:48
EclipticCoord precess(double const epochTo) const
precess to new epoch
Definition: Coord.cc:1107
virtual Fk5Coord toFk5() const
Convert outself from Topocentric to Fk5 (use current epoch)
Definition: Coord.cc:1155
AngleUnit const hours
Definition: Angle.h:93
virtual TopocentricCoord toTopocentric(Observatory const &obs, lsst::daf::base::DateTime const &obsDate) const
Convert ourself to Altitude/Azimuth: alt, az.
Definition: Coord.cc:903
std::pair< lsst::afw::geom::Angle, lsst::afw::geom::Angle > getTangentPlaneOffset(Coord const &c) const
Get the offset on the tangent plane.
Definition: Coord.cc:748
Fk5Coord precess(double const epochTo) const
Precess ourselves from whence we are to a new epoch.
Definition: Coord.cc:947
virtual TopocentricCoord toTopocentric(Observatory const &obs, lsst::daf::base::DateTime const &date) const
Convert ourself from Topocentric to Topocentric ... a no-op.
Definition: Coord.cc:1163
virtual TopocentricCoord toTopocentric() const
Convert ourself from Topocentric to Topocentric with no observatory or date arguments.
Definition: Coord.cc:1185
virtual EclipticCoord toEcliptic(double const epoch) const
Convert ourself from Ecliptic to Ecliptic ... a no-op (but precess to new epoch)
Definition: Coord.cc:1072
double get(DateSystem system=MJD, Timescale scale=TAI) const
Definition: DateTime.cc:419
Coord()
Default constructor for the Coord base class.
Definition: Coord.cc:441
std::pair< lsst::afw::geom::Angle, lsst::afw::geom::Angle > getOffsetFrom(Coord const &c) const
Compute the offset from a coordinate.
Definition: Coord.cc:713
lsst::afw::coord::IcrsCoord IcrsCoord
Definition: misc.h:38
virtual IcrsCoord toIcrs() const
Convert ourself to ICRS: RA, Dec (basically J2000)
Definition: Coord.cc:853
void rotate(Coord const &axis, lsst::afw::geom::Angle const theta)
Rotate our current coords about a pole.
Definition: Coord.cc:529
A class used to convert scalar POD types such as double to Angle.
Definition: Angle.h:71
void _verifyValues() const
Make sure the values we&#39;ve got are in the range 0 &lt;= x &lt; 2PI.
Definition: Coord.cc:447
virtual void reset(lsst::afw::geom::Angle const longitude, lsst::afw::geom::Angle const latitude)
Definition: Coord.h:86
A coordinate class intended to represent absolute positions.
Definition: PSF.h:39
Matrix alpha
lsst::afw::geom::Angle _latitude
Definition: Coord.h:147
virtual Fk5Coord toFk5(double const epoch) const
Convert ourself to Fk5: RA, Dec (precess to new epoch)
Definition: Coord.cc:778
AngleUnit const degrees
Definition: Angle.h:92
virtual Fk5Coord toFk5() const
Convert ourself from galactic to Fk5 (no epoch specified)
Definition: Coord.cc:1049
A class to handle topocentric (AltAz) coordinates (inherits from Coord)
Definition: Coord.h:329
virtual Fk5Coord toFk5() const
Convert ourself to Fk5: RA, Dec (use current epoch)
Definition: Coord.cc:784
double w
Definition: CoaddPsf.cc:57
virtual EclipticCoord toEcliptic() const
Convert ourself to Ecliptic: lambda, beta (use current epoch)
Definition: Coord.cc:896
EigenVector const & asEigen() const
Return a fixed-size Eigen representation of the coordinate object.
lsst::afw::geom::Angle getDec() const
Definition: Coord.h:215
Coord::Ptr makeCoord(CoordSystem const system, lsst::afw::geom::Angle const ra, lsst::afw::geom::Angle const dec, double const epoch)
Factory function to create a Coord of arbitrary type with decimal RA,Dec.
Definition: Coord.cc:1211
double getEpoch() const
Definition: Coord.h:93
std::string angleToDmsString(lsst::afw::geom::Angle const deg)
a Function to convert a coordinate in decimal degrees to a string with form dd:mm:ss ...
Definition: Coord.cc:293
lsst::afw::geom::Angle offset(lsst::afw::geom::Angle const phi, lsst::afw::geom::Angle const arcLen)
offset our current coords along a great circle defined by an angle wrt a declination parallel ...
Definition: Coord.cc:566
std::string angleToHmsString(lsst::afw::geom::Angle const deg)
a function to convert decimal degrees to a string with form hh:mm:ss.s
Definition: Coord.cc:300
Interface for DateTime class.
virtual void reset(lsst::afw::geom::Angle const longitude, lsst::afw::geom::Angle const latitude)
special reset() overload to make sure no epoch can be set
Definition: Coord.cc:1033
Point< double, 3 > Point3D
Definition: Point.h:287
int x
A class to handle Fk5 coordinates (inherits from Coord)
Definition: Coord.h:191
virtual GalacticCoord toGalactic() const
Convert ourself to Galactic: l, b.
Definition: Coord.cc:868
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A class to handle Ecliptic coordinates (inherits from Coord)
Definition: Coord.h:278
virtual std::string getClassName() const
Definition: Coord.h:101
tuple m
Definition: lsstimport.py:48
Extent< int, N > floor(Extent< double, N > const &input)
lsst::afw::geom::Angle Angle
Definition: misc.h:39
lsst::afw::geom::Angle getRa() const
Definition: Coord.h:214
lsst::afw::geom::Angle getLongitude() const
The main access method for the longitudinal coordinate.
Definition: Coord.h:113
double const HALFPI
Definition: Angle.h:20
virtual IcrsCoord toIcrs() const
Convert ourself to ICRS: RA, Dec (basically J2000)
Definition: Coord.cc:792
AngleUnit const arcseconds
Definition: Angle.h:95
lsst::afw::geom::Angle getLongitude() const
The main access method for the longitudinal coordinate.
Definition: Observatory.cc:88
afw::table::Key< double > b
Point< double, 2 > Point2D
Definition: Point.h:286
CoordSystem makeCoordEnum(std::string const system)
A utility function to get the enum value of a coordinate system from a string name.
Definition: Coord.cc:117
virtual Fk5Coord toFk5() const
Fk5 converter for IcrsCoord. (no epoch specified)
Definition: Coord.cc:1010
std::ostream & operator<<(std::ostream &stream, Exception const &e)
Definition: Exception.cc:124
virtual GalacticCoord toGalactic() const
Convert ourself from Galactic to Galactic ... a no-op.
Definition: Coord.cc:1056
Coord::Ptr convert(CoordSystem system) const
Convert to a specified Coord type.
Definition: Coord.cc:637
lsst::afw::geom::Point3D getVector() const
Return our contents in a position vector.
Definition: Coord.cc:491
lsst::afw::geom::Angle getLatitude() const
The main access method for the latitudinal coordinate.
Definition: Coord.h:120
lsst::afw::geom::Angle getLatitude() const
The main access method for the longitudinal coordinate.
Definition: Observatory.cc:99
A class to handle Icrs coordinates (inherits from Coord)
Definition: Coord.h:157
lsst::afw::geom::Angle dmsStringToAngle(std::string const dms)
Convert a dd:mm:ss string to Angle.
Definition: Coord.cc:344
virtual EclipticCoord toEcliptic() const
Convert ourself from Ecliptic to Ecliptic ... a no-op (use the current epoch)
Definition: Coord.cc:1078
lsst::afw::geom::Point2D getPosition(lsst::afw::geom::AngleUnit unit=lsst::afw::geom::degrees) const
Return our contents in a Point2D object.
Definition: Coord.cc:477
lsst::afw::geom::Angle eclipticPoleInclination(double const epoch)
get the inclination of the ecliptic pole (obliquity) at epoch
Definition: Coord.cc:354
Include files required for standard LSST Exception handling.
Coord transform(Coord const &poleFrom, Coord const &poleTo) const
Tranform our current coords to another spherical polar system.
Definition: Coord.cc:507
lsst::afw::geom::Angle angularSeparation(Coord const &c) const
compute the angular separation between two Coords
Definition: Coord.cc:689
virtual TopocentricCoord toTopocentric(Observatory const &obs, lsst::daf::base::DateTime const &obsDate) const
Convert ourself to Altitude/Azimuth: alt, az.
Definition: Coord.cc:819
Functions to handle coordinates.