kstars/kstars/geolocation.cpp

   1 /***************************************************************************
   2                           geolocation.cpp  -  K Desktop Planetarium
   3                              -------------------
   4     begin                : Sun Feb 11 2001
   5     copyright            : (C) 2001-2005 by Jason Harris
   6     email                : jharris@30doradus.org
   7     copyright            : (C) 2003-2005 by Pablo de Vicente
   8     email                : p.devicente@wanadoo.es
   9  ***************************************************************************/
  10
  11 /***************************************************************************
  12  *                                                                         *
  13  *   This program is free software; you can redistribute it and/or modify  *
  14  *   it under the terms of the GNU General Public License as published by  *
  15  *   the Free Software Foundation; either version 2 of the License, or     *
  16  *   (at your option) any later version.                                   *
  17  *                                                                         *
  18  ***************************************************************************/
  19
  20 #include <qstring.h>
  21
  22 #include "geolocation.h"
  23 #include "timezonerule.h"
  24
  25 GeoLocation::GeoLocation(){
  26         GeoLocation( 0.0, 0.0 );
  27         TZrule = NULL;
  28 }
  29
  30 GeoLocation::GeoLocation( const GeoLocation &g ) {
  31         Longitude = g.Longitude;
  32         Latitude  = g.Latitude;
  33         Name      = g.Name;
  34         Province  = g.Province;
  35         Country   = g.Country;
  36         TimeZone  = g.TimeZone;
  37         TZrule    = g.TZrule;
  38         Height    = g.Height;
  39         indexEllipsoid = g.indexEllipsoid;
  40         setEllipsoid ( indexEllipsoid );
  41         geodToCart();
  42 }
  43
  44 GeoLocation::GeoLocation( GeoLocation *g ) {
  45         Longitude = g->Longitude;
  46         Latitude  = g->Latitude;
  47         Name      = g->Name;
  48         Province  = g->Province;
  49         Country   = g->Country;
  50         TimeZone  = g->TimeZone;
  51         TZrule    = g->TZrule;
  52         Height    = g->Height;
  53         indexEllipsoid = g->indexEllipsoid;
  54         setEllipsoid ( indexEllipsoid );
  55         geodToCart();
  56 }
  57
  58 GeoLocation::GeoLocation( dms lng, dms lat,
  59                                 QString name, QString province, QString country, double tz, TimeZoneRule *tzrule, int iEllips, double hght ) {
  60         Longitude = lng;
  61         Latitude = lat;
  62         Name = name;
  63         Province = province;
  64         Country = country;
  65         TimeZone = tz;
  66         TZrule = tzrule;
  67         Height = hght;
  68         indexEllipsoid = iEllips;
  69         setEllipsoid ( indexEllipsoid );
  70         geodToCart();
  71 }
  72
  73 GeoLocation::GeoLocation( double lng, double lat,
  74                                 QString name, QString province, QString country, double tz, TimeZoneRule *tzrule, int iEllips, double hght ) {
  75         Longitude.set( lng );
  76         Latitude.set( lat );
  77         Name = name;
  78         Province = province;
  79         Country = country;
  80         TimeZone = tz;
  81         TZrule = tzrule;
  82         Height = hght;
  83         indexEllipsoid = iEllips;
  84         setEllipsoid ( indexEllipsoid );
  85         geodToCart();
  86 }
  87
  88 GeoLocation::GeoLocation( double x, double y, double z, QString name, QString province, QString country, double TZ, TimeZoneRule *tzrule, int iEllips ) {
  89         PosCartX = x;
  90         PosCartY = y;
  91         PosCartZ = z;
  92         Name = name;
  93         Province = province;
  94         Country = country;
  95         TimeZone = TZ;
  96         TZrule = tzrule;
  97         indexEllipsoid = iEllips;
  98         setEllipsoid ( indexEllipsoid );
  99         cartToGeod();
 100 }
 101
 102 QString GeoLocation::fullName() const {
 103         QString s;
 104         if ( province().isEmpty() ) {
 105                 s = translatedName() + ", " + translatedCountry();
 106         } else {
 107                 s = translatedName() + ", " + translatedProvince() + ", " + translatedCountry();
 108         }
 109
 110         return s;
 111 }
 112
 113 void GeoLocation::reset( GeoLocation *g ) {
 114         indexEllipsoid = g->ellipsoid();
 115         setEllipsoid ( indexEllipsoid );
 116         setLong( g->lng()->Degrees() );
 117         setLat( g->lat()->Degrees() );
 118         Name      = g->name();
 119         Province  = g->province();
 120         Country   = g->country();
 121         TimeZone  = g->TZ();
 122         TZrule    = g->tzrule();
 123         Height    = g->height();
 124 }
 125
 126
 127 void GeoLocation::setEllipsoid(int index) {
 128         static const double A[] = { 6378140.0, 6378137.0, 6378137.0, 6378137.0, 6378136.0 };
 129         static const double F[] = { 0.0033528131779, 0.0033528106812, 0.0033528131779, 0.00335281066474, 0.0033528131779 };
 130
 131         axis = A[index];
 132         flattening = F[index];
 133 }
 134
 135 void GeoLocation::changeEllipsoid(int index) {
 136
 137         setEllipsoid(index);
 138         cartToGeod();
 139
 140 }
 141
 142 void GeoLocation::cartToGeod(void)
 143 {
 144         static const double RIT = 2.7778e-6;
 145         double e2, rpro, lat1, xn, s1, sqrtP2, latd, sinl;
 146
 147         e2 = 2*flattening-flattening*flattening;
 148
 149         sqrtP2 = sqrt(PosCartX*PosCartX+PosCartY*PosCartY);
 150
 151         rpro = PosCartZ/sqrtP2;
 152         latd = atan(rpro/(1-e2));
 153         lat1 = 0.;
 154
 155         while ( fabs( latd-lat1 ) > RIT ) {
 156                 lat1 = latd;
 157                 s1 = sin(lat1);
 158                 xn = axis/(sqrt(1-e2*s1*s1));
 159                 latd = atan( rpro*(1+e2*xn*s1/PosCartZ) );
 160         }
 161
 162         sinl = sin(latd);
 163         xn = axis/( sqrt(1-e2*sinl*sinl) );
 164
 165         Height = sqrtP2/cos(latd)-xn;
 166         Longitude.setRadians( atan2(PosCartY,PosCartX) );
 167         Latitude.setRadians(latd);
 168 }
 169
 170 void GeoLocation::geodToCart (void) {
 171         double e2, xn;
 172         double sinLong, cosLong, sinLat, cosLat;
 173
 174         e2 = 2*flattening-flattening*flattening;
 175
 176         Longitude.SinCos(sinLong,cosLong);
 177         Latitude.SinCos(sinLat,cosLat);
 178
 179         xn = axis/( sqrt(1-e2*sinLat*sinLat) );
 180         PosCartX = (xn+Height)*cosLat*cosLong;
 181         PosCartY = (xn+Height)*cosLat*sinLong;
 182         PosCartZ = (xn*(1-e2)+Height)*sinLat;
 183 }
 184
 185 void GeoLocation::TopocentricVelocity(double vtopo[], dms gst) {
 186
 187         double Wearth = 7.29211510e-5;     // rads/s
 188         dms angularVEarth;
 189
 190         dms time= GSTtoLST(gst);
 191         // angularVEarth.setRadians(time.Hours()*Wearth*3600.);
 192         double se, ce;
 193         // angularVEarth.SinCos(se,ce);
 194         time.SinCos(se,ce);
 195
 196         double d0 = sqrt(PosCartX*PosCartX+PosCartY*PosCartY);
 197         // km/s
 198         vtopo[0] = - d0 * Wearth * se /1000.;
 199         vtopo[1] = d0 * Wearth * ce /1000.;
 200         vtopo[2] = 0.;
 201 }