- added Calc::distance_in_meter, Calc::distance_in_miles
[activitydiary:activitydiary.git] / src / Calc.cpp
1 /*
2  *    ActivityDiary - The smart way to manage your GPS tracks
3  *
4  *    Copyright (C) 2009
5  *
6  *                - Tom Patzig <tom@activitydiary.org>
7  *                - Thomas Goettlicher <thomas@activitydiary.org>
8  *
9  *    This file is part of ActivityDiary.
10  *
11  *    ActivityDiary is free software: you can redistribute it and/or modify
12  *    it under the terms of the GNU General Public License as published by
13  *    the Free Software Foundation, either version 3 of the License, or
14  *    (at your option) any later version.
15  *
16  *    ActivityDiary is distributed in the hope that it will be useful,
17  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
18  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19  *    GNU General Public License for more details.
20  *
21  *    You should have received a copy of the GNU General Public License
22  *    along with ActivityDiary.  If not, see <http://www.gnu.org/licenses/>.
23 */
24
25 #include "Calc.h"
26 #include <QDebug>
27
28
29 Calc::Calc()
30 {
31 }
32
33 Calc::~Calc()
34 {
35 }
36
37
38 float Calc::projectF(float lat)
39 {
40     //??
41     lat = degToRad(lat);
42     return log(tan(lat) + (1.0 / cos(lat)));
43 }
44
45
46 float Calc::projectMercToLat(float v) 
47 {
48     //??
49     return 180.0 / M_PI * atan(sinh(v));
50 }
51
52
53 QRectF Calc::calcLatLonFromTileXY(QPoint p, long zoom)
54 {
55     //??
56     float unit = 1.0 / pow(2 , zoom);
57
58     double relY1 = p.y() * unit;
59     double relY2 = relY1 + unit;
60     double LimitY = projectF(85.0511);
61     double RangeY = 2.0 * LimitY;
62     relY1 = LimitY - RangeY * relY1;
63     relY2 = LimitY - RangeY * relY2;
64     float Lat1 = projectMercToLat(relY1);
65     float Lat2 = projectMercToLat(relY2);
66     unit = 360.0 / pow(2, zoom);
67     float Lon1 = -180.0 + p.x() * unit;
68     float Lon2 = Lon1 + unit;
69
70     // Lat2, Long1, Lat1, Long1 + $Unit)); # S,W,N,E
71
72     return QRectF(Lat1,Lon1,Lat2-Lat1,Lon2-Lon1);
73 }
74
75
76 QPoint Calc::calcTileXY(float lat, float lon, long zoom)
77 {
78     //??
79
80     float xf = (lon + 180.0) / 360.0 * pow (2.0, zoom);
81     float yf = (1.0 - log(tan(lat * M_PI / 180.0) + 1.0 / cos(lat * M_PI / 180.0)) / M_PI) / 2.0 * pow(2, zoom);
82
83     return QPoint(xf,yf);
84 }
85
86
87
88 double Calc::degToRad(double d)
89 {
90     //OK
91     return d / 180.0 * M_PI;
92 }
93
94 double Calc::radToDeg(double r)
95 {
96     //OK
97     return (r * 180.0) / M_PI;
98 }
99
100
101
102 // helper function - converts a latitude at a certain zoom into a y pixel
103 int Calc::latitudeToYAtZoom(double lat, int zoom ) 
104 {
105     //OK
106     double arc = (double) EARTHCIRCUM / ((1 << zoom) * 256.0);
107     double sinLat = sin(degToRad(lat));
108     double metersY = (double) EARTHRADIUS / 2.0 * log((1.0 + sinLat) / (1.0 - sinLat));
109     return round (( ((double) EARTHCIRCUM/2.0)  - metersY ) / arc) ;
110 }
111
112 // helper function - converts a y pixel at a certain zoom to a latitude
113 double Calc::yAtZoomToLatitude(int y, int zoom ) 
114
115 {
116     //OK
117     double arc = (double) EARTHCIRCUM / ((1 << zoom) * 256.0);
118
119     double arg = exp ( ( (EARTHCIRCUM/2.0) - ( y * arc ) ) / (EARTHRADIUS / 2.0) );
120
121     double sinlat = ( (arg - 1.0 ) / (arg + 1.0) );
122
123     return ( radToDeg ( asin( sinlat )));
124 }
125
126 //helper function - converts a longitude at a certain zoom into a x pixel
127 int Calc::longitudeToXAtZoom(double lon, int zoom )
128 {
129     //OK
130     double arc = (double) EARTHCIRCUM / ((1 << zoom) * 256.0);
131     double metersX = (double) EARTHRADIUS * degToRad(lon);
132     return round ((((double) EARTHCIRCUM/2.0)  + metersX) / arc);
133 }
134
135 //helper function - converts x pixel at a cetrain zoom  to a longitude 
136 double Calc::xAtZoomToLongitude(int x, int zoom )
137 {
138     //OK
139     double arc = (double) EARTHCIRCUM / ((1 << zoom) * 256.0);
140
141     double lon = (( (double) x * arc ) - (EARTHCIRCUM/2) ) / (double) EARTHRADIUS;
142     return radToDeg(lon);
143 }
144
145 float Calc::distance_in_miles(float lat1, float lon1, float lat2, float lon2)
146 {
147 //  float distance = (acos(sin(lat1/180*M_PI)*sin(lat2/180*M_PI) + cos(lat1/180*M_PI)*cos(lat2/180*M_PI)*cos(lon1/180*M_PI-lon2/180*M_PI))*180*60/M_PI);
148
149     float distance = (acos(sin(degToRad(lat1)) * sin(degToRad(lat2)) + 
150                     cos(degToRad(lat1)) * cos(degToRad(lat2)) * cos(degToRad(lon1)-degToRad(lon2))) * 180*60/M_PI);
151
152     // miles in km; return 0 if distance is NaN
153     return ( distance >= 0 ? distance : 0 );
154 }
155
156 float Calc::distance_in_km(float lat1, float lon1, float lat2, float lon2)
157 {
158     return ( distance_in_miles(lat1,lon1,lat2,lon2) * MILES_PER_KM );
159 }
160
161 float Calc::distance_in_meter(float lat1, float lon1, float lat2, float lon2)
162 {
163     return ( distance_in_miles(lat1,lon1,lat2,lon2) * MILES_PER_M );
164 }
165
166 /*
167 float Calc::acos(float f)
168 {
169     return (atan2(sqrt(abso(1 - pow(f,2))), f) );
170 }
171 */
172
173 float Calc::abso(float x)
174 {
175     return x >= 0 ? x : -x;
176 }