# source:XIOS/trunk/extern/remap/src/triple.cpp@688

Last change on this file since 688 was 688, checked in by mhnguyen, 9 years ago

Integrating remap library into XIOS

+) Change name of some files of remap library to be compatible with XIOS
+) Implement function to fill in automatically boundary longitude and latitude

Test
+) On Curie
+) test_remap correct

File size: 3.1 KB
Line
1#include "triple.hpp"
2
3namespace sphereRemap {
4
5extern const Coord ORIGIN(0.0, 0.0, 0.0);
6
7std::ostream& operator<<(std::ostream& os, const Coord& c) {
8    return os << "(" << c.x << ", " << c.y << ", " << c.z << ")";
9}
10
11std::ostream& operator<<(std::ostream& os, const Vector& vec)
12{
13    return os << "(" << vec.u << ", " << vec.v << ", " << vec.w << ")";
14}
15
16double norm(const Coord &x)
17{
18        return sqrt(x.x*x.x + x.y*x.y + x.z*x.z);
19}
20
21Coord proj(const Coord &x)
22{
23        double n=norm(x);
24        return Coord(x.x/n, x.y/n, x.z/n);
25}
26
27double scalarprod(const Coord &a, const Coord &b)
28{
29        return a.x*b.x + a.y*b.y + a.z*b.z;
30}
31
32Coord crossprod(const Coord &a, const Coord &b)
33{
34        return Coord(a.y*b.z - a.z*b.y,
35                     a.z*b.x - a.x*b.z,
36                     a.x*b.y - a.y*b.x);
37}
38
39/* arc distance (along great circle) */
40double arcdist(const Coord &x, const Coord &y)
41{ // et angles aigus non orientes
42        Coord n = crossprod(x, y);
43        double a = asin(norm(n));
44        if (squaredist(x,y) > 2.0) a = M_PI - a;
45        return a;
46}
47
48/* alternate way to compute distance along great circle */
49double ds(const Coord &x, const Coord &y)
50{
51        return 2*asin(0.5*eucldist(x,y));
52}
53
54/* dp is both small circle distances */
55double dp(const Coord &x, const Coord &y, const Coord &pole)
56{
57        return 2*asin(0.5 * eucldist(x,y) / sin(ds(pole,x)));
58}
59
60//* angle between `p` and `q` */
61/*double angle(const Vector &p, const Vector &q)
62{
63        return acos(scalarprod(p, q));
64}*/
65
66void lonlat(const Coord &a, double &lon, double &lat)
67{
68        lon = atan2(a.y, a.x) * 180.0/M_PI;
69        lat = (M_PI_2 - acos(a.z)) * 180.0/M_PI;
70}
71
72Coord xyz(double lon, double lat)
73{
74        double phi   = M_PI/180.0*lon;
75        double theta = M_PI_2 - M_PI/180.0*lat;
76        return Coord(sin(theta)*cos(phi),
77                     sin(theta)*sin(phi),
78                     cos(theta));
79}
80
81/** computes the midpoint on spherical arc between a and b */
82Coord midpoint(const Coord &a, const Coord &b)
83{
84        return proj(a + b);
85}
86
87/** computes the midpoint on *small circle* between a and b */
88Coord midpointSC(const Coord& p, const Coord& q)
89{
90        double phi1 = atan2(p.y, p.x);
91        double phi2 = atan2(q.y, q.x);
92        if (phi1*phi2 < 0)
93                phi1 += (phi1 < phi2) ? 2*M_PI : -2*M_PI;
94        double theta = acos(p.z);
95        double phi   = 0.5*(phi1 + phi2);
96        return Coord(sin(theta)*cos(phi),
97                     sin(theta)*sin(phi),
98                     cos(theta));
99}
100
101/* rotates us by angle theta around u (r is rotatiion matrix) */
102void Coord::rot(const Coord &u, double theta)
103{
104        double x = this->x;
105        double y = this->y;
106        double z = this->z;
107
108        double ux2 = u.x*u.x, uy2 = u.y*u.y, uz2 = u.z*u.z;
109        double k = 1 - cos(theta);
110
111        double r00 = ux2      + (1-ux2)*cos(theta),  r01 = u.x*u.y*k - u.z*sin(theta),  r02 = u.x*u.z*k + u.y*sin(theta);
112        double r10 = u.x*u.y*k + u.z*sin(theta),  r11 = uy2      + (1-uy2)*cos(theta),  r12 = u.y*u.z*k - u.x*sin(theta);
113        double r20 = u.x*u.z*k - u.y*sin(theta),  r21 = u.y*u.z*k + u.x*sin(theta),  r22 = uz2 + (1-uz2)*cos(theta);
114
115        this->x = r00*x + r01*y + r02*z;
116        this->y = r10*x + r11*y + r12*z;
117        this->z = r20*x + r21*y + r22*z;
118}
119
120double angle(const Coord &a, const Coord &b, const Coord &pole)
121{
122        return scalarprod(crossprod(a, b), pole);
123}
124
125}
Note: See TracBrowser for help on using the repository browser.