navigation.c

Go to the documentation of this file.

/**
\file    navigation.c
\brief   GNSS core 'c' function library: navigation functions.
\author  Glenn D. MacGougan (GDM)
\date    2007-11-29
\since   2005-08-23

\b "LICENSE INFORMATION" \n
Copyright (c) 2007, refer to 'author' doxygen tags \n
All rights reserved. \n

Redistribution and use in source and binary forms, with or without
modification, are permitted provided the following conditions are met: \n

- Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer. \n
- Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer in the
  documentation and/or other materials provided with the distribution. \n
- The name(s) of the contributor(s) may not be used to endorse or promote 
  products derived from this software without specific prior written 
  permission. \n

THIS SOFTWARE IS PROVIDED BY THE CONTRIBUTORS ``AS IS'' AND ANY EXPRESS 
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 
OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 
SUCH DAMAGE.
*/

#include <math.h>
#include "constants.h"
#include "geodesy.h"
#include "navigation.h"

void NAVIGATION_ComputeDerivativesOf_Range_WithRespectTo_XYZ( 
  const double x,     // User X coordinate WGS84 ECEF      [m]
  const double y,     // User Y coordinate WGS84 ECEF      [m]
  const double z,     // User Z coordinate WGS84 ECEF      [m]
  const double satX,  // Satellite X coordinate WGS84 ECEF [m]
  const double satY,  // Satellite Y coordinate WGS84 ECEF [m]
  const double satZ,  // Satellite Z coordinate WGS84 ECEF [m]  
  double* dPdx,       // Derivative of P wrt X             []
  double* dPdy,       // Derivative of P wrt Y             []
  double* dPdz,       // Derivative of P wrt Z             []    
  double* range )     // computed user to satellite range  [m]  
{
  double dx;
  double dy;
  double dz;
         
  dx = satX - x;
  dy = satY - y;
  dz = satZ - z;
  
  *range  = sqrt( dx*dx + dy*dy + dz*dz );
  
  // dPdx = d/dx( sqrt( (x_s-x)^2 + (y_s-y)^2 + (z_s-z)^2 ) = -dx/range
  *dPdx   = -dx / (*range);
  *dPdy   = -dy / (*range); 
  *dPdz   = -dz / (*range);      
}




void NAVIGATION_ComputeDerivativesOf_Range_WithRespectToLatitudeLongitudeHeight( 
  const double latitude,  // User geodetic latitude                   [rad]
  const double longitude, // User geodetic longtiude                  [rad]
  const double height,    // User geodetic height                     [m]
  const double satX,      // Satellite X coordinate WGS84 ECEF        [m]
  const double satY,      // Satellite Y coordinate WGS84 ECEF        [m]
  const double satZ,      // Satellite Z coordinate WGS84 ECEF        [m]  
  double* dlat,           // d(P)/d(lat) but not in units of [m/rad], [m/m]
  double* dlon,           // d(P)/d(lon) but not in units of [m/rad], [m/m]
  double* dhgt,           // d(P)/d(hgt)                              [m/m]
  double* range )         // computed user to satellite range         [m]
{
  double sinlat;
  double coslat; 
  double sinlon;
  double coslon;
  double dx;
  double dy;
  double dz;  
  double x; // User X position WGS84 ECEF [m]
  double y; // User Y position WGS84 ECEF [m]
  double z; // User Z position WGS84 ECEF [m]
  
  sinlat = sin(latitude);
  coslat = cos(latitude);
  sinlon = sin(longitude);
  coslon = cos(longitude);

  GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates(
    GEODESY_REFERENCE_ELLIPSE_WGS84,
    latitude,
    longitude,
    height,
    &x,
    &y, 
    &z );
  
  dx = satX - x;
  dy = satY - y;
  dz = satZ - z;

  // d/dx(P) = d/dx( sqrt( (x_s-x)^2 + (y_s-y)^2 + (z_s-z)^2 ) = -dx/range
  // likewise for d/dy, d/dz
  dx *= -1.0;
  dy *= -1.0;
  dz *= -1.0;

  *range  = sqrt( dx*dx + dy*dy + dz*dz );

  /*
  GEODESY_ComputePrimeVerticalRadiusOfCurvature( GEODESY_REFERENCE_ELLIPSE_WGS84, latitude, &N );
  GEODESY_ComputeMeridianRadiusOfCurvature( GEODESY_REFERENCE_ELLIPSE_WGS84, latitude, &M );  

  The derivatives of range with respect to lat and lon in [rad/m]
  lead to numeric instabilities in matrix inversion. Values in [m/rad] are huge!
  *dlat = (M + hgt) * ( -dx*sinlat*coslon - dy*sinlat*sinlon + dz*coslat ) / (*range);
  *dlon = (N + hgt) * ( -dx*coslat*sinlon + dy*coslat*coslon ) / (*range);  
  to scale dlat into a change in [m/m]
  *dlat /= (M + hgt)
  to scale dlon into a change in [m/m] 
  *dlon /= (N + hgt)*coslat
  
  This is effectively the same as rotating the derivatives with repect to xyz to the local frame
  This method requires a good approximation for the initial position when using this parametrization.
  */
  
  *dlat = ( -dx*sinlat*coslon - dy*sinlat*sinlon + dz*coslat ) / (*range);
  *dlon = ( -dx*sinlon        + dy*coslon                    ) / (*range);  
  *dhgt = (  dx*coslat*coslon + dy*coslat*sinlon + dz*sinlat ) / (*range);  
}







