source: subversion/applications/utils/srtm2shp/tomerc.cpp @ 30195

Last change on this file since 30195 was 6182, checked in by nick, 12 years ago

new and internationally-capable srtm2shp added

File size: 1.3 KB
Line 
1#include <cmath>
2#include "Map.h"
3#include "tomerc.h"
4
5EarthPoint lltomerc(const EarthPoint& ep)
6{
7        EarthPoint ep2;
8        ep2.x = lonToMerc(ep.x);
9        ep2.y = latToMerc(ep.y);
10        return ep2;
11}
12
13EarthPoint merctoll(const EarthPoint& ep)
14{
15        EarthPoint ep2;
16        ep2.x = mercToLon(ep.x);
17        ep2.y = mercToLat(ep.y);
18        return ep2;
19}
20
21
22// www.upcnet.es/~jgc2/docs/mercator.html
23double lonToMerc(double lon)
24{
25    lon *= (M_PI/180);
26        double a = 6378137;
27        return a*lon;
28}
29
30double latToMerc(double lat)
31{
32        lat *= (M_PI/180);
33        double a = 6378137;
34        double b = 6356752.3142;
35        double f=  (a-b)/a;
36        double e = sqrt(2*f-pow(f,2));
37
38        double res = a*log( tan(M_PI/4 + lat/2) *
39                                pow(( (1-e*sin(lat)) / (1 +e *sin(lat))) , e/2) );
40
41        return res;
42}
43
44double mercToLon(double merc_x)
45{
46        return (merc_x/6378137.0)*(180.0/M_PI);
47}
48
49double mercToLat(double merc_y)
50{
51        double a = 6378137.0;
52        double b = 6356752.3142;
53        double t = 1.0 - b/a;
54        double es=2*t -t*t;
55        double e = sqrt(es);
56        double lat = phi2(exp(-merc_y/a),e);
57        return lat*(180.0/M_PI);
58}
59
60double phi2(double ts, double e)
61{
62        double eccnth=0.5*e;
63        double phi = (M_PI/2) - 2.0*atan(ts);
64        double dphi;
65        double con;
66        int i=15;
67        do
68        {
69                con=e*sin(phi);
70                dphi = (M_PI/2) - 2.0*atan(ts*pow((1.0-con)/(1.0+con),eccnth))-phi;
71                phi+=dphi;
72        }
73        while(abs(dphi)>0.0000000001 && --i);
74        return phi;
75}
Note: See TracBrowser for help on using the repository browser.