source: subversion/applications/lib/ccoord/OSRef.cpp @ 18255

Last change on this file since 18255 was 2197, checked in by nick, 13 years ago

Support for generating grid square letters added

File size: 4.3 KB
Line 
1#include "LatLng.h"
2#include "OSRef.h"
3#include "util.h"
4#include "RefEll.h"
5#include <cmath>
6
7#include <iostream>
8#include <sstream>
9#include <iomanip>
10
11using std::setfill;
12using std::setw;
13
14// Based on Jonathan Stott's JCoord.
15// Licenced under the GNU GPL.
16
17LatLng OSRef::toLatLng()
18
19  {
20    double OSGB_F0 = 0.9996012717;
21    double N0 = -100000.0;
22    double E0 = 400000.0;
23    double phi0 = toRadians(49.0);
24    double lambda0 = toRadians(-2.0);
25        RefEll AIRY_1830 (6377563.396, 6356256.909);
26    double a = AIRY_1830.getMaj();
27    double b = AIRY_1830.getMin();
28    double eSquared = AIRY_1830.getEcc();
29    double phi = 0.0;
30    double lambda = 0.0;
31    double E = this->easting;
32    double N = this->northing;
33    double n = (a - b) / (a + b);
34    double M = 0.0;
35    double phiPrime = ((N - N0) / (a * OSGB_F0)) + phi0;
36    do {
37      M =
38          (b * OSGB_F0)
39              * (((1 + n + ((5.0 / 4.0) * n * n) + ((5.0 / 4.0) * n * n * n)) * (phiPrime - phi0))
40                  - (((3 * n) + (3 * n * n) + ((21.0 / 8.0) * n * n * n))
41                      * sin(phiPrime - phi0) * cos(phiPrime + phi0))
42                  + ((((15.0 / 8.0) * n * n) + ((15.0 / 8.0) * n * n * n))
43                      * sin(2.0 * (phiPrime - phi0)) * 
44                      cos(2.0 * (phiPrime + phi0))) - (((35.0 / 24.0) * n * n
45                                                                * n)
46                  * sin(3.0 * (phiPrime - phi0)) * 
47                  cos(3.0 * (phiPrime + phi0))));
48      phiPrime += (N - N0 - M) / (a * OSGB_F0);
49    } while ((N - N0 - M) >= 0.001);
50    double v =
51        a * OSGB_F0
52            * pow(1.0 - eSquared * sinSquared(phiPrime), -0.5);
53    double rho =
54        a * OSGB_F0 * (1.0 - eSquared)
55            * pow(1.0 - eSquared * sinSquared(phiPrime), -1.5);
56    double etaSquared = (v / rho) - 1.0;
57    double VII = tan(phiPrime) / (2 * rho * v);
58    double VIII =
59        (tan(phiPrime) / (24.0 * rho * pow(v, 3.0)))
60            * (5.0 + (3.0 * tanSquared(phiPrime)) + etaSquared - (9.0 * 
61                tanSquared(phiPrime) * etaSquared));
62    double IX =
63        (tan(phiPrime) / (720.0 * rho * pow(v, 5.0)))
64            * (61.0 + (90.0 * tanSquared(phiPrime)) + (45.0 * 
65                tanSquared(phiPrime) * tanSquared(phiPrime)));
66    double X = sec(phiPrime) / v;
67    double XI =
68        (sec(phiPrime) / (6.0 * v * v * v))
69            * ((v / rho) + (2 * tanSquared(phiPrime)));
70    double XII =
71        (sec(phiPrime) / (120.0 * pow(v, 5.0)))
72            * (5.0 + (28.0 * tanSquared(phiPrime)) + (24.0 * 
73                tanSquared(phiPrime) * tanSquared(phiPrime)));
74    double XIIA =
75        (sec(phiPrime) / (5040.0 * pow(v, 7.0)))
76            * (61.0
77                + (662.0 * tanSquared(phiPrime))
78                + (1320.0 * tanSquared(phiPrime) * 
79                    tanSquared(phiPrime)) + (720.0 * tanSquared(phiPrime)
80                * tanSquared(phiPrime) * tanSquared(phiPrime)));
81    phi =
82        phiPrime - (VII * pow(E - E0, 2.0))
83            + (VIII * pow(E - E0, 4.0)) - (IX * pow(E - E0, 6.0));
84    lambda =
85        lambda0 + (X * (E - E0)) - (XI * pow(E - E0, 3.0))
86            + (XII * pow(E - E0, 5.0)) - (XIIA * pow(E - E0, 7.0));
87
88    return LatLng(toDegrees(phi), toDegrees(lambda));
89  }
90
91  /**
92   * Return a String representation of this OSGB grid reference using the
93   * six-figure notation in the form XY123456
94   *
95   * @return a String representing this OSGB grid reference in six-figure
96   *         notation
97   * @since 1.0
98   */
99        std::string OSRef::toSixFigureString() 
100        {
101    int hundredkmE = (int) floor(easting / 100000);
102    int hundredkmN = (int) floor(northing / 100000);
103        char firstLetter;
104    if (hundredkmN < 5) {
105      if (hundredkmE < 5) {
106        firstLetter = 'S';
107      } else {
108        firstLetter = 'T';
109      }
110    } else if (hundredkmN < 10) {
111      if (hundredkmE < 5) {
112        firstLetter = 'N';
113      } else {
114        firstLetter = 'O';
115      }
116    } else {
117      firstLetter = 'H';
118    }
119
120    int index = 65 + ((4 - (hundredkmN % 5)) * 5) + (hundredkmE % 5);
121    // int ti = index;
122    if (index >= 73)
123      index++;
124        char secondLetter = (char) index;
125
126    int e = (int) floor((easting - (100000 * hundredkmE)) / 100);
127    int n = (int) floor((northing - (100000 * hundredkmN)) / 100);
128
129        std::ostringstream sstream;
130        sstream << firstLetter << secondLetter << setfill('0') << 
131                        setw(3) << e << setfill('0') << setw(3) <<  n;
132        return sstream.str();
133  }
Note: See TracBrowser for help on using the repository browser.