Program Listing for File raytrace.h

Return to documentation for file (environment/plasma/tec/raytrace.h)

#pragma once

#include <functional>
#include <limits>
#include <vector>

#include "lupnt/environment/plasma/core/definitions.h"
#include "lupnt/environment/plasma/gcpm/constants_gcpm.h"
#include "lupnt/environment/plasma/gcpm/conversions.h"
#include "lupnt/environment/plasma/gcpm/gcpm_interface.h"

namespace pecsim {

  // Constants
  const double freq_L1 = 1575.42e6;  // L1 frequency
  const double freq_L2 = 1227.60e6;  // L2 frequency
  const double freq_L5 = 1176.45e6;  // L5 frequency

  struct RayTraceConfig {
    double freq_Hz = std::numeric_limits<double>::quiet_NaN();  // Frequency in Hz
    double step_size = 10.0;                                    // Step size for ray tracing [km]
    bool correction = true;                  // Apply correction to the ray tracing
    bool fine_correction = false;            // Apply fine correction to the ray tracing
    double cutoff_r = 4 * RE;                // Cutoff radius for ray tracing [km]
    double gradn_dx = 1.0;                   // Gradient step size for refractive index [km]
    std::string integ_method = "Euler";      // Integration method ("RK4" or "Euler")
    std::string correction_method = "grid";  // Correction method ("grid" or "newton")
    double kp = -1;                          // Kp index for the ionosphere model
    bool use_fortran_gcpm = true;            // Use Fortran for ray tracing
    double corr_tol = 1.0;                   // Correction tolerance [m] for ray tracing
    bool compute_higher_order = true;        // Compute second-order delays
    bool use_adaptive_step = true;           // Use adaptive step size for ray tracing
    bool straight_ray = false;               // Assume straight ray path
  };

  struct PathProfile {
    // per each timestep
    VecXd s;             // Path length
    VecXd tec_section;   // Total Electron density for each section [TECU]
    VecXd r;             // Geocentric radial distance
    MatXd pos_eci;       // Position in Cartesian coordinates
    VecXd az_dir;        // Azimuth direction
    VecXd el_dir;        // Elevation direction
    VecXd dist_to_line;  // Distance to the straight line

    // final values
    double sf;                   // Final distance of the ray [km]
    Vec3d dir_start;             // Direction vector at the start
    Vec3d dir_end;               // Direction vector at the end
    Vec3d final_pos;             // Final position in Cartesian coordinate
    Vec3d corr_final_pos_err;    // Final position error in Cartesian coordinate due
                                 // to insufficient correction
    double corr_final_time_err;  // Final time error [s] due to insufficient
                                 // correction
    double t_tx;                 // Transmission epoch (seconds since J2000)
    double t_rx;                 // Reception epoch (seconds since J2000)
    double prop_time_total;      // Total propagation time [s] (technically t_rx -
                                 // t_tx, but computed with care in float limits)
    double dist_bend_m;          // Distance bent [m]
    double dist_straight_km;     // Straight distance [km]
    double total_delay_m;        // Total delay [m]
    double tec_delay_m;          // TEC delay [m]
    double tec_delay_bend_m;     // additional TEC due to ray bending [m] (part of
                                 // tec_delay_m)
    double second_delay_m;       // Second order delay [m]
    double third_delay_m;        // Third order delay [m] (if computed)
    double tecu;                 // Total electron content [TECU]
    double tecu_bend;            // TEC bend [TECU]
    double max_sep_line_m;       // Maximum separation from the straight line [m]
  };

  bool is_method_rk4(const std::string& method);

  Vec3d azel_to_unitvec(double az, double el);

  Vec2d unitvec_to_azel(const Vec3d& x);

  double compute_ne(double t_j2000, const Vec3d& pos_geo, RayTraceConfig config,
                    bool debug = false);

  Vec3d get_iono_params(double t_j2000, double kp = -1.0);

  Vec3d compute_B(double t_j2000, const Vec3d& pos_geo, RayTraceConfig config, bool debug = false);

  double refractive_index_neB(double ne_m3, double freq_Hz, double B, double costheta,
                              bool compute_higher_order = true);

  double refractive_index(double t, const Vec3d& x, const Vec3d& shat, RayTraceConfig config,
                          bool debug = false);

  Vec3d gradient_n(double t, const Vec3d& x, RayTraceConfig config, bool debug = false);

  VecXd ray_derivative(double s, const VecXd& z, double t_epoch_tx, RayTraceConfig config,
                       VecXd& store_vec, bool use_precomputed_vals, bool debug = false);

  void integ_step(VecXd& z, double& s, double h, double t_epoch_tx, RayTraceConfig config,
                  VecXd& store_vec, bool use_precomputed_vals, bool debug = false);

  double adjust_stepsize(double r, const VecXd& store_vec);

  VecXd propagate_ray(double t_prop, const Vec3d& dir, double sf, const Vec3d& x1, const Vec3d& x2,
                      double t_rx, RayTraceConfig config, MatXd& store_mat,
                      bool use_precomputed_densities = false);

  PathProfile trace_ray(double epoch_utc_rx, const Vec3d& x1, const Vec3d& x2,
                        RayTraceConfig config, bool debug_prop = false, bool debug_corr = false);

  PathProfile compute_path_profile(double t_prop, const Vec3d& dir, double sf, const Vec3d& x1,
                                   const Vec3d& x2, double t_rx, RayTraceConfig config,
                                   MatXd& store_mat, bool use_precomputed_densities = false,
                                   bool debug = false);

  double compute_tec_straight(double t_tx, const Vec3d& initial_pos, const Vec3d& final_pos,
                              RayTraceConfig config);

  double compute_tec_section(int i, double t_tx, const Vec3d& initial_pos, const Vec3d& dir,
                             const RayTraceConfig& config);

  void correct_proptime(double& prop_time, Vec3d& dir, double sf, const Vec3d& x1, const Vec3d& x2,
                        double t_rx, RayTraceConfig config, VecXd& zf);

  void correct_dir_neldermead(double& t_prop, Vec3d& dir, double sf, const Vec3d& x1,
                              const Vec3d& x2, double t_rx, RayTraceConfig config, VecXd& zf,
                              MatXd& store_mat,
                              bool coarse_mode = false,  // Use precomputed dzds
                              bool debug = false);
  void correct_dir_newton(double& t_prop, Vec3d& dir, double sf, const Vec3d& x1, const Vec3d& x2,
                          double t_rx, RayTraceConfig config, VecXd& zf, bool debug = false);

}  // namespace pecsim