NR UE: introduce orbit propagation to improve N_UE_TA_adj accuracy

Instead of trying to estimate N_UE_TA_adj with a second degree polynomial,
we introduce orbit propagation to calculate N_UE_TA_adj with higher accuracy.

When we receive the satellite ephemeris data consisting of the satellite position vector
and the satellite velocity vector, we calculate parameters for orbit propagation, assuming
a circular orbit in the plane created by the satellite position vector and the satellite
velocity vector. This is done in the function prepare_ue_sat_ta().

Every millisecond, we use these parameters together with the time since epoch to calculate
the satellite position on that orbit and from this the reound-trip-time between UE and the
satellite, what corresponds to N_UE_TA_adj. This is done in the function apply_ntn_timing_advance().

executables/nr-ue.c

@@ -723,6 +723,7 @@ static inline int get_readBlockSize(uint16_t slot, NR_DL_FRAME_PARMS *fp) {
   return rem_samples + next_slot_first_symbol;
 }
 
+#define SPEED_OF_LIGHT 299792458
 static inline void apply_ntn_timing_advance(PHY_VARS_NR_UE *UE, const NR_DL_FRAME_PARMS *fp, int abs_subframe_tx)
 {
   const fapi_nr_dl_ntn_config_command_pdu *ntn_config_params = &UE->ntn_config_message->ntn_config_params;
@@ -732,21 +733,44 @@ static inline void apply_ntn_timing_advance(PHY_VARS_NR_UE *UE, const NR_DL_FRAM
                                + ntn_config_params->epoch_hfn * 10240;
   const int ms_since_epoch = abs_subframe_tx - abs_subframe_epoch;
 
-  const double total_ta_ms = ntn_config_params->ntn_total_time_advance_ms;
-  const double total_ta_drift = ntn_config_params->ntn_total_time_advance_drift; // µs/s
-  const double total_ta_drift_variant = ntn_config_params->ntn_total_time_advance_drift_variant; // µs/s²
+  const double omega = ntn_config_params->omega;
+  const double cos_wt = cos(omega * ms_since_epoch);
+  const double sin_wt = sin(omega * ms_since_epoch);
 
-  UE->timing_advance = (total_ta_ms
-                     + (total_ta_drift / 1000.0) * (ms_since_epoch / 1000.0)
-                     + (total_ta_drift_variant / 1000.0) * (ms_since_epoch / 1000.0) * (ms_since_epoch / 1000.0))
+  const position_t pos_sat_0 = ntn_config_params->pos_sat_0;
+  const position_t pos_sat_90 = ntn_config_params->pos_sat_90;
+  const position_t pos_sat = {pos_sat_0.X * cos_wt + pos_sat_90.X * sin_wt,
+                              pos_sat_0.Y * cos_wt + pos_sat_90.Y * sin_wt,
+                              pos_sat_0.Z * cos_wt + pos_sat_90.Z * sin_wt};
+
+  position_t pos_ue = {0};
+  get_position_coordinates(UE->Mod_id, &pos_ue);
+
+  // calculate directional vector from SAT to UE
+  const position_t dir_sat_ue = {pos_ue.X - pos_sat.X, pos_ue.Y - pos_sat.Y, pos_ue.Z - pos_sat.Z};
+
+  // calculate distance between SAT and UE
+  const double distance = sqrt(dir_sat_ue.X * dir_sat_ue.X + dir_sat_ue.Y * dir_sat_ue.Y + dir_sat_ue.Z * dir_sat_ue.Z);
+
+  // calculate round-trip-time (factor 2) between SAT and UE in ms (factor 1000)
+  const double N_UE_TA_adj = 2000 * distance / SPEED_OF_LIGHT;
+
+  const double N_common_ta_adj = ntn_config_params->N_common_ta_adj;
+  const double N_common_ta_drift = ntn_config_params->N_common_ta_drift;
+  const double N_common_ta_drift_variant = ntn_config_params->N_common_ta_drift_variant;
+
+  UE->timing_advance = (N_UE_TA_adj + N_common_ta_adj
+                      + N_common_ta_drift * ms_since_epoch / 1e6
+                      + N_common_ta_drift_variant * ((int64_t)ms_since_epoch * ms_since_epoch) / 1e9)
                      * fp->samples_per_subframe;
 
   LOG_D(PHY,
-        "total_ta_ms = %f ms, total_ta_drift = %f µs/s, total_ta_drift_variant = %f µs/s², ms_since_epoch = %d ms, "
+        "N_UE_TA_adj = %f ms, N_common_ta_adj = %f ms, N_common_ta_drift = %f µs/s, N_common_ta_drift_variant = %f µs/s², ms_since_epoch = %d ms, "
         "computed timing_advance = %d samples\n",
-        total_ta_ms,
-        total_ta_drift,
-        total_ta_drift_variant,
+        N_UE_TA_adj,
+        N_common_ta_adj,
+        N_common_ta_drift,
+        N_common_ta_drift_variant,
         ms_since_epoch,
         UE->timing_advance);
 }

executables/position_interface.c

@@ -20,6 +20,7 @@
  */
 
 #include "position_interface.h"
+#include "executables/nr-uesoftmodem.h"
 
 static void read_position_coordinates(char *sectionName, position_t *position)
 {

executables/position_interface.h

@@ -22,11 +22,6 @@
 #ifndef __POSITION_INTERFACE__H__
 #define __POSITION_INTERFACE__H__
 
-#include <executables/nr-softmodem-common.h>
-#include <executables/softmodem-common.h>
-#include "PHY/defs_nr_UE.h"
-#include "executables/nr-uesoftmodem.h"
-
 /* position configuration parameters name */
 #define CONFIG_STRING_POSITION_X "x"
 #define CONFIG_STRING_POSITION_Y "y"

nfapi/open-nFAPI/nfapi/public_inc/fapi_nr_ue_interface.h

@@ -24,6 +24,7 @@
 #include "PHY/impl_defs_top.h"
 #include "PHY/impl_defs_nr.h"
 #include "common/utils/nr/nr_common.h"
+#include "executables/position_interface.h"
 
 #define NFAPI_UE_MAX_NUM_CB 8
 #define NFAPI_MAX_NUM_UL_PDU 255
@@ -566,15 +567,22 @@ typedef struct {
   int epoch_sfn;
   int epoch_subframe;
 
+  // orbital angular velocity in rad/ms
+  double omega;
+  // satellite position at epoch time
+  position_t pos_sat_0;
+  // satellite position at 90° orbit
+  position_t pos_sat_90;
+
+  // N_common_ta_adj represents common round-trip-time between gNB and SAT received in SIB19 (ms)
+  double N_common_ta_adj;
+  // drift rate of common ta in µs/s
+  double N_common_ta_drift;
+  // change rate of common ta drift in µs/s²
+  double N_common_ta_drift_variant;
+
   // cell scheduling offset expressed in terms of 15kHz SCS
   long cell_specific_k_offset;
-
-  // ntn_total_time_advance_ms represents the complete round-trip-time between gNB and UE via SAT
-  double ntn_total_time_advance_ms;
-  // drift rate of ntn_total_time_advance_ms in µs/s
-  double ntn_total_time_advance_drift;
-  // change rate of ntn_total_time_advance_ms drift in µs/s²
-  double ntn_total_time_advance_drift_variant;
 } fapi_nr_dl_ntn_config_command_pdu;
 
 typedef struct {

openair1/SCHED_NR_UE/fapi_nr_ue_l1.c

@@ -98,10 +98,17 @@ static void configure_ntn_params(PHY_VARS_NR_UE *ue, fapi_nr_dl_ntn_config_comma
   ue->ntn_config_message->ntn_config_params.epoch_hfn = ntn_params_message->epoch_hfn;
   ue->ntn_config_message->ntn_config_params.epoch_sfn = ntn_params_message->epoch_sfn;
   ue->ntn_config_message->ntn_config_params.epoch_subframe = ntn_params_message->epoch_subframe;
+
+  ue->ntn_config_message->ntn_config_params.omega = ntn_params_message->omega;
+  ue->ntn_config_message->ntn_config_params.pos_sat_0 = ntn_params_message->pos_sat_0;
+  ue->ntn_config_message->ntn_config_params.pos_sat_90 = ntn_params_message->pos_sat_90;
+
+  ue->ntn_config_message->ntn_config_params.N_common_ta_adj = ntn_params_message->N_common_ta_adj;
+  ue->ntn_config_message->ntn_config_params.N_common_ta_drift = ntn_params_message->N_common_ta_drift;
+  ue->ntn_config_message->ntn_config_params.N_common_ta_drift_variant = ntn_params_message->N_common_ta_drift_variant;
+
   ue->ntn_config_message->ntn_config_params.cell_specific_k_offset = ntn_params_message->cell_specific_k_offset;
-  ue->ntn_config_message->ntn_config_params.ntn_total_time_advance_ms = ntn_params_message->ntn_total_time_advance_ms;
-  ue->ntn_config_message->ntn_config_params.ntn_total_time_advance_drift = ntn_params_message->ntn_total_time_advance_drift;
-  ue->ntn_config_message->ntn_config_params.ntn_total_time_advance_drift_variant = ntn_params_message->ntn_total_time_advance_drift_variant;
+
   ue->ntn_config_message->update = true;
 }
 

openair2/LAYER2/NR_MAC_UE/config_ue.c

@@ -31,7 +31,6 @@
  */
 
 #define _GNU_SOURCE
-#define SPEED_OF_LIGHT 299792458
 
 #include "mac_defs.h"
 #include "NR_MAC_UE/mac_proto.h"
@@ -41,7 +40,6 @@
 #include "executables/softmodem-common.h"
 #include "SCHED_NR/phy_frame_config_nr.h"
 #include "oai_asn1.h"
-#include "executables/position_interface.h"
 
 #define ASIGN_P_VAL(dst, src) \
   do {                        \
@@ -266,55 +264,54 @@ static void config_common_ue_sa(NR_UE_MAC_INST_t *mac, NR_ServingCellConfigCommo
 
 }
 
-// computes round-trip-time between ue and sat based on SIB19 ephemeris data
-static void calculate_ue_sat_ta(const position_t *position_params,
-                                const NR_PositionVelocity_r17_t *sat_pos,
-                                ntn_timing_advance_componets_t *ntn_ta)
+// prepare data for orbit propagation based on SIB19 ephemeris data to be able to compute the round-trip-time between ue and sat
+static void prepare_ue_sat_ta(const NR_PositionVelocity_r17_t *sat_pos, ntn_timing_advance_componets_t *ntn_ta)
 {
-  // get UE position coordinates
-  const position_t pos_ue = *position_params;
-
   // get sat position coordinates
   const position_t pos_sat = {sat_pos->positionX_r17 * 1.3, sat_pos->positionY_r17 * 1.3, sat_pos->positionZ_r17 * 1.3};
 
-  // calculate directional vector from SAT to UE
-  const position_t dir_sat_ue = {pos_ue.X - pos_sat.X, pos_ue.Y - pos_sat.Y, pos_ue.Z - pos_sat.Z};
-
-  // calculate distance between SAT and UE
-  double distance = sqrt(dir_sat_ue.X * dir_sat_ue.X + dir_sat_ue.Y * dir_sat_ue.Y + dir_sat_ue.Z * dir_sat_ue.Z);
-
   // get sat velocity vector
   const position_t vel_sat = {sat_pos->velocityVX_r17 * 0.06, sat_pos->velocityVY_r17 * 0.06, sat_pos->velocityVZ_r17 * 0.06};
 
-  // calculate SAT velocity towards UE
-  double velocity = (vel_sat.X * dir_sat_ue.X + vel_sat.Y * dir_sat_ue.Y + vel_sat.Z * dir_sat_ue.Z) / distance;
-
-  // calculate SAT acceleration towards UE
+  // calculate orbital radius
   const double radius_2 = pos_sat.X * pos_sat.X + pos_sat.Y * pos_sat.Y + pos_sat.Z * pos_sat.Z;
+  const double radius = sqrt(radius_2);
+
+  // calculate sat velocity magnitude
   const double vel_sat_2 = vel_sat.X * vel_sat.X + vel_sat.Y * vel_sat.Y + vel_sat.Z * vel_sat.Z;
-  const double acceleration =
-      -(pos_sat.X * dir_sat_ue.X + pos_sat.Y * dir_sat_ue.Y + pos_sat.Z * dir_sat_ue.Z) * vel_sat_2 / (radius_2 * distance);
+  const double vel_mag = sqrt(vel_sat_2);
 
-  LOG_D(NR_MAC, "Satellite velocity towards UE: %f m/s, acceleration towards UE: %f m/s²\n", velocity, acceleration);
+  // calculate angular velocity in rad/ms
+  const double omega = vel_mag / (radius * 1000);
 
-  ntn_ta->N_UE_TA_adj = (2 * distance / SPEED_OF_LIGHT) * 1e3; // in ms
-  ntn_ta->N_UE_TA_drift = (2 * -velocity / SPEED_OF_LIGHT) * 1e6; // in µs/s
-  ntn_ta->N_UE_TA_drift_variant = (2 * acceleration / SPEED_OF_LIGHT) * 1e6; // in µs/s²
+  // calculate sat position in 90° orbit
+  position_t pos_sat_90 = pos_sat;
+  if (vel_mag) {
+    const double scaling = radius / vel_mag;
+    pos_sat_90 = (position_t){vel_sat.X * scaling, vel_sat.Y * scaling, vel_sat.Z * scaling};
+  }
+
+  LOG_I(NR_MAC,
+        "Satellite angular velocity = %e rad/ms, sat_pos = {%f, %f, %f}, sat_pos_90 = {%f, %f, %f}\n",
+        omega,
+        pos_sat.X,
+        pos_sat.Y,
+        pos_sat.Z,
+        pos_sat_90.X,
+        pos_sat_90.Y,
+        pos_sat_90.Z);
+
+  ntn_ta->omega = omega;
+  ntn_ta->pos_sat_0 = pos_sat;
+  ntn_ta->pos_sat_90 = pos_sat_90;
 }
 
 // populate ntn_ta structure from mac
-static void configure_ntn_ta(module_id_t module_id,
-                             ntn_timing_advance_componets_t *ntn_ta,
-                             const NR_NTN_Config_r17_t *ntn_Config_r17,
-                             int hfn,
-                             int frame)
+static void configure_ntn_ta(ntn_timing_advance_componets_t *ntn_ta, const NR_NTN_Config_r17_t *ntn_Config_r17, int hfn, int frame)
 {
   if (!ntn_Config_r17)
     return;
 
-  position_t position_params = {0};
-  get_position_coordinates(module_id, &position_params);
-
   // epochTime_r17 must be present (this is assured by function `eval_epoch_time()`)
   const NR_EpochTime_r17_t *epoch_time_r17 = ntn_Config_r17->epochTime_r17;
   AssertFatal(epoch_time_r17, "epoch_time_r17 should not be NULL here\n");
@@ -331,17 +328,17 @@ static void configure_ntn_ta(module_id_t module_id,
     if (ephemeris_info->present == NR_EphemerisInfo_r17_PR_positionVelocity_r17) {
       const NR_PositionVelocity_r17_t *position_velocity = ephemeris_info->choice.positionVelocity_r17;
       AssertFatal(position_velocity, "position_velocity should not be NULL here\n");
-      calculate_ue_sat_ta(&position_params, position_velocity, ntn_ta);
+      prepare_ue_sat_ta(position_velocity, ntn_ta);
     } else {
       LOG_W(NR_MAC, "NR UE currently supports only ephemerisInfo_r17 of type positionVelocity_r17\n");
-      ntn_ta->N_UE_TA_adj = 0;
-      ntn_ta->N_UE_TA_drift = 0;
-      ntn_ta->N_UE_TA_drift_variant = 0;
+      ntn_ta->omega = 0;
+      ntn_ta->pos_sat_0 = (position_t){0, 0, 0};
+      ntn_ta->pos_sat_90 = (position_t){0, 0, 0};
     }
   } else { // Need R - Release if not present
-    ntn_ta->N_UE_TA_adj = 0;
-    ntn_ta->N_UE_TA_drift = 0;
-    ntn_ta->N_UE_TA_drift_variant = 0;
+    ntn_ta->omega = 0;
+    ntn_ta->pos_sat_0 = (position_t){0, 0, 0};
+    ntn_ta->pos_sat_90 = (position_t){0, 0, 0};
   }
 
   // handle cellSpecificKoffset_r17
@@ -374,18 +371,14 @@ static void configure_ntn_ta(module_id_t module_id,
 
   LOG_I(NR_MAC,
         "SIB19 Rxd. Epoch HFN: %d, Epoch SFN: %d, Epoch Subframe: %d, k_offset: %ldms, "
-        "N_Common_Ta: %fms, drift: %fµs/s, variant %fµs/s², "
-        "N_UE_TA: %fms, drift: %fµs/s, variant %fµs/s²\n",
+        "N_Common_Ta: %fms, drift: %fµs/s, variant %fµs/s²\n",
         ntn_ta->epoch_hfn,
         ntn_ta->epoch_sfn,
         ntn_ta->epoch_subframe,
         ntn_ta->cell_specific_k_offset,
         ntn_ta->N_common_ta_adj,
         ntn_ta->N_common_ta_drift,
-        ntn_ta->N_common_ta_drift_variant,
-        ntn_ta->N_UE_TA_adj,
-        ntn_ta->N_UE_TA_drift,
-        ntn_ta->N_UE_TA_drift_variant);
+        ntn_ta->N_common_ta_drift_variant);
 }
 
 static void config_common_ue(NR_UE_MAC_INST_t *mac, NR_ServingCellConfigCommon_t *scc, int cc_idP, int hfn, int frame)
@@ -563,7 +556,7 @@ static void config_common_ue(NR_UE_MAC_INST_t *mac, NR_ServingCellConfigCommon_t
   // NTN Config
   if (scc->ext2) {
     UPDATE_IE(mac->sc_info.ntn_Config_r17, scc->ext2->ntn_Config_r17, NR_NTN_Config_r17_t);
-    configure_ntn_ta(mac->ue_id, &mac->ntn_ta, mac->sc_info.ntn_Config_r17, hfn, frame);
+    configure_ntn_ta(&mac->ntn_ta, mac->sc_info.ntn_Config_r17, hfn, frame);
   } else {
     asn1cFreeStruc(asn_DEF_NR_NTN_Config_r17, mac->sc_info.ntn_Config_r17);
   }
@@ -1964,7 +1957,7 @@ void nr_rrc_mac_config_other_sib(module_id_t module_id, NR_SIB19_r17_t *sib19, i
   if (sib19) {
     // update ntn_Config_r17 with received values
     UPDATE_IE(mac->sc_info.ntn_Config_r17, sib19->ntn_Config_r17, NR_NTN_Config_r17_t);
-    configure_ntn_ta(mac->ue_id, &mac->ntn_ta, mac->sc_info.ntn_Config_r17, hfn, frame);
+    configure_ntn_ta(&mac->ntn_ta, mac->sc_info.ntn_Config_r17, hfn, frame);
   }
   if (mac->state == UE_RECEIVING_SIB && can_start_ra)
     mac->state = UE_PERFORMING_RA;

openair2/LAYER2/NR_MAC_UE/mac_defs.h

@@ -62,6 +62,8 @@
 #include "NR_MeasConfig.h"
 #include "NR_ServingCellConfigCommonSIB.h"
 
+/* position_t */
+#include "executables/position_interface.h"
 
 // ==========
 // NR UE defs
@@ -547,18 +549,20 @@ typedef struct ntn_timing_advance_components {
   int epoch_sfn;
   int epoch_subframe;
 
+  // orbital angular velocity in rad/ms
+  double omega;
+  // satellite position at epoch time
+  position_t pos_sat_0;
+  // satellite position at 90° orbit
+  position_t pos_sat_90;
+
   // N_common_ta_adj represents common round-trip-time between gNB and SAT received in SIB19 (ms)
   double N_common_ta_adj;
   // drift rate of common ta in µs/s
   double N_common_ta_drift;
   // change rate of common ta drift in µs/s²
   double N_common_ta_drift_variant;
-  // N_UE_TA_adj calculated round-trip-time between UE and SAT (ms)
-  double N_UE_TA_adj;
-  // drift rate of N_UE_TA in µs/s
-  double N_UE_TA_drift;
-  // change rate of N_UE_TA drift in µs/s²
-  double N_UE_TA_drift_variant;
+
   // cell scheduling offset expressed in terms of 15kHz SCS
   long cell_specific_k_offset;
 
@@ -681,17 +685,7 @@ static inline long GET_DURATION_RX_TO_TX(const ntn_timing_advance_componets_t *n
 
 static inline double get_total_TA_ms(const ntn_timing_advance_componets_t *ntn_ta)
 {
-  return ntn_ta->N_common_ta_adj + ntn_ta->N_UE_TA_adj;
-}
-
-static inline double get_total_TA_drift(const ntn_timing_advance_componets_t *ntn_ta)
-{
-  return ntn_ta->N_common_ta_drift + ntn_ta->N_UE_TA_drift;
-}
-
-static inline double get_total_TA_drift_variant(const ntn_timing_advance_componets_t *ntn_ta)
-{
-  return ntn_ta->N_common_ta_drift_variant + ntn_ta->N_UE_TA_drift_variant;
+  return ntn_ta->N_common_ta_adj;
 }
 
 /*@}*/

openair2/LAYER2/NR_MAC_UE/nr_ue_scheduler.c

@@ -2852,10 +2852,16 @@ static void schedule_ntn_config_command(fapi_nr_dl_config_request_t *dl_config,
   ntn_config_command_pdu->epoch_hfn = mac->ntn_ta.epoch_hfn;
   ntn_config_command_pdu->epoch_sfn = mac->ntn_ta.epoch_sfn;
   ntn_config_command_pdu->epoch_subframe = mac->ntn_ta.epoch_subframe;
+
+  ntn_config_command_pdu->omega = mac->ntn_ta.omega;
+  ntn_config_command_pdu->pos_sat_0 = mac->ntn_ta.pos_sat_0;
+  ntn_config_command_pdu->pos_sat_90 = mac->ntn_ta.pos_sat_90;
+
+  ntn_config_command_pdu->N_common_ta_adj = mac->ntn_ta.N_common_ta_adj;
+  ntn_config_command_pdu->N_common_ta_drift = mac->ntn_ta.N_common_ta_drift;
+  ntn_config_command_pdu->N_common_ta_drift_variant = mac->ntn_ta.N_common_ta_drift_variant;
+
   ntn_config_command_pdu->cell_specific_k_offset = mac->ntn_ta.cell_specific_k_offset;
-  ntn_config_command_pdu->ntn_total_time_advance_ms = get_total_TA_ms(&mac->ntn_ta);
-  ntn_config_command_pdu->ntn_total_time_advance_drift = get_total_TA_drift(&mac->ntn_ta);
-  ntn_config_command_pdu->ntn_total_time_advance_drift_variant = get_total_TA_drift_variant(&mac->ntn_ta);
 
   dl_config->dl_config_list[dl_config->number_pdus].pdu_type = FAPI_NR_DL_NTN_CONFIG_PARAMS;
   dl_config->number_pdus += 1;