source: XIOS/dev/dev_olga/extern/remap/src/triple.hpp @ 1158

Last change on this file since 1158 was 1158, checked in by oabramkina, 4 years ago

Two server levels: merging with trunk r1137.
There are bugs.

File size: 2.8 KB
Line 
1#ifndef  __TRIPLE_H__
2#define __TRIPLE_H__
3
4#include <cmath>
5#include <iostream>
6
7namespace sphereRemap {
8
9#define EPS 1e-15
10
11/* coordinate on sphere */
12class Coord
13{
14public:
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
28class Vector
29{
30public:
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
46extern const Coord ORIGIN;
47
48std::ostream& operator<<(std::ostream& os, const Coord& c);
49std::ostream& operator<<(std::ostream& os, const Vector& c);
50
51inline double lengthSquared(Vector& vec)
52{
53        return vec.u*vec.u + vec.v*vec.v + vec.w*vec.w;
54}
55       
56inline 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
63inline 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 */
71inline Coord operator*(const Coord &lhs, double rhs)
72{
73        return Coord(lhs.x * rhs, 
74                     lhs.y * rhs, 
75                     lhs.z * rhs);
76}
77
78inline 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
85inline double eucldist(const Coord &x, const Coord &y)
86{
87        return sqrt(squaredist(x, y));
88}
89
90
91double norm(const Coord &x);
92
93Coord proj(const Coord &x);
94
95double scalarprod(const Coord &a, const Coord &b);
96
97Coord crossprod(const Coord &a, const Coord &b);
98
99double arcdist(const Coord &x, const Coord &y);
100
101double ds(const Coord &x, const Coord &y);
102
103double dp(const Coord &x, const Coord &y, const Coord &pole);
104
105//double angle(const Vector&, const Vector&)
106
107void lonlat(const Coord &a, double &lon, double &lat);
108
109Coord xyz(double lon, double lat);
110
111/** Computes the midpoint on spherical arc between two points on the sphere */
112Coord midpoint(const Coord &a, const Coord &b);
113
114/** Computes the midpoint on *small circle* between two points on the sphere */
115Coord midpointSC(const Coord &a, const Coord &b);
116
117void rot(const Coord &u, const double th, Coord &x);
118
119void rotsg(const Coord &a, Coord &b, const Coord &c,
120           const Coord &d, Coord &x);
121
122double angle(const Coord &a, const Coord &b, const Coord &pole);
123
124// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b
125double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ;
126
127void print_Coord(Coord &p);
128
129}
130
131#endif
Note: See TracBrowser for help on using the repository browser.