[688] | 1 | #ifndef __TRIPLE_H__ |
---|
| 2 | #define __TRIPLE_H__ |
---|
| 3 | |
---|
| 4 | #include <cmath> |
---|
| 5 | #include <iostream> |
---|
| 6 | |
---|
| 7 | namespace sphereRemap { |
---|
| 8 | |
---|
| 9 | #define EPS 1e-15 |
---|
| 10 | |
---|
| 11 | /* coordinate on sphere */ |
---|
| 12 | class Coord |
---|
| 13 | { |
---|
| 14 | public: |
---|
| 15 | Coord() : x(0), y(0), z(0) {}; |
---|
| 16 | Coord(double x, double y, double z) : x(x), y(y), z(z) {}; |
---|
| 17 | |
---|
| 18 | bool operator==(const Coord& rhs) const |
---|
| 19 | { |
---|
| 20 | return x == rhs.x && y == rhs.y && z == rhs.z; |
---|
| 21 | } |
---|
| 22 | |
---|
| 23 | void rot(const Coord &u, double theta); |
---|
| 24 | |
---|
| 25 | double x, y, z; |
---|
| 26 | }; |
---|
| 27 | |
---|
| 28 | class Vector |
---|
| 29 | { |
---|
| 30 | public: |
---|
| 31 | Vector(double u = 0, double v = 0, double w = 0) : u(u), v(v), w(w) {}; |
---|
| 32 | |
---|
| 33 | /* vector from origin to point on sphere */ |
---|
| 34 | Vector(const Coord& to) : u(to.x), v(to.y), w(to.z) {}; |
---|
| 35 | |
---|
| 36 | Vector& operator+=(const Vector& rhs) |
---|
| 37 | { |
---|
| 38 | u += rhs.u; |
---|
| 39 | v += rhs.v; |
---|
| 40 | w += rhs.w; |
---|
| 41 | } |
---|
| 42 | |
---|
| 43 | double u, v, w; |
---|
| 44 | }; |
---|
| 45 | |
---|
| 46 | extern const Coord ORIGIN; |
---|
| 47 | |
---|
| 48 | std::ostream& operator<<(std::ostream& os, const Coord& c); |
---|
| 49 | std::ostream& operator<<(std::ostream& os, const Vector& c); |
---|
| 50 | |
---|
| 51 | inline double lengthSquared(Vector& vec) |
---|
| 52 | { |
---|
| 53 | return vec.u*vec.u + vec.v*vec.v + vec.w*vec.w; |
---|
| 54 | } |
---|
| 55 | |
---|
| 56 | inline Coord operator+(const Coord &lhs, const Coord &rhs) |
---|
| 57 | { |
---|
| 58 | return Coord(lhs.x + rhs.x, |
---|
| 59 | lhs.y + rhs.y, |
---|
| 60 | lhs.z + rhs.z); |
---|
| 61 | } |
---|
| 62 | |
---|
| 63 | inline Coord operator-(const Coord &lhs, const Coord &rhs) |
---|
| 64 | { |
---|
| 65 | return Coord(lhs.x - rhs.x, |
---|
| 66 | lhs.y - rhs.y, |
---|
| 67 | lhs.z - rhs.z); |
---|
| 68 | } |
---|
| 69 | |
---|
| 70 | /** vector scalar multiplication */ |
---|
| 71 | inline Coord operator*(const Coord &lhs, double rhs) |
---|
| 72 | { |
---|
| 73 | return Coord(lhs.x * rhs, |
---|
| 74 | lhs.y * rhs, |
---|
| 75 | lhs.z * rhs); |
---|
| 76 | } |
---|
| 77 | |
---|
| 78 | inline double squaredist(const Coord &x, const Coord &y) |
---|
| 79 | { |
---|
| 80 | return (y.x-x.x)*(y.x-x.x) |
---|
| 81 | + (y.y-x.y)*(y.y-x.y) |
---|
| 82 | + (y.z-x.z)*(y.z-x.z); |
---|
| 83 | } |
---|
| 84 | |
---|
| 85 | inline double eucldist(const Coord &x, const Coord &y) |
---|
| 86 | { |
---|
| 87 | return sqrt(squaredist(x, y)); |
---|
| 88 | } |
---|
| 89 | |
---|
| 90 | |
---|
| 91 | double norm(const Coord &x); |
---|
| 92 | |
---|
| 93 | Coord proj(const Coord &x); |
---|
| 94 | |
---|
| 95 | double scalarprod(const Coord &a, const Coord &b); |
---|
| 96 | |
---|
| 97 | Coord crossprod(const Coord &a, const Coord &b); |
---|
| 98 | |
---|
| 99 | double arcdist(const Coord &x, const Coord &y); |
---|
| 100 | |
---|
| 101 | double ds(const Coord &x, const Coord &y); |
---|
| 102 | |
---|
| 103 | double dp(const Coord &x, const Coord &y, const Coord &pole); |
---|
| 104 | |
---|
| 105 | //double angle(const Vector&, const Vector&) |
---|
| 106 | |
---|
| 107 | void lonlat(const Coord &a, double &lon, double &lat); |
---|
| 108 | |
---|
| 109 | Coord xyz(double lon, double lat); |
---|
| 110 | |
---|
| 111 | /** Computes the midpoint on spherical arc between two points on the sphere */ |
---|
| 112 | Coord midpoint(const Coord &a, const Coord &b); |
---|
| 113 | |
---|
| 114 | /** Computes the midpoint on *small circle* between two points on the sphere */ |
---|
| 115 | Coord midpointSC(const Coord &a, const Coord &b); |
---|
| 116 | |
---|
| 117 | void rot(const Coord &u, const double th, Coord &x); |
---|
| 118 | |
---|
| 119 | void rotsg(const Coord &a, Coord &b, const Coord &c, |
---|
| 120 | const Coord &d, Coord &x); |
---|
| 121 | |
---|
| 122 | double angle(const Coord &a, const Coord &b, const Coord &pole); |
---|
[1016] | 123 | |
---|
| 124 | // return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b |
---|
| 125 | double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ; |
---|
| 126 | |
---|
[688] | 127 | void print_Coord(Coord &p); |
---|
| 128 | |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | #endif |
---|