int NAVIGATION_PerformClosedFormPositionSolution_FromPseuodrangeMeasurements(
  double p1,          //!< 1st raw pseudorange measurement [m]
  double p2,          //!< 2nd raw pseudorange measurement [m]
  double p3,          //!< 3rd raw pseudorange measurement [m]
  double p4,          //!< 4th raw pseudorange measurement [m]
  double prc_satclk1, //!< 1st satellite clock corrections for psuedoranges [m]
  double prc_satclk2, //!< 2nd satellite clock corrections for psuedoranges [m]
  double prc_satclk3, //!< 3rd satellite clock corrections for psuedoranges [m]
  double prc_satclk4, //!< 4th satellite clock corrections for psuedoranges [m]
  double x1,          //!< 1st satellite X coordinates, WGS84 ECEF [m]
  double x2,          //!< 2nd satellite X coordinates, WGS84 ECEF [m]
  double x3,          //!< 3rd satellite X coordinates, WGS84 ECEF [m]
  double x4,          //!< 4th satellite X coordinates, WGS84 ECEF [m]
  double y1,          //!< 1st satellite Y coordinates, WGS84 ECEF [m]
  double y2,          //!< 2nd satellite Y coordinates, WGS84 ECEF [m]
  double y3,          //!< 3rd satellite Y coordinates, WGS84 ECEF [m]
  double y4,          //!< 4th satellite Y coordinates, WGS84 ECEF [m]
  double z1,          //!< 1st satellite Z coordinates, WGS84 ECEF [m]
  double z2,          //!< 2nd satellite Z coordinates, WGS84 ECEF [m]
  double z3,          //!< 3rd satellite Z coordinates, WGS84 ECEF [m]
  double z4,          //!< 4th satellite Z coordinates, WGS84 ECEF [m]
  double* latitude,   //!< The computed geodetic latitude [rad]
  double* longitude,  //!< The computed geodetic longitude [rad]
  double* height,     //!< The computed geodetic height [m]
  double* rx_clock_bias //!< The computed receiver clock bias [m]
  )
{
  // difference values
  double x12, x13, x14; // 'xij' = 'xi' - 'xj' [m]
  double y12, y13, y14;
  double z12, z13, z14;
  double p21, p31, p41; // 'pij' = 'pi' - 'pj' [m]

  double k1, k2, k3; // coefficients
  double c1, c2, c3;
  double f1, f2, f3;
  double m1, m2, m3;

  double d1;   // clock bias, solution 1 [m]
  double d2;   // clock bias, solution 2 [m]
  double detA; // determinant of A
  double tmp1;
  double tmp2;
  double tmp3;

  double A[3][3];
  double D[3][3]; // D is A^-1 * detA

  typedef struct
  {
    double x;
    double y;
    double z;    
  } struct_SOLN;

  struct_SOLN s1; 
  struct_SOLN s2;

  // apply the satellite clock corrections
  p1 = p1 + prc_satclk1;
  p2 = p2 + prc_satclk2;
  p3 = p3 + prc_satclk3;
  p4 = p4 + prc_satclk4;  
  
  x12 = x1 - x2;
  x13 = x1 - x3;
  x14 = x1 - x4;

  y12 = y1 - y2;
  y13 = y1 - y3;
  y14 = y1 - y4;
  
  z12 = z1 - z2;
  z13 = z1 - z3;
  z14 = z1 - z4;

  p21 = p2 - p1;
  p31 = p3 - p1;
  p41 = p4 - p1;

  k1 = x12*x12 + y12*y12 + z12*z12 - p21*p21;
  k2 = x13*x13 + y13*y13 + z13*z13 - p31*p31;
  k3 = x14*x14 + y14*y14 + z14*z14 - p41*p41;

  A[0][0] = 2.0*x12;
  A[1][0] = 2.0*x13;
  A[2][0] = 2.0*x14;

  A[0][1] = 2.0*y12;
  A[1][1] = 2.0*y13;
  A[2][1] = 2.0*y14;

  A[0][2] = 2.0*z12;
  A[1][2] = 2.0*z13;
  A[2][2] = 2.0*z14;

  tmp1 = A[1][1]*A[2][2] - A[2][1]*A[1][2];
  tmp2 = A[1][0]*A[2][2] - A[2][0]*A[1][2];
  tmp3 = A[1][0]*A[2][1] - A[2][0]*A[1][1];

  detA = A[0][0]*tmp1 - A[0][1]*tmp2 + A[0][2]*tmp3;

  D[0][0] =  tmp1;
  D[1][0] = -tmp2;
  D[2][0] =  tmp3;

  D[0][1] = -A[0][1]*A[2][2] + A[2][1]*A[0][2];
  D[1][1] =  A[0][0]*A[2][2] - A[2][0]*A[0][2];
  D[2][1] = -A[0][0]*A[2][1] + A[2][0]*A[0][1];

  D[0][2] =  A[0][1]*A[1][2] - A[1][1]*A[0][2];
  D[1][2] = -A[0][0]*A[1][2] + A[1][0]*A[0][2];
  D[2][2] =  A[0][0]*A[1][1] - A[1][0]*A[0][1];

  c1 = (D[0][0]*p21 + D[0][1]*p31 + D[0][2]*p41) * 2.0 / detA;
  c2 = (D[1][0]*p21 + D[1][1]*p31 + D[1][2]*p41) * 2.0 / detA;
  c3 = (D[2][0]*p21 + D[2][1]*p31 + D[2][2]*p41) * 2.0 / detA;

  f1 = (D[0][0]*k1 + D[0][1]*k2 + D[0][2]*k3) / detA;
  f2 = (D[1][0]*k1 + D[1][1]*k2 + D[1][2]*k3) / detA;
  f3 = (D[2][0]*k1 + D[2][1]*k2 + D[2][2]*k3) / detA;

  m1 = c1*c1 + c2*c2 + c3*c3 - 1.0;
  m2 = -2.0*( c1*f1 + c2*f2 + c3*f3 );
  m3 = f1*f1 + f2*f2 + f3*f3;

  tmp1 = m2*m2 - 4.0*m1*m3;
  if( tmp1 < 0 )
  {
    // not good, there is no solution
    return FALSE;
  }

  d1 = ( -m2 + sqrt( tmp1 ) ) * 0.5 / m1;
  d2 = ( -m2 - sqrt( tmp1 ) ) * 0.5 / m1;

  // two solutions
  s1.x = c1*d1 - f1 + x1;
  s1.y = c2*d1 - f2 + y1;
  s1.z = c3*d1 - f3 + z1;

  s2.x = c1*d2 - f1 + x1;
  s2.y = c2*d2 - f2 + y1;
  s2.z = c3*d2 - f3 + z1;

  tmp1 = sqrt( s1.x*s1.x + s1.y*s1.y + s1.z*s1.z );
  tmp2 = sqrt( s2.x*s2.x + s2.y*s2.y + s2.z*s2.z );

  // choose the correct solution based
  // on whichever solutino is closest to 
  // the Earth's surface
  tmp1 = fabs( tmp1 - 6371000.0 );
  tmp2 = fabs( tmp2 - 6371000.0 );

  if( tmp2 < tmp1 )
  {
    s1 = s2;
    d1 = d2;
  }

  *rx_clock_bias = d1;

  GEODESY_ConvertEarthFixedCartesianToGeodeticCurvilinearCoordinates(
    GEODESY_REFERENCE_ELLIPSE_WGS84,
    s1.x, s1.y, s1.z,
    latitude, longitude, height );

  if( *height < -1500.0 || *height > 18000.0 )
  {
    // height is out of the likely range for terrestrial users
    return FALSE;
  }
  
  return TRUE;
}