1 /*
2 SPDX-FileCopyrightText: 2001-2005 Jason Harris <jharris@30doradus.org>
3 SPDX-FileCopyrightText: 2003-2005 Pablo de Vicente <p.devicente@wanadoo.es>
4
5 SPDX-License-Identifier: GPL-2.0-or-later
6 */
7
8 #include "geolocation.h"
9
10 #include "timezonerule.h"
11
GeoLocation(const dms & lng,const dms & lat,const QString & name,const QString & province,const QString & country,double tz,TimeZoneRule * tzrule,double elevation,bool readOnly,int iEllips)12 GeoLocation::GeoLocation(const dms &lng, const dms &lat, const QString &name, const QString &province, const QString &country,
13 double tz, TimeZoneRule *tzrule, double elevation, bool readOnly, int iEllips) :
14 Longitude(lng), Latitude(lat)
15 {
16 Name = name;
17 Province = province;
18 Country = country;
19 TimeZone = tz;
20 TZrule = tzrule;
21 Elevation = elevation;
22 indexEllipsoid = iEllips;
23 ReadOnly = readOnly;
24 setEllipsoid(indexEllipsoid);
25 geodToCart();
26 }
27
GeoLocation(double x,double y,double z,const QString & name,const QString & province,const QString & country,double TZ,TimeZoneRule * tzrule,double elevation,bool readOnly,int iEllips)28 GeoLocation::GeoLocation(double x, double y, double z, const QString &name, const QString &province,
29 const QString &country, double TZ, TimeZoneRule *tzrule, double elevation, bool readOnly, int iEllips)
30 {
31 PosCartX = x;
32 PosCartY = y;
33 PosCartZ = z;
34 Name = name;
35 Province = province;
36 Country = country;
37 TimeZone = TZ;
38 TZrule = tzrule;
39 Elevation = elevation;
40 indexEllipsoid = iEllips;
41 ReadOnly = readOnly;
42 setEllipsoid(indexEllipsoid);
43 cartToGeod();
44 }
45
fullName() const46 QString GeoLocation::fullName() const
47 {
48 if (country().isEmpty())
49 return translatedName();
50 else if (province().isEmpty())
51 {
52 return QString("%1, %2").arg(translatedName(), translatedCountry());
53 }
54 else
55 {
56 return QString("%1, %2, %3").arg(translatedName(), translatedProvince(), translatedCountry());
57 }
58 }
59
setEllipsoid(int i)60 void GeoLocation::setEllipsoid(int i)
61 {
62 static const double A[] = { 6378140.0, 6378137.0, 6378137.0, 6378137.0, 6378136.0 };
63 static const double F[] = { 0.0033528131779, 0.0033528106812, 0.0033528131779, 0.00335281066474, 0.0033528131779 };
64
65 Q_ASSERT(i >= 0 && (unsigned int)i < sizeof(A) / sizeof(A[0]) && "Index must be in bounds");
66 axis = A[i];
67 flattening = F[i];
68 }
69
changeEllipsoid(int index)70 void GeoLocation::changeEllipsoid(int index)
71 {
72 setEllipsoid(index);
73 cartToGeod();
74 }
75
translatedName() const76 QString GeoLocation::translatedName() const
77 {
78 QString context;
79 if (province().isEmpty())
80 {
81 context = QString("City in %1").arg(country());
82 }
83 else
84 {
85 context = QString("City in %1 %2").arg(province(), country());
86 }
87 return Name.isEmpty() ? QString() : i18nc(context.toUtf8().data(), Name.toUtf8().data());
88 }
89
translatedProvince() const90 QString GeoLocation::translatedProvince() const
91 {
92 return Province.isEmpty() ?
93 QString() :
94 i18nc(QString("Region/state in " + country()).toUtf8().data(), Province.toUtf8().data());
95 }
96
translatedCountry() const97 QString GeoLocation::translatedCountry() const
98 {
99 return Country.isEmpty() ? QString() : i18nc("Country name", Country.toUtf8().data());
100 }
101
cartToGeod()102 void GeoLocation::cartToGeod()
103 {
104 static const double RIT = 2.7778e-6;
105 double e2, rpro, lat1, xn, sqrtP2, latd, sinl;
106
107 e2 = 2 * flattening - flattening * flattening;
108
109 sqrtP2 = sqrt(PosCartX * PosCartX + PosCartY * PosCartY);
110
111 rpro = PosCartZ / sqrtP2;
112 latd = atan2(rpro, (1 - e2));
113 lat1 = 0.;
114
115 while (fabs(latd - lat1) > RIT)
116 {
117 double s1 = sin(latd);
118
119 lat1 = latd;
120 xn = axis / (sqrt(1 - e2 * s1 * s1));
121 latd = atan2(static_cast<long double>(rpro) * (1 + e2 * xn * s1), PosCartZ);
122 }
123
124 sinl = sin(latd);
125 xn = axis / (sqrt(1 - e2 * sinl * sinl));
126
127 Elevation = sqrtP2 / cos(latd) - xn;
128 Longitude.setRadians(atan2(PosCartY, PosCartX));
129 Latitude.setRadians(latd);
130 }
131
geodToCart()132 void GeoLocation::geodToCart()
133 {
134 double e2, xn;
135 double sinLong, cosLong, sinLat, cosLat;
136
137 e2 = 2 * flattening - flattening * flattening;
138
139 Longitude.SinCos(sinLong, cosLong);
140 Latitude.SinCos(sinLat, cosLat);
141
142 xn = axis / (sqrt(1 - e2 * sinLat * sinLat));
143 PosCartX = (xn + Elevation) * cosLat * cosLong;
144 PosCartY = (xn + Elevation) * cosLat * sinLong;
145 PosCartZ = (xn * (1 - e2) + Elevation) * sinLat;
146 }
147
TopocentricVelocity(double vtopo[],const dms & gst)148 void GeoLocation::TopocentricVelocity(double vtopo[], const dms &gst)
149 {
150 double Wearth = 7.29211510e-5; // rads/s
151 dms angularVEarth;
152
153 dms time = GSTtoLST(gst);
154 // angularVEarth.setRadians(time.Hours()*Wearth*3600.);
155 double se, ce;
156 // angularVEarth.SinCos(se,ce);
157 time.SinCos(se, ce);
158
159 double d0 = sqrt(PosCartX * PosCartX + PosCartY * PosCartY);
160 // km/s
161 vtopo[0] = -d0 * Wearth * se / 1000.;
162 vtopo[1] = d0 * Wearth * ce / 1000.;
163 vtopo[2] = 0.;
164 }
165
LMST(double jd)166 double GeoLocation::LMST(double jd)
167 {
168 int divresult;
169 double ut, tu, gmst, theta;
170
171 ut = (jd + 0.5) - floor(jd + 0.5);
172 jd -= ut;
173 tu = (jd - 2451545.) / 36525.;
174
175 gmst = 24110.54841 + tu * (8640184.812866 + tu * (0.093104 - tu * 6.2e-6));
176 divresult = (int)((gmst + 8.6400e4 * 1.00273790934 * ut) / 8.6400e4);
177 gmst = (gmst + 8.6400e4 * 1.00273790934 * ut) - (double)divresult * 8.6400e4;
178 theta = 2. * dms::PI * gmst / (24. * 60. * 60.);
179 divresult = (int)((theta + Longitude.radians()) / (2. * dms::PI));
180 return ((theta + Longitude.radians()) - (double)divresult * 2. * dms::PI);
181 }
182
isReadOnly() const183 bool GeoLocation::isReadOnly() const
184 {
185 return ReadOnly;
186 }
187
setReadOnly(bool value)188 void GeoLocation::setReadOnly(bool value)
189 {
190 ReadOnly = value;
191 }
192
UTtoLT(const KStarsDateTime & ut) const193 KStarsDateTime GeoLocation::UTtoLT(const KStarsDateTime &ut) const
194 {
195 KStarsDateTime lt = ut.addSecs(int(3600. * TZ()));
196 // 2017-08-11 (Jasem): setUtcOffset is deprecated
197 //lt.setUtcOffset(int(3600. * TZ()));
198 lt.setTimeSpec(Qt::LocalTime);
199
200 return lt;
201 }
202
LTtoUT(const KStarsDateTime & lt) const203 KStarsDateTime GeoLocation::LTtoUT(const KStarsDateTime <) const
204 {
205 KStarsDateTime ut = lt.addSecs(int(-3600. * TZ()));
206 ut.setTimeSpec(Qt::UTC);
207 // 2017-08-11 (Jasem): setUtcOffset is deprecated
208 //ut.setUtcOffset(0);
209
210 return ut;
211 }
212
distanceTo(const dms & longitude,const dms & latitude)213 double GeoLocation::distanceTo(const dms &longitude, const dms &latitude)
214 {
215 const double R = 6378.135;
216
217 double diffLongitude = (Longitude - longitude).radians();
218 double diffLatitude = (Latitude - latitude).radians();
219
220 double a = sin(diffLongitude / 2) * sin(diffLongitude / 2) + cos(Longitude.radians()) * cos(longitude.radians()) *
221 sin(diffLatitude / 2) * sin(diffLatitude / 2);
222 double c = 2 * atan2(sqrt(a), sqrt(1 - a));
223
224 double distance = R * c;
225
226 return distance;
227 }
228