.. _program_listing_file_environment_plasma_gcpm_conversions.cc: Program Listing for File conversions.cc ======================================= |exhale_lsh| :ref:`Return to documentation for file ` (``environment/plasma/gcpm/conversions.cc``) .. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS .. code-block:: cpp #include "lupnt/environment/plasma/gcpm/conversions.h" #include "lupnt/environment/plasma/core/math_utils.h" #include "lupnt/environment/plasma/gcpm/constants_gcpm.h" namespace pecsim { void pol_to_cart(double lat, double lon, double r, std::array& pos_sm) { double clat = cos(lat); pos_sm[0] = r * clat * cos(lon); pos_sm[1] = r * clat * sin(lon); pos_sm[2] = r * sin(lat); } void cart_to_pol(const std::array& pos_geo, double& blatr, double& blong, double& rtemp) { double x = pos_geo[0]; double y = pos_geo[1]; double z = pos_geo[2]; double row = sqrt(x * x + y * y); rtemp = sqrt(x * x + y * y + z * z); blatr = atan2(z, row); blong = atan2(y, x); blong = blong + (1.0 - sign(1.0, blong)) * PI; } // transformation from SM to GEO coordinates // SM coordinates are the coordinates defined in the Solar Magnetic (SM) // coordinate system and used in the GCPM model Geographic corrdinates are the // coordinates defined in ECEF and are inputs for the IRI model void sm_to_geo(const std::array& itime, const std::array& x_in, std::array& x_out) { std::array temp1, temp2, temp3; t4(itime, x_in, temp1, -1); t3(itime, temp1, temp2, -1); t2(itime, temp2, temp3, -1); t1(itime, temp3, x_out, 1); } void geo_to_sm(const std::array& itime, const std::array& x_in, std::array& x_out) { std::array temp1, temp2, temp3; t1(itime, x_in, temp1, -1); // GEO -> GEI t2(itime, temp1, temp2, 1); // GEI -> GSE t3(itime, temp2, temp3, 1); // GSE -> GSM t4(itime, temp3, x_out, 1); // GSM -> SM } void gei_to_sm(const std::array& itime, const std::array& x_in, std::array& x_out) { std::array temp1, temp2; t2(itime, x_in, temp1, 1); // GEI -> GSE t3(itime, temp1, temp2, 1); // GSE -> GSM t4(itime, temp2, x_out, 1); // GSM -> SM } // function subroutine for computing time in Julian centuries double to_mjd(const std::array& itime, double& iyr, double& iday, double& ut) { iyr = itime[0] / 1000; iday = itime[0] - iyr * 1000; ut = itime[1] / 1000.0 / 3600.0; // UT in hours double fracday = ut / 24.0; double rmjd = 45.0 + double(iyr - 1859) * 365.0 + (float((iyr - 1861) / 4) + 1.0) + double(iday) - 1.0 + fracday; double T0 = (rmjd - 51544.5) / 36525.0; return T0; } // transformation matrix t1= // GEI -> GEO // Rotates by the Greenwich sidereal time (θ) about the z‑axis to “lock” the // inertial frame to the rotating Earth. void t1(const std::array& itime, const std::array& x_in, std::array& x_out, int iverse) { double iyr, iday, ut; double temp = to_mjd(itime, iyr, iday, ut); double theta = (100.461 + 36000.770 * temp + 15.04107 * ut) * DEG2RAD; rotate_z(double(iverse) * theta, x_in, x_out); } // transformation matrix t2=* // GEI -> GSE // Applies rotations by the obliquity (ε) and the Sun’s ecliptic longitude (λ) // to reorient the inertial frame to one where the x‑axis points to the Sun and // the z‑axis is normal to the ecliptic. void t2(const std::array& itime, const std::array& x_in, std::array& x_out, int iverse) { double iyr, iday, ut; double tt0 = to_mjd(itime, iyr, iday, ut); double epsilon = (23.439 - 0.013 * tt0) * DEG2RAD; double m = (357.528 + 35999.05 * tt0 + 0.04107 * ut) * DEG2RAD; double cgamma = 280.46 + 36000.772 * tt0 + 0.04107 * ut; double lamdas = (cgamma + (1.915 - 0.0048 * tt0) * sin(m) + 0.02 * sin(2.0 * m)) * DEG2RAD; std::array temp; if (iverse == 1) { rotate_x(epsilon, x_in, temp); rotate_z(lamdas, temp, x_out); } else { rotate_z(double(iverse) * lamdas, x_in, temp); rotate_x(double(iverse) * epsilon, temp, x_out); } } // rotation matrix t3=<-psi,X> // GSE -> GSM // Rotates about the x‑axis by an angle (ψ) so that the magnetic dipole // projection falls into the desired plane. void t3(const std::array& itime, const std::array& x_in, std::array& x_out, int iverse) { double iyr, iday, ut; std::array q_c; get_q_c(itime, q_c); double psi = 0.0; if (q_c[2] == 0.0) { psi = -sign(1.5707963, q_c[1]); } else { psi = -atan(q_c[1] / q_c[2]); } rotate_x(double(iverse) * psi, x_in, x_out); } // transformation matrix t4=<-mu,Y> // GSM ->SM // Rotates about the y‑axis by an angle (μ) to align the magnetic dipole with // the z‑axis of the SM frame. void t4(const std::array& itime, const std::array& x_in, std::array& x_out, int iverse) { double iyr, iday, ut; std::array q_c; get_q_c(itime, q_c); double mu = -atan(q_c[0] / sqrt(q_c[1] * q_c[1] + q_c[2] * q_c[2])); rotate_y(double(iverse) * mu, x_in, x_out); } // transformation matrix t5=* // GEO -> MAG // Rotates the geographic (Earth-fixed) coordinates by a geomagnetic longitude // (λ) and a latitude offset (φ – 90°) to produce a magnetic (geomagnetic) // coordinate system. void t5(const std::array& itime, const std::array& x_in, std::array& x_out, int iverse) { double pid2 = 3.1415927 / 2.0; double iyr = itime[0] / 1000; double iday = itime[0] - iyr * 1000; double ut = itime[1] / 1000.0 / 3600.0; double fracday = ut / 24.0; double rmjd = 45.0 + double(iyr - 1859) * 365.0 + (float((iyr - 1861) / 4) + 1.0) + double(iday) - 1.0 + fracday; double factor = (rmjd - 46066.0) / 365.25; double phi = (78.8 + 4.283e-2 * factor) * DEG2RAD; double lamda = (289.1 - 1.413e-2 * factor) * DEG2RAD; std::array temp; if (iverse == 1) { rotate_z(lamda, x_in, temp); rotate_y(phi - pid2, temp, x_out); } else { rotate_y(iverse * (phi - pid2), x_in, temp); rotate_z(iverse * lamda, temp, x_out); } } void rotate_x(double angle, const std::array& x_in, std::array& x_out) { double c = cos(angle); double s = sin(angle); x_out[0] = x_in[0]; x_out[1] = c * x_in[1] + s * x_in[2]; x_out[2] = -s * x_in[1] + c * x_in[2]; } void rotate_y(double angle, const std::array& x_in, std::array& x_out) { double c = cos(angle); double s = sin(angle); x_out[0] = c * x_in[0] + s * x_in[2]; x_out[1] = x_in[1]; x_out[2] = -s * x_in[0] + c * x_in[2]; } void rotate_z(double angle, const std::array& x_in, std::array& x_out) { double c = cos(angle); double s = sin(angle); x_out[0] = c * x_in[0] + s * x_in[1]; x_out[1] = -s * x_in[0] + c * x_in[1]; x_out[2] = x_in[2]; } // subroutine for computing the dipolar axis in GSE coordinates void get_q_c(const std::array& itime, std::array& q_c) { double iyr = itime[0] / 1000; double iday = itime[0] - iyr * 1000; double ut = itime[1] / 1000.0 / 3600.0; double fracday = ut / 24.0; double rmjd = 45.0 + double(iyr - 1859) * 365.0 + (float((iyr - 1861) / 4) + 1.0) + double(iday) - 1.0 + fracday; double factor = (rmjd - 46066.0) / 365.25; double phi = (78.8 + 4.283e-2 * factor) * DEG2RAD; double lamda = (289.1 - 1.413e-2 * factor) * DEG2RAD; std::array temp, q_g; pol_to_cart(phi, lamda, 1.0, q_g); t1(itime, q_g, temp, -1); // GEO -> GEI t2(itime, temp, q_c, 1); // GEI -> GSE } } // namespace pecsim