geo/tcGeo.h

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 #ifndef _tcGeo_h_
  21 #define _tcGeo_h_
  22
  23 /**
  24  * @file tcGeo.h
  25  * @brief Geographical coordinates.
  26  */
  27
  28 #include <Vector.h>
  29 #include <Matrix.h>
  30
  31 #include <QString>
  32 #include <QObject>
  33
  34 #include <cmath>
  35
  36 /// Use doubles for angles.
  37 typedef double tcGeoAngle;
  38
  39 class tcGeoAngleDMS
  40 {
  41   public:
  42     tcGeoAngleDMS(tcGeoAngle angle, bool latitude)
  43     : latitude(latitude)
  44     {
  45       positive = (angle >= 0.0);
  46       if (!positive)
  47       {
  48         angle = -angle;
  49       }
  50       angle *= 180.0 / M_PI;
  51       degrees = (int)floor(angle);
  52       angle -= degrees;
  53       angle *= 60.0;
  54       arcmins = (int)floor(angle);
  55       angle -= arcmins;
  56       arcsecs = angle * 60.0;
  57     }
  58
  59     operator QString () const
  60     {
  61       return QObject::tr("%1\xB0%2'%3\"%4", "degrees, arcminutes, arcseconds, N/S/E/W")
  62                 .arg(degrees)
  63                 .arg(arcmins, 2, 10, QLatin1Char('0'))
  64                 .arg(arcsecs, 5, 'f', 2, QLatin1Char('0'))
  65                 .arg(latitude
  66                      ? (positive ? QObject::tr("N", "north") : QObject::tr("S", "south"))
  67                      : (positive ? QObject::tr("E", "east")  : QObject::tr("W", "west"))
  68                      );
  69     }
  70
  71     bool latitude;
  72     bool positive;
  73     int degrees;
  74     int arcmins;
  75     double arcsecs;
  76 };
  77
  78 /// Geographical coordinates.
  79 class tcGeo
  80 {
  81   public:
  82
  83     /*
  84      * Constructors + destructor
  85      */
  86
  87     /// Primary constructor.
  88     tcGeo()
  89     : m_longitude(0.0)
  90     , m_latitude(0.0)
  91     {
  92     }
  93
  94     /// Primary constructor.
  95     tcGeo(const tcGeoAngle lon, const tcGeoAngle lat)
  96     : m_longitude(lon)
  97     , m_latitude(lat)
  98     {
  99     }
 100
 101     /// Primary constructor.
 102     explicit tcGeo(const maths::Vector<2,double> lonlat)
 103     : m_longitude(lonlat[0])
 104     , m_latitude(lonlat[1])
 105     {
 106     }
 107
 108     /// Construct from a vector.
 109     template <typename T>
 110     tcGeo(maths::Vector<3,T> vec, bool normalized)
 111     {
 112       if (!normalized)
 113       {
 114         vec.normalize();
 115       }
 116       m_longitude = atan2(vec[0], -vec[1]);
 117       m_latitude = asin(vec[2]);
 118     }
 119
 120     /*
 121      * Conversions
 122      */
 123
 124     /// Convert to a 2d position vector with same values.
 125     operator maths::Vector<2,double> () const
 126     {
 127       return maths::Vector<2,double>(m_longitude, m_latitude);
 128     }
 129
 130     /// Convert to a 3d direction vector.
 131     operator maths::Vector<3,float> () const
 132     {
 133       float z = sinf(m_latitude);
 134       float xy = cosf(m_latitude);
 135       return maths::Vector<3,float>(xy*sinf(m_longitude), -xy*cosf(m_longitude), z);
 136     }
 137
 138     /// Convert to a 3d direction vector.
 139     operator maths::Vector<3,double> () const
 140     {
 141       //return this->operator maths::Matrix<3,double>() * maths::Vector<3,double>(0.0, 0.0, 1.0);
 142       double z = sin(m_latitude);
 143       double xy = cos(m_latitude);
 144       return maths::Vector<3,double>(xy*sin(m_longitude), -xy*cos(m_longitude), z);
 145     }
 146
 147     /// Convert to a rotation matrix.
 148     operator maths::Matrix<3,double> () const
 149     {
 150       return maths::Matrix<3,double>(
 151         maths::RotateRadMatrix44('x', m_latitude-M_PI/2)
 152         *
 153         maths::RotateRadMatrix44('z', -m_longitude)
 154         );
 155       double sina = sin(m_longitude);
 156       double cosa = cos(m_longitude);
 157       double sine = sin(m_latitude-M_PI/2);
 158       double cose = cos(m_latitude-M_PI/2);
 159       /*
 160       [[cos(a),                sin(a),                0],
 161        [-sin(a) cos(-1/2 PI+e),cos(a) cos(-1/2 PI+e), sin(-1/2 PI+e)],
 162        [sin(a) sin(-1/2 PI+e), -cos(a) sin(-1/2 PI+e),cos(-1/2 PI+e)]]
 163       */
 164       return maths::Matrix<3,double>(
 165         // First column
 166         maths::Vector<3,double>(cosa, -sina*cose, sina*sine),
 167         // Second column
 168         maths::Vector<3,double>(sina, cosa*cose, -cosa*sine),
 169         // Third column
 170         maths::Vector<3,double>(0.0, sine, cose)
 171       );
 172     }
 173
 174     /*
 175      * Accessors
 176      */
 177
 178     /// Get the longitude.
 179     tcGeoAngle lon() const
 180     {
 181       return m_longitude;
 182     }
 183     /// Get the azimuth (longitude).
 184     tcGeoAngle azim() const
 185     {
 186       return m_longitude;
 187     }
 188
 189     /// Get the latitude.
 190     tcGeoAngle lat() const
 191     {
 192       return m_latitude;
 193     }
 194     /// Get the elevation (latitude).
 195     tcGeoAngle elev() const
 196     {
 197       return m_latitude;
 198     }
 199
 200     /// Get a textual description.
 201     QString describe() const
 202     {
 203       tcGeoAngleDMS dmsLon(m_longitude, false);
 204       tcGeoAngleDMS dmsLat(m_latitude,  true);
 205       return QString("%1 %2").arg(dmsLat).arg(dmsLon);
 206     }
 207
 208     /*
 209      * Mutators
 210      */
 211
 212     /// Set the longitude.
 213     void setLon(const tcGeoAngle lon)
 214     {
 215       m_longitude = lon;
 216     }
 217     /// Set the azimuth (longitude).
 218     void setAzim(const tcGeoAngle lon)
 219     {
 220       m_longitude = lon;
 221     }
 222
 223     /// Set the latitude.
 224     void setLat(const tcGeoAngle lat)
 225     {
 226       m_latitude = lat;
 227     }
 228     /// Set the elevation (latitude).
 229     void setElev(const tcGeoAngle lat)
 230     {
 231       m_latitude = lat;
 232     }
 233
 234     /// Set the longitude and latitude.
 235     void setLonLat(const tcGeoAngle lon, const tcGeoAngle lat)
 236     {
 237       m_longitude = lon;
 238       m_latitude = lat;
 239     }
 240     /// Set the azimuth (longitude) and elevation (latitude).
 241     void setAzimElev(const tcGeoAngle lon, const tcGeoAngle lat)
 242     {
 243       m_longitude = lon;
 244       m_latitude = lat;
 245     }
 246
 247     /*
 248      * Operators
 249      */
 250
 251     tcGeo operator + (const tcGeo& other) const
 252     {
 253       return tcGeo(m_longitude + other.m_longitude,
 254                    m_latitude  + other.m_latitude);
 255     }
 256     tcGeo operator - (const tcGeo& other) const
 257     {
 258       return tcGeo(m_longitude - other.m_longitude,
 259                    m_latitude  - other.m_latitude);
 260     }
 261     tcGeo operator * (const tcGeoAngle factor) const
 262     {
 263       return tcGeo(m_longitude * factor,
 264                    m_latitude  * factor);
 265     }
 266     tcGeo operator / (const tcGeoAngle factor) const
 267     {
 268       return tcGeo(m_longitude / factor,
 269                    m_latitude  / factor);
 270     }
 271     maths::Vector<2,double> operator / (const tcGeo& other) const
 272     {
 273       return maths::Vector<2,double>(m_longitude / other.m_longitude,
 274                                      m_latitude  / other.m_latitude);
 275     }
 276     tcGeo operator * (const maths::Vector<2,double>& other) const
 277     {
 278       return tcGeo(m_longitude * other[0],
 279                    m_latitude  * other[1]);
 280     }
 281     template <typename T>
 282     maths::Vector<3,T> operator * (const maths::Vector<3,T>& other) const
 283     {
 284       return (maths::Vector<3,T>)((maths::Vector<3,double>)other * operator maths::Matrix<3,double>());
 285       float sinLon = sinf(m_longitude);
 286       float cosLon = cosf(m_longitude);
 287       float sinLat = sinf(m_latitude);
 288       float cosLat = cosf(m_latitude);
 289       return maths::Vector<3,T>(other[0]*cosLon + other[2]*sinLon*cosLat,
 290                                 other[0]*sinLon - other[2]*cosLon*cosLat,
 291                                 other[1]*cosLat + other[2]*sinLat);
 292     }
 293
 294     float angle() const
 295     {
 296       return sqrt(m_longitude*m_longitude + m_latitude*m_latitude);
 297     }
 298
 299   private:
 300
 301     /*
 302      * Variables
 303      */
 304
 305     /** East-West geographical coordinate.
 306      * Units are radians.
 307      */
 308     tcGeoAngle m_longitude;
 309
 310     /** North-South geographical coordinate.
 311      * Units are radians.
 312      */
 313     tcGeoAngle m_latitude;
 314
 315 };
 316
 317 #endif // _tcGeo_h_
 318