From dca9dab25e9cf34cdd08c9bf9116fdf8e6533e36 Mon Sep 17 00:00:00 2001
From: Torsten Mayer-Guerr <mayer-guerr@tugraz.at>
Date: Tue, 22 Aug 2023 09:37:32 +0200
Subject: [PATCH] gnssReceiver: corrected local transformationen (north, east,
 up) for LEOs.

---
 source/gnss/gnssObservation.cpp                      |  6 +++---
 source/gnss/gnssReceiver.h                           |  8 ++++----
 .../gnssReceiverGeneratorLowEarthOrbiter.cpp         | 12 ++++++------
 .../gnssReceiverGeneratorStationNetwork.cpp          |  6 +++---
 4 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/source/gnss/gnssObservation.cpp b/source/gnss/gnssObservation.cpp
index e7de3d66..1575033c 100644
--- a/source/gnss/gnssObservation.cpp
+++ b/source/gnss/gnssObservation.cpp
@@ -45,7 +45,7 @@ static void positionVelocityTime(const GnssReceiver &receiver, const GnssTransmi
 
     // line of sight from transmitter to receiver
     k              = normalize(posRecv - posTrans);
-    kRecv          = receiver.local2antennaFrame(idEpoch).transform(receiver.global2localFrame(idEpoch).transform(rotCrf2Trf.rotate(-k))); // line of sight in receiver antenna system (north, east, up)
+    kRecv          = receiver.global2antennaFrame(idEpoch).transform(rotCrf2Trf.rotate(-k)); // line of sight in receiver antenna system (north, east, up)
     azimutRecv     = kRecv.lambda();
     elevationRecv  = kRecv.phi();
     kTrans         = transmitter.celestial2antennaFrame(idEpoch, timeTrans).transform(k);
@@ -128,7 +128,7 @@ Bool GnssObservation::init(const GnssReceiver &receiver, const GnssTransmitter &
     // phase wind-up
     // Carrier phase wind-up in GPS reflectometry, Georg Beyerle, Springer Verlag 2008
     // -------------------------------------------------------------------------------
-    const Transform3d crf2arfRecv  = receiver.local2antennaFrame(idEpoch) * receiver.global2localFrame(idEpoch) * rotCrf2Trf;
+    const Transform3d crf2arfRecv  = receiver.global2antennaFrame(idEpoch) * rotCrf2Trf;
     const Transform3d crf2arfTrans = transmitter.celestial2antennaFrame(idEpoch, timeTrans);
     const Vector3d Tx = crf2arfRecv.transform(crossProduct(crossProduct(k, crf2arfTrans.inverseTransform(Vector3d(1,0,0))), k));
     const Vector3d Ty = crf2arfRecv.transform(crossProduct(crossProduct(k, crf2arfTrans.inverseTransform(Vector3d(0,1,0))), k));
@@ -299,7 +299,7 @@ void GnssObservationEquation::compute(const GnssObservation &observation, const
                          k, kRecvAnt, kTrans);
     const Double   rDotTrans  = inner(k, velocityTrans)/LIGHT_VELOCITY;
     const Double   rDotRecv   = inner(k, velocityRecv) /LIGHT_VELOCITY;
-    const Vector3d kRecvLocal = receiver_.local2antennaFrame(idEpoch).inverseTransform(kRecvAnt);
+    const Vector3d kRecvLocal = receiver_.global2localFrame(idEpoch).transform(rotCrf2Trf.rotate(-k));
     azimutRecvLocal    = kRecvLocal.lambda();
     elevationRecvLocal = kRecvLocal.phi();
 
diff --git a/source/gnss/gnssReceiver.h b/source/gnss/gnssReceiver.h
index 522224d3..4f11db83 100644
--- a/source/gnss/gnssReceiver.h
+++ b/source/gnss/gnssReceiver.h
@@ -50,7 +50,7 @@ class GnssReceiver : public GnssTransceiver
   std::vector<Double>       clk;
   std::vector<Vector3d>     pos, vel;  // regularized marker pos in global system
   std::vector<Vector3d>     offset;    // pos to ARF in local system
-  std::vector<Transform3d>  global2local, local2antenna;
+  std::vector<Transform3d>  global2local, global2antenna;
   std::vector<GnssTrackPtr> tracks;   // tracking phase observations
   Double                    observationSampling;
   Bool                      integerAmbiguities;
@@ -97,11 +97,11 @@ class GnssReceiver : public GnssTransceiver
   /** @brief velocity in TRF or CRF [m/s]. */
   Vector3d velocity(UInt idEpoch) const {return vel.at(idEpoch);}
 
-  /** @brief Rotation from terrestrial reference frame (TRF) or celestial reference frame (CRF) to local/body frame (north, east, up or vehicle system). */
+  /** @brief Rotation from terrestrial reference frame (TRF) or celestial reference frame (CRF) to local horizont system (north, east, up). */
   Transform3d global2localFrame(UInt idEpoch) const {return global2local.at(idEpoch);}
 
-  /** @brief Rotation from local/body frame to left-handed antenna system (usually north, east, up). */
-  Transform3d local2antennaFrame(UInt idEpoch) const {return local2antenna.at(idEpoch);}
+  /** @brief Rotation from terrestrial reference frame (TRF) or celestial reference frame (CRF) to left-handed antenna system (usually north, east, up). */
+  Transform3d global2antennaFrame(UInt idEpoch) const {return global2antenna.at(idEpoch);}
 
   /** @brief Transformation matrix for observed (composed) types from original transmitted types.
   * E.g. C2DG = C1CG - C1WG + C2WG.
diff --git a/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorLowEarthOrbiter.cpp b/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorLowEarthOrbiter.cpp
index a707ee9c..1902d803 100644
--- a/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorLowEarthOrbiter.cpp
+++ b/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorLowEarthOrbiter.cpp
@@ -129,7 +129,7 @@ void GnssReceiverGeneratorLowEarthOrbiter::init(const std::vector<Time> &times,
         recv->vel.resize(times.size());
         recv->offset.resize(times.size());
         recv->global2local.resize(times.size());
-        recv->local2antenna.resize(times.size());
+        recv->global2antenna.resize(times.size());
 
         OrbitArc      orbit      = InstrumentFile::read(fileNameOrbit);
         StarCameraArc starCamera = InstrumentFile::read(fileNameStarCamera);
@@ -155,11 +155,11 @@ void GnssReceiverGeneratorLowEarthOrbiter::init(const std::vector<Time> &times,
             continue;
           }
 
-          recv->pos.at(idEpoch)           = orbit.at(i).position;
-          recv->vel.at(idEpoch)           = orbit.at(i).velocity;
-          recv->offset.at(idEpoch)        = antenna->position - platform.referencePoint(times.at(idEpoch));
-          recv->global2local.at(idEpoch)  = inverse(starCamera.at(i).rotary);
-          recv->local2antenna.at(idEpoch) = antenna->local2antennaFrame;
+          recv->pos.at(idEpoch)            = orbit.at(i).position;
+          recv->vel.at(idEpoch)            = orbit.at(i).velocity;
+          recv->global2local.at(idEpoch)   = inverse(localNorthEastUp(recv->pos.at(idEpoch), Ellipsoid()));
+          recv->global2antenna.at(idEpoch) = antenna->local2antennaFrame * inverse(starCamera.at(i).rotary);
+          recv->offset.at(idEpoch)         = recv->global2local.at(idEpoch).transform(starCamera.at(i).rotary.rotate(antenna->position - platform.referencePoint(times.at(idEpoch))));
           idEpoch++;
         }
         recv->preprocessingInfo("init()");
diff --git a/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorStationNetwork.cpp b/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorStationNetwork.cpp
index aca3ee61..b3f8df99 100644
--- a/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorStationNetwork.cpp
+++ b/source/gnss/gnssReceiverGenerator/gnssReceiverGeneratorStationNetwork.cpp
@@ -190,14 +190,14 @@ void GnssReceiverGeneratorStationNetwork::init(const std::vector<Time> &times, c
             recv->vel.resize(times.size());
             recv->offset.resize(times.size());
             recv->global2local.resize(times.size(), inverse(localNorthEastUp(recv->platform.approxPosition, Ellipsoid())));
-            recv->local2antenna.resize(times.size());
+            recv->global2antenna.resize(times.size());
             for(UInt idEpoch=0; idEpoch<times.size(); idEpoch++)
             {
               auto antenna = recv->platform.findEquipment<PlatformGnssAntenna>(times.at(idEpoch));
               if(antenna && antenna->antennaDef && antenna->accuracyDef)
               {
-                recv->offset.at(idEpoch)        = antenna->position - recv->platform.referencePoint(times.at(idEpoch));
-                recv->local2antenna.at(idEpoch) = antenna->local2antennaFrame;
+                recv->offset.at(idEpoch)         = antenna->position - recv->platform.referencePoint(times.at(idEpoch));
+                recv->global2antenna.at(idEpoch) = antenna->local2antennaFrame * recv->global2local.at(idEpoch);
               }
               else
                 recv->disable(idEpoch, "missing antenna/accuracy patterns");