diff --git a/.gitignore b/.gitignore index 63b9617..cc838e6 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,11 @@ ## ## Get latest from https://github.com/github/gitignore/blob/master/VisualStudio.gitignore +# Additional +FOORT/cfgs/* +FOORT/Output/* +Benchmarking/Seppe/* + # FOORT compilation files FOORT/*.o FOORT/FOORT diff --git a/FOORT/src/CMakeLists.txt b/FOORT/src/CMakeLists.txt index c560469..b063489 100644 --- a/FOORT/src/CMakeLists.txt +++ b/FOORT/src/CMakeLists.txt @@ -7,6 +7,8 @@ add_library(config_reader_lib STATIC ConfigReader.cpp ConfigReader.h) add_library(diagnostics_lib STATIC Diagnostics.cpp Diagnostics.h Geometry.h) add_library(diagnostics_emission_lib STATIC DiagnosticsEmission.cpp DiagnosticsEmission.h Geometry.h) add_library(geodesic_lib STATIC Geodesic.cpp Geodesic.h Geometry.h) +add_library(grid_lib STATIC Grid.cpp Grid.h) +add_library(interpolator_lib STATIC Interpolator.cpp Interpolator.h) add_library(input_output_lib STATIC InputOutput.cpp InputOutput.h Geometry.h) add_library(integrators_lib STATIC Integrators.cpp Integrators.h Geometry.h) add_library(mesh_lib STATIC Mesh.cpp Mesh.h Geometry.h) @@ -24,7 +26,7 @@ target_link_libraries(geodesic_lib PUBLIC diagnostics_lib input_output_lib integ target_link_libraries(input_output_lib stdc++fs) target_link_libraries(integrators_lib PUBLIC geodesic_lib metric_lib) target_link_libraries(mesh_lib PUBLIC diagnostics_lib input_output_lib utilities_lib) -target_link_libraries(metric_lib PUBLIC input_output_lib integrators_lib spline_lib) +target_link_libraries(metric_lib PUBLIC grid_lib input_output_lib integrators_lib interpolator_lib spline_lib) target_link_libraries(terminations_lib PUBLIC geodesic_lib input_output_lib) target_link_libraries(utilities_lib PUBLIC diagnostics_lib geodesic_lib integrators_lib metric_lib terminations_lib viewscreen_lib) target_link_libraries(viewscreen_lib PUBLIC metric_lib mesh_lib) diff --git a/FOORT/src/Config.cpp b/FOORT/src/Config.cpp index 4dcdb65..3cc9ee3 100644 --- a/FOORT/src/Config.cpp +++ b/FOORT/src/Config.cpp @@ -294,6 +294,28 @@ std::unique_ptr Config::GetMetric(const ConfigCollection &theCfg) // All settings complete; create Metric object! TheMetric = std::unique_ptr(new BosonStarMetric(Phi_infinity, num_lines, rLogScale, Phi_filename, m_filename)); } + else if (MetricName == "rotatingbosonstar") + { + // The rotating boson star with solitonic potential + + // First setting to look up: using a logarithmic r coordinate or not. + // Don't need to output message if setting not found + bool rLogScale{false}; + bool flipAngularMomentum{false}; + std::string MetricFolder{"RotatingBosonStar/data_Will/"}; + int NumX{500}; + int NumTh{399}; + real L{1.}; + MetricSettings.LookupValue("RLogScale", rLogScale); + MetricSettings.LookupValue("FlipAngularMomentum", flipAngularMomentum); + MetricSettings.LookupValue("MetricFolder", MetricFolder); + MetricSettings.LookupValue("NumX", NumX); + MetricSettings.LookupValue("NumTh", NumTh); + MetricSettings.LookupValue("L", L); + + // All settings complete; create Metric object! + TheMetric = std::unique_ptr(new RotatingBosonStarMetric(rLogScale, MetricFolder, NumX, NumTh, L, flipAngularMomentum)); + } //// METRIC ADD POINT B //// // Add an else if clause to check for your new Metric object! // To look for additional options in the metric configuration, use @@ -602,26 +624,32 @@ void Config::InitializeDiagnostics(const ConfigCollection &theCfg, DiagBitflag & //// Fluid four-velocity model selection and initialization //// - // Default fluid model and default parameters - real defaultxi{1.0}; - real defaultbetar{1.0}; - real defaultbetaphi{1.0}; - std::unique_ptr theFluidModel{new GeneralCircularRadialFluid(defaultxi, defaultbetar, defaultbetaphi, theMetric)}; + // Fluid four-velocity model parameters (with sensible defaults) + real subKeplerianparam{1.0}; + real betaR{1.0}; + real betaPhi{1.0}; + real iscolowerbound{0.05}; + real iscoupperbound{1000.0}; // Read in fluid velocity model - std::string fluidmodelstring{""}; + std::string fluidmodelstring{"GeneralCircularRadial"}; AllDiagSettings["EquatorialEmission"].LookupValue("FluidVelocityModel", fluidmodelstring); if (fluidmodelstring == "GeneralCircularRadial") { - real subKeplerianparam{defaultxi}; - real betaR{defaultbetar}; - real betaPhi{defaultbetaphi}; AllDiagSettings["EquatorialEmission"].LookupValue("xi", subKeplerianparam); AllDiagSettings["EquatorialEmission"].LookupValue("betar", betaR); AllDiagSettings["EquatorialEmission"].LookupValue("betaphi", betaPhi); - - theFluidModel = std::unique_ptr{new GeneralCircularRadialFluid(subKeplerianparam, betaR, betaPhi, theMetric)}; + AllDiagSettings["EquatorialEmission"].LookupValue("ISCOLowerBound", iscolowerbound); + AllDiagSettings["EquatorialEmission"].LookupValue("ISCOUpperBound", iscoupperbound); } + else + { + ScreenOutput("Unknown FluidVelocityModel \"" + fluidmodelstring + "\". Using default GeneralCircularRadial model.", + Output_Other_Default); + } + + std::unique_ptr theFluidModel{ + new GeneralCircularRadialFluid(subKeplerianparam, betaR, betaPhi, theMetric, iscolowerbound, iscoupperbound)}; // Other fluid velocity models can be checked for here... // Set EquatorialEmissionDiagnostic options struct diff --git a/FOORT/src/DiagnosticsEmission.cpp b/FOORT/src/DiagnosticsEmission.cpp index 5de7457..cd48882 100644 --- a/FOORT/src/DiagnosticsEmission.cpp +++ b/FOORT/src/DiagnosticsEmission.cpp @@ -180,7 +180,7 @@ std::string GeneralCircularRadialFluid::getFullDescriptionStr() const if (m_ISCOexists && m_theMetric->getrLogScale()) trueISCOradius = exp(m_ISCOr); - return "Circular/radial flow (sub-Keplerian parameter xi = " + std::to_string(m_subKeplerParam) + ", beta_r = " + std::to_string(m_betaR) + ", beta_phi = " + std::to_string(m_betaPhi) + "; " + (m_ISCOexists ? "ISCO = " + std::to_string(trueISCOradius) : "no ISCO found") + ")"; + return "Circular/radial flow (sub-Keplerian parameter xi = " + std::to_string(m_subKeplerParam) + ", beta_r = " + std::to_string(m_betaR) + ", beta_phi = " + std::to_string(m_betaPhi) + "; " + (m_ISCOexists ? "ISCO = " + std::to_string(trueISCOradius) : "no ISCO found") + "; " + "ISCO lower bound = " + std::to_string(m_ISCOlowerbound) + ", ISCO upper bound = " + std::to_string(m_ISCOupperbound) + ")"; } /** @@ -352,14 +352,20 @@ OneIndex GeneralCircularRadialFluid::GetRadialVelocityd(const Point &p) const */ void GeneralCircularRadialFluid::FindISCO() { - real lowerbound{0.0}; - real upperbound{1000.0}; + real lowerbound{m_ISCOlowerbound}; + real upperbound{m_ISCOupperbound}; const SphericalHorizonMetric *sphermetric = dynamic_cast(m_theMetric); if (sphermetric) { lowerbound = sphermetric->getrLogScale() ? log(sphermetric->getHorizonRadius()) : sphermetric->getHorizonRadius(); upperbound = sphermetric->getrLogScale() ? log(10.0 * sphermetric->getHorizonRadius()) : 10.0 * sphermetric->getHorizonRadius(); } + else if (m_theMetric->getrLogScale()) + { + // If not a spherical horizon metric, but we are using log(r) coordinates, then set the lower bound to 0.0 + lowerbound = log(m_ISCOlowerbound); // ln(0.05) = -3.912023005428146 + upperbound = log(m_ISCOupperbound); // ln(1000.0) = 6.907755278982137 + } // Perform binary search for ISCO // Allow only 1000 iterations max diff --git a/FOORT/src/DiagnosticsEmission.h b/FOORT/src/DiagnosticsEmission.h index 8c242ac..513592a 100644 --- a/FOORT/src/DiagnosticsEmission.h +++ b/FOORT/src/DiagnosticsEmission.h @@ -95,8 +95,10 @@ struct FluidVelocityModel struct GeneralCircularRadialFluid final : public FluidVelocityModel { // Constructor with three parameters and Metric pointer (which is passed to base class constructor) - GeneralCircularRadialFluid(real subKeplerParam, real betar, real betaphi, const Metric *const theMetric) : m_subKeplerParam{fmin(fmax(subKeplerParam, 0.0), 1.0)}, m_betaR{fmin(fmax(betar, 0.0), 1.0)}, - m_betaPhi{fmin(fmax(betaphi, 0.0), 1.0)}, FluidVelocityModel(theMetric) + GeneralCircularRadialFluid(real subKeplerParam, real betar, real betaphi, const Metric *const theMetric, + real ISCO_lowerbound, real ISCO_upperbound) : m_subKeplerParam{fmin(fmax(subKeplerParam, 0.0), 1.0)}, m_betaR{fmin(fmax(betar, 0.0), 1.0)}, + m_betaPhi{fmin(fmax(betaphi, 0.0), 1.0)}, FluidVelocityModel(theMetric), + m_ISCOlowerbound{ISCO_lowerbound}, m_ISCOupperbound{ISCO_upperbound} { // Do some checks on three params, which must lie between 0.0 and 1.0 (note that they are adjusted as such in // initializer above) @@ -150,6 +152,10 @@ struct GeneralCircularRadialFluid final : public FluidVelocityModel bool m_ISCOexists{false}; //! ISCO radius real m_ISCOr{-1.0}; + //! Lower bound for ISCO radius search + real m_ISCOlowerbound; + //! Upper bound for ISCO radius search + real m_ISCOupperbound; //! ISCO t momentum real m_ISCOpt{}; //! ISCO phi momentum diff --git a/FOORT/src/Grid.cpp b/FOORT/src/Grid.cpp new file mode 100644 index 0000000..296c810 --- /dev/null +++ b/FOORT/src/Grid.cpp @@ -0,0 +1,46 @@ +#include +#include +#include +#include "Grid.h" + +//! A class for a 2D grid that contains the values of the necessary functions. +Grid::Grid(int N_row, int N_col) +{ + this->N_row = N_row; + this->N_col = N_col; + this->data = new double[N_row * N_col]; + size = N_row * N_col; +} + +void Grid::initialize_from_file(std::string file) +{ + std::ifstream inputFile(file); + + if (!inputFile) + { + std::cerr << "Error opening file!" << std::endl; + } + + std::string line; + int i = 0; + int j = 0; + + while (std::getline(inputFile, line)) + { + std::istringstream iss(line); + + double val; + while (iss >> val) + { + data[i * N_col + j] = val; + j += 1; + } + j = 0; + i += 1; + } + std::cout << "Grid initialized from " << file << std::endl + << std::endl; + + // close the file + inputFile.close(); +} diff --git a/FOORT/src/Grid.h b/FOORT/src/Grid.h new file mode 100644 index 0000000..5f60d3e --- /dev/null +++ b/FOORT/src/Grid.h @@ -0,0 +1,32 @@ +#include +#include + +#ifndef GRID_H +#define GRID_H + +//! A class for a 2D grid that contains the values of the necessary functions. +class Grid +{ +public: + //! Number of rows in the grid + int N_row; + //! Number of columns in the grid + int N_col; + //! Pointer to the data + double *data; + //! Number of rows times the number of columns + int size; + Grid(int N_row, int N_col); + + //! Destructor + ~Grid() { delete[] this->data; } + + //! Overload the () operator to access the data + double &operator()(int i, int j) { return this->data[i * N_col + j]; } + //! Overload the () operator to access the data (const version) + const double &operator()(int i, int j) const { return this->data[i * N_col + j]; } + + void initialize_from_file(std::string file); +}; + +#endif // GRID_H diff --git a/FOORT/src/Interpolator.cpp b/FOORT/src/Interpolator.cpp new file mode 100644 index 0000000..19f780a --- /dev/null +++ b/FOORT/src/Interpolator.cpp @@ -0,0 +1,249 @@ +#include +#include +#include +#include +#include +#include "Interpolator.h" +#include "Grid.h" + +Interpolator::Interpolator(int dim_x, int dim_th, std::string x_file, std::string th_file) +{ + this->m_x = new double[dim_x]; + this->m_theta = new double[dim_th]; + + std::ifstream x_input(x_file); + std::ifstream th_input(th_file); + + if (!th_input || !x_input) + { + std::cerr << "Error opening file!" << std::endl; + } + + for (int i = 0; i < dim_x; ++i) + { + x_input >> this->m_x[i]; + } + for (int i = 0; i < dim_th; ++i) + { + th_input >> this->m_theta[i]; + } +} + +//! Returns the fourth order lagrangian polynomial evaluated at the point of interest +double Interpolator::fourth_order_lagange(double x_points[4], double y_points[4], double x) +{ + double x1{x_points[0]}, x2{x_points[1]}, x3{x_points[2]}, x4{x_points[3]}; + + double value = (x - x2) * (x - x3) * (x - x4) * y_points[0] / ((x1 - x2) * (x1 - x3) * (x1 - x4)) + (x - x1) * (x - x3) * (x - x4) * y_points[1] / ((x2 - x1) * (x2 - x3) * (x2 - x4)) + (x - x1) * (x - x2) * (x - x4) * y_points[2] / ((x3 - x1) * (x3 - x2) * (x3 - x4)) + (x - x1) * (x - x2) * (x - x3) * y_points[3] / ((x4 - x1) * (x4 - x2) * (x4 - x3)); + return value; +} + +//! Interpolates the data of the grid +//! @param grid The grid to be interpolated +double Interpolator::interpolate(Grid *grid, double p_x, double p_theta) +{ + // find indices of theta and x + int x_ind = 0; + int th_ind = 0; + for (int k = 0; k < grid->N_col; ++k) + { + if (p_x <= m_x[k]) + { + x_ind = k - 1; + break; + } + } + for (int k = 0; k < grid->N_row; ++k) + { + if (p_theta <= m_theta[k]) + { + th_ind = k - 1; + break; + } + } + + //! find the interpolated value + //! extract the coordinates at which we have values + double x1 = m_x[x_ind - 1]; + double x2 = m_x[x_ind]; + double x3 = m_x[x_ind + 1]; + double x4 = m_x[x_ind + 2]; + double x[4] = {x1, x2, x3, x4}; + + double y1 = m_theta[th_ind - 1]; + double y2 = m_theta[th_ind]; + double y3 = m_theta[th_ind + 1]; + double y4 = m_theta[th_ind + 2]; + double y[4] = {y1, y2, y3, y4}; + + double vals1[4]{(*grid)(th_ind - 1, x_ind - 1), (*grid)(th_ind - 1, x_ind), (*grid)(th_ind - 1, x_ind + 1), (*grid)(th_ind - 1, x_ind + 2)}; + double vals2[4]{(*grid)(th_ind, x_ind - 1), (*grid)(th_ind, x_ind), (*grid)(th_ind, x_ind + 1), (*grid)(th_ind, x_ind + 2)}; + double vals3[4]{(*grid)(th_ind + 1, x_ind - 1), (*grid)(th_ind + 1, x_ind), (*grid)(th_ind + 1, x_ind + 1), (*grid)(th_ind + 1, x_ind + 2)}; + double vals4[4]{(*grid)(th_ind + 2, x_ind - 1), (*grid)(th_ind + 2, x_ind), (*grid)(th_ind + 2, x_ind + 1), (*grid)(th_ind + 2, x_ind + 2)}; + + double final_vals[4]; + final_vals[0] = fourth_order_lagange(x, vals1, p_x); + final_vals[1] = fourth_order_lagange(x, vals2, p_x); + final_vals[2] = fourth_order_lagange(x, vals3, p_x); + final_vals[3] = fourth_order_lagange(x, vals4, p_x); + + double f = fourth_order_lagange(y, final_vals, p_theta); + + return f; +} + +//////////////// BicubicSplineInterpolator implementation (production-grade, C2 continuity) ////////////////// + +double BicubicSplineInterpolator::interpolate(double p_x, double p_theta, + bool allow_extrapolation) const +{ + // Find cell containing the point + auto x_it = std::lower_bound(m_x.begin(), m_x.end(), p_x); + auto th_it = std::lower_bound(m_theta.begin(), m_theta.end(), p_theta); + + size_t i = std::max(1, std::distance(m_x.begin(), x_it)) - 1; + size_t j = std::max(1, std::distance(m_theta.begin(), th_it)) - 1; + + i = std::min(i, m_x.size() - 2); + j = std::min(j, m_theta.size() - 2); + + // Normalized coordinates within cell + double tx = (p_x - m_x[i]) / (m_x[i + 1] - m_x[i]); + double ty = (p_theta - m_theta[j]) / (m_theta[j + 1] - m_theta[j]); + + // Evaluate bicubic polynomial using precomputed coefficients + double result = 0.0; + for (int px = 0; px < 4; ++px) + { + for (int py = 0; py < 4; ++py) + { + double coeff = m_coeffs[j * (m_x.size() - 1) + i][py * 4 + px]; + result += coeff * std::pow(tx, px) * std::pow(ty, py); + } + } + + return result; +} + +void BicubicSplineInterpolator::compute_spline_coefficients(const Grid *grid) +{ + // Reserve space for all cells + size_t n_cells = (m_x.size() - 1) * (m_theta.size() - 1); + m_coeffs.resize(n_cells); + + // For each cell, compute the 16 bicubic coefficients + // This requires function values, derivatives, and cross-derivatives + // (Simplified version - full implementation would compute proper derivatives) + + for (size_t j = 0; j < m_theta.size() - 1; ++j) + { + for (size_t i = 0; i < m_x.size() - 1; ++i) + { + size_t cell_idx = j * (m_x.size() - 1) + i; + + // Get the 4 corner values + double f00 = (*grid)(j, i); + double f10 = (*grid)(j, i + 1); + double f01 = (*grid)(j + 1, i); + double f11 = (*grid)(j + 1, i + 1); + + // Estimate derivatives using finite differences + // IMPORTANT: Need to scale derivatives by cell width for bicubic Hermite formula + // which expects df/dt (derivative wrt normalized parameter), not df/dx + double fx00 = estimate_dx(grid, j, i) * (m_x[i + 1] - m_x[i]); + double fx10 = estimate_dx(grid, j, i + 1) * (m_x[i + 1] - m_x[i]); + double fx01 = estimate_dx(grid, j + 1, i) * (m_x[i + 1] - m_x[i]); + double fx11 = estimate_dx(grid, j + 1, i + 1) * (m_x[i + 1] - m_x[i]); + + double fy00 = estimate_dy(grid, j, i) * (m_theta[j + 1] - m_theta[j]); + double fy10 = estimate_dy(grid, j, i + 1) * (m_theta[j + 1] - m_theta[j]); + double fy01 = estimate_dy(grid, j + 1, i) * (m_theta[j + 1] - m_theta[j]); + double fy11 = estimate_dy(grid, j + 1, i + 1) * (m_theta[j + 1] - m_theta[j]); + + double fxy00 = estimate_dxdy(grid, j, i) * (m_x[i + 1] - m_x[i]) * (m_theta[j + 1] - m_theta[j]); + double fxy10 = estimate_dxdy(grid, j, i + 1) * (m_x[i + 1] - m_x[i]) * (m_theta[j + 1] - m_theta[j]); + double fxy01 = estimate_dxdy(grid, j + 1, i) * (m_x[i + 1] - m_x[i]) * (m_theta[j + 1] - m_theta[j]); + double fxy11 = estimate_dxdy(grid, j + 1, i + 1) * (m_x[i + 1] - m_x[i]) * (m_theta[j + 1] - m_theta[j]); + + // Solve for bicubic coefficients (matrix inversion) + compute_cell_coefficients(m_coeffs[cell_idx], + f00, f10, f01, f11, + fx00, fx10, fx01, fx11, + fy00, fy10, fy01, fy11, + fxy00, fxy10, fxy01, fxy11); + } + } +} + +double BicubicSplineInterpolator::estimate_dx(const Grid *grid, size_t j, size_t i) const +{ + if (i == 0) + return ((*grid)(j, i + 1) - (*grid)(j, i)) / (m_x[i + 1] - m_x[i]); + if (i == m_x.size() - 1) + return ((*grid)(j, i) - (*grid)(j, i - 1)) / (m_x[i] - m_x[i - 1]); + // Central difference + return ((*grid)(j, i + 1) - (*grid)(j, i - 1)) / (m_x[i + 1] - m_x[i - 1]); +} + +double BicubicSplineInterpolator::estimate_dy(const Grid *grid, size_t j, size_t i) const +{ + if (j == 0) + return ((*grid)(j + 1, i) - (*grid)(j, i)) / (m_theta[j + 1] - m_theta[j]); + if (j == m_theta.size() - 1) + return ((*grid)(j, i) - (*grid)(j - 1, i)) / (m_theta[j] - m_theta[j - 1]); + return ((*grid)(j + 1, i) - (*grid)(j - 1, i)) / (m_theta[j + 1] - m_theta[j - 1]); +} + +double BicubicSplineInterpolator::estimate_dxdy(const Grid *grid, size_t j, size_t i) const +{ + // Simplified cross-derivative estimation + if (i == 0 || i == m_x.size() - 1 || j == 0 || j == m_theta.size() - 1) + return 0.0; + + double dx = m_x[i + 1] - m_x[i - 1]; + double dy = m_theta[j + 1] - m_theta[j - 1]; + return ((*grid)(j + 1, i + 1) - (*grid)(j + 1, i - 1) - + (*grid)(j - 1, i + 1) + (*grid)(j - 1, i - 1)) / + (4.0 * dx * dy); +} + +void BicubicSplineInterpolator::compute_cell_coefficients(std::array &coeffs, + double f00, double f10, double f01, double f11, + double fx00, double fx10, double fx01, double fx11, + double fy00, double fy10, double fy01, double fy11, + double fxy00, double fxy10, double fxy01, double fxy11) const +{ + // Bicubic interpolation coefficient matrix + // This is a standard formula - see Numerical Recipes + static const int A[16][16] = { + {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {-3, 3, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {2, -2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, -2, -1, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 1, 1, 0, 0}, + {-3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, -2, 0, -1, 0}, + {9, -9, -9, 9, 6, 3, -6, -3, 6, -6, 3, -3, 4, 2, 2, 1}, + {-6, 6, 6, -6, -3, -3, 3, 3, -4, 4, -2, 2, -2, -2, -1, -1}, + {2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0}, + {-6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1}, + {4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1}}; + + double x[16] = {f00, f10, f01, f11, + fx00, fx10, fx01, fx11, + fy00, fy10, fy01, fy11, + fxy00, fxy10, fxy01, fxy11}; + + for (int i = 0; i < 16; ++i) + { + coeffs[i] = 0.0; + for (int j = 0; j < 16; ++j) + { + coeffs[i] += A[i][j] * x[j]; + } + } +} diff --git a/FOORT/src/Interpolator.h b/FOORT/src/Interpolator.h new file mode 100644 index 0000000..de1d08c --- /dev/null +++ b/FOORT/src/Interpolator.h @@ -0,0 +1,106 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Grid.h" + +#ifndef INTERPOLATOR_H +#define INTERPOLATOR_H + +class Interpolator +{ +public: + double *m_x; // Pointer to the x values + double *m_theta; // Pointer to the theta values + + //! Constructor + Interpolator(int dim_x, int dim_th, std::string x_file, std::string th_file); + + //! Destructor + ~Interpolator() + { + delete[] this->m_x; + delete[] this->m_theta; + } + + //! Interpolates the data of the grid + //! @param grid The grid to be interpolated + double interpolate(Grid *grid, double p_x, double p_theta); + +private: + //! Returns the fourth order lagrangian polynomial evaluated at the point of interest + double fourth_order_lagange(double x_points[4], double y_points[4], double x); +}; + +class BicubicSplineInterpolator +{ +public: + BicubicSplineInterpolator(const std::string x_file = "RotatingBosonStar/x.txt", + const std::string th_file = "RotatingBosonStar/theta.txt") + { + // Load x and theta coordinates from files + std::ifstream x_input(x_file); + std::ifstream th_input(th_file); + + if (!th_input || !x_input) + { + std::cerr << "Error opening file!" << std::endl; + throw std::runtime_error("File not found"); + } + + double val; + while (x_input >> val) + { + m_x.push_back(val); + } + while (th_input >> val) + { + m_theta.push_back(val); + } + + // Explicitly close files + x_input.close(); + th_input.close(); + } + + // BicubicSplineInterpolator(const std::vector &x_coords, + // const std::vector &theta_coords, + // const Grid *grid) + // : m_x(x_coords), m_theta(theta_coords) + // { + + // compute_spline_coefficients(grid); + // } + + void set_grid(const Grid *grid) + { + // Precompute spline coefficients for better performance + // This is more complex but provides C2 continuity + compute_spline_coefficients(grid); + std::cout << "Interpolation domain set: x in [" << m_x.front() << ", " << m_x.back() << "], theta in [" << m_theta.front() << ", " << m_theta.back() << "].\n"; + } + double interpolate(double p_x, double p_theta, + bool allow_extrapolation = false) const; + +private: + std::vector m_x; + std::vector m_theta; + std::vector> m_coeffs; // 16 coefficients per cell + + void compute_spline_coefficients(const Grid *grid); + double estimate_dx(const Grid *grid, size_t j, size_t i) const; + double estimate_dy(const Grid *grid, size_t j, size_t i) const; + double estimate_dxdy(const Grid *grid, size_t j, size_t i) const; + void compute_cell_coefficients(std::array &coeffs, + double f00, double f10, double f01, double f11, + double fx00, double fx10, double fx01, double fx11, + double fy00, double fy10, double fy01, double fy11, + double fxy00, double fxy10, double fxy01, double fxy11) const; +}; + +#endif // INTERPOLATOR_H diff --git a/FOORT/src/Metric.cpp b/FOORT/src/Metric.cpp index 9b807b8..f249d5e 100644 --- a/FOORT/src/Metric.cpp +++ b/FOORT/src/Metric.cpp @@ -5,9 +5,12 @@ #include // needed for sqrt() and sin() etc (only on Linux) #include // needed for std::find +#include // needed for std::setprecision, for error messages -#include "Spline.h" // needed for spline interpolation -#include // needed for string stream +#include "Spline.h" // needed for spline interpolation +#include "Grid.h" // needed for the grid class, for the Rotating Boson star metric +#include "Interpolator.h" // needed for the Interpolator class +#include // needed for string stream /** * @file Metric.h @@ -1223,4 +1226,150 @@ std::string BosonStarMetric::getFullDescriptionStr() const return "Boson star (Phi infinity = " + std::to_string(m_Phi_infinity) + ", num lines = " + std::to_string(m_num_lines) + ")"; } +// RotatingBosonStarMetric functions (implementation by Seppe Staelens) + +/** + * @brief Construct a new Rotating Boson Star Metric object + * @param rLogScale whether we are using a logarithmic radial scale + * @param FlipAngularMomentum whether to flip the sign of Omega (reverse rotation direction) + */ +RotatingBosonStarMetric::RotatingBosonStarMetric(bool rLogScale, std::string MetricFolder, + int num_x, int num_th, real L, bool FlipAngularMomentum) : Metric(rLogScale), + m_grid_f(new Grid(num_th, num_x)), + m_grid_l(new Grid(num_th, num_x)), + m_grid_g(new Grid(num_th, num_x)), + m_grid_Omega(new Grid(num_th, num_x)), + m_fInterpolator(new BicubicSplineInterpolator(MetricFolder + "x.txt", + MetricFolder + "theta.txt")), + m_lInterpolator(new BicubicSplineInterpolator(MetricFolder + "x.txt", + MetricFolder + "theta.txt")), + m_gInterpolator(new BicubicSplineInterpolator(MetricFolder + "x.txt", + MetricFolder + "theta.txt")), + m_OmegaInterpolator(new BicubicSplineInterpolator(MetricFolder + "x.txt", + MetricFolder + "theta.txt")), + m_L(L), + m_OmegaSign(FlipAngularMomentum ? -1 : 1) + +{ + // Make sure we are in four spacetime dimensions + if constexpr (dimension != 4) + { + ScreenOutput("Rotating boson star is only defined in four dimensions!", OutputLevel::Level_0_WARNING); + } + // Rotating Boson star has a Killing vector along t and phi, so we initialize the symmetries accordingly + m_Symmetries = {0, 3}; + + // Read the different metric functions + m_grid_f->initialize_from_file(MetricFolder + "f.txt"); + m_fInterpolator->set_grid(m_grid_f); + + m_grid_l->initialize_from_file(MetricFolder + "l.txt"); + m_lInterpolator->set_grid(m_grid_l); + + m_grid_g->initialize_from_file(MetricFolder + "g.txt"); + m_gInterpolator->set_grid(m_grid_g); + + m_grid_Omega->initialize_from_file(MetricFolder + "omega.txt"); + m_OmegaInterpolator->set_grid(m_grid_Omega); +} + +RotatingBosonStarMetric::~RotatingBosonStarMetric() +{ + delete m_grid_f; + delete m_grid_l; + delete m_grid_g; + delete m_grid_Omega; + delete m_fInterpolator; + delete m_lInterpolator; + delete m_gInterpolator; + delete m_OmegaInterpolator; +} + +/** + * @brief RotatingBosonStar metric getter, indices down + * @param p Point at which to evaluate the metric + * @return TwoIndex + */ +TwoIndex RotatingBosonStarMetric::getMetric_dd(const Point &p) const +{ + // spherical coordinates + // If logscale is turned on, then the first coordinate is actually u = log(r), so r = e^u + real r = m_rLogScale ? exp(p[1]) : p[1]; + r += 1e-9; // to avoid r = 0 + real x = m_L * r / (1. + r); + real theta = p[2]; + real sint = sin(theta); + + // real f = m_interpolator->interpolate(m_grid_f, x, theta); + // real l = m_interpolator->interpolate(m_grid_l, x, theta); + // real g = m_interpolator->interpolate(m_grid_g, x, theta); + // real Omega = m_interpolator->interpolate(m_grid_Omega, x, theta); + + real f = m_fInterpolator->interpolate(x, theta); + real l = m_lInterpolator->interpolate(x, theta); + real g = m_gInterpolator->interpolate(x, theta); + real Omega = m_OmegaSign * m_OmegaInterpolator->interpolate(x, theta); + + // Covariant metric elements + real g00 = -(f - l * Omega * Omega * sint * sint / f); + real g11 = l * g / f; + real g22 = g11 * r * r; + real g33 = l * r * r * sint * sint / f; + real g03 = -l * r * Omega * sint * sint / f; + // If the log scale is set on, the true coordinate we are calculating the metric in is u = log(r), so dr = r du + if (m_rLogScale) + { + g11 *= (r * r); + } + return TwoIndex{{{g00, 0, 0, g03}, {0, g11, 0, 0}, {0, 0, g22, 0}, {g03, 0, 0, g33}}}; +} + +/** + * @brief RotatingBosonStar metric getter, indices up + * @param p Point at which to evaluate the metric + * @return TwoIndex + */ +TwoIndex RotatingBosonStarMetric::getMetric_uu(const Point &p) const +{ + // spherical coordinates + // If logscale is turned on, then the first coordinate is actually u = log(r), so r = e^u + real r = m_rLogScale ? exp(p[1]) : p[1]; + r += 1e-9; // to avoid r = 0 + real x = m_L * r / (1. + r); + real theta = p[2]; + real sint = sin(theta); + + // real f = m_interpolator->interpolate(m_grid_f, x, theta); + // real l = m_interpolator->interpolate(m_grid_l, x, theta); + // real g = m_interpolator->interpolate(m_grid_g, x, theta); + // real Omega = m_interpolator->interpolate(m_grid_Omega, x, theta); + real f = m_fInterpolator->interpolate(x, theta); + real l = m_lInterpolator->interpolate(x, theta); + real g = m_gInterpolator->interpolate(x, theta); + real Omega = m_OmegaSign * m_OmegaInterpolator->interpolate(x, theta); + + // Contravariant metric elements + real g00 = -1 / f; + real g11 = f / (l * g); + real g22 = g11 / (r * r); + real g33 = (f - l * Omega * Omega * sint * sint / f) / (l * r * r * sint * sint); + real g03 = -Omega / (r * f); + // If the log scale is set on, the true coordinate we are calculating the metric in is u = log(r), so , so dr = r du + if (m_rLogScale) + { + g11 *= 1.0 / (r * r); + } + + return TwoIndex{{{g00, 0, 0, g03}, {0, g11, 0, 0}, {0, 0, g22, 0}, {g03, 0, 0, g33}}}; +} + +/** + * @brief RotatingBosonStar metric description string getter + * @return std::string + */ +std::string RotatingBosonStarMetric::getFullDescriptionStr() const +{ + return std::string("Rotating boson star") + (m_OmegaSign == -1 ? " (flipped angular momentum)" : ""); +} + //// (New Metric classes can define their member functions here) diff --git a/FOORT/src/Metric.h b/FOORT/src/Metric.h index 7918252..d27e848 100644 --- a/FOORT/src/Metric.h +++ b/FOORT/src/Metric.h @@ -3,9 +3,11 @@ #include "Geometry.h" // Needed for basic tensor objects etc. -#include "Spline.h" -#include // for strings -#include // needed for the (non-fixed size) vector of symmetries in the metric +#include "Spline.h" // needed for spline interpolation, for the Boson star metric +#include "Grid.h" // needed for the grid class, for the Rotating Boson star metric +#include "Interpolator.h" // needed for the Interpolator class, for the Rotating Boson star metric +#include // for strings +#include // needed for the (non-fixed size) vector of symmetries in the metric /** * @file Metric.h @@ -330,6 +332,41 @@ class BosonStarMetric final : public Metric void read_data(); }; +class RotatingBosonStarMetric final : public Metric +{ +public: + // Simple (default) constructor is all that is needed + RotatingBosonStarMetric(bool rLogScale = false, std::string MetricFolder = "RotatingBosonStar/data_Will/", int num_x = 500, int num_th = 399, real L = 1., bool FlipAngularMomentum = false); + ~RotatingBosonStarMetric(); + + // The override of the basic metric getter functions + TwoIndex getMetric_dd(const Point &p) const final; + TwoIndex getMetric_uu(const Point &p) const final; + // The override of the description string getter + std::string getFullDescriptionStr() const final; + +private: + //! The grids with the metric functions + Grid *m_grid_f; + Grid *m_grid_l; + Grid *m_grid_g; + Grid *m_grid_Omega; + +public: + //! The grid interpolators + BicubicSplineInterpolator *m_fInterpolator; + BicubicSplineInterpolator *m_lInterpolator; + BicubicSplineInterpolator *m_gInterpolator; + BicubicSplineInterpolator *m_OmegaInterpolator; + + const real m_L; + // Sign multiplier for Omega: +1 for normal rotation, -1 for flipped rotation + const int m_OmegaSign; + + //! The interpolator for the metric functions + // Interpolator *m_interpolator; +}; + //// METRIC ADD POINT A //// // Declare your new Metric class here, publically inheriting from the base class Metric // (or SphericalHorizonMetric if your Metric has a horizon, or SingularityMetric if your Metric has other, arbitrary singularities) diff --git a/FOORT/src/ViewScreen.h b/FOORT/src/ViewScreen.h index c051ff6..8296a2a 100644 --- a/FOORT/src/ViewScreen.h +++ b/FOORT/src/ViewScreen.h @@ -67,6 +67,10 @@ class ViewScreen OutputLevel::Level_0_WARNING); } + std::cout << "ViewScreen initialized with position " << toString(m_Pos) << ", looking direction " << toString(m_Direction) + << ", screen size " << toString(m_ScreenSize) << ", screen center " << toString(m_ScreenCenter) + << ", geodesic type " << (m_GeodType == GeodesicType::Null ? "Null" : (m_GeodType == GeodesicType::Timelike ? "Timelike" : "Spacelike")) + << ", and metric: " << m_theMetric->getFullDescriptionStr() << "\n"; // Construct the vielbein now ConstructVielbein(); } diff --git a/FOORT/test/CMakeLists.txt b/FOORT/test/CMakeLists.txt index 4c5cf64..9c6812e 100644 --- a/FOORT/test/CMakeLists.txt +++ b/FOORT/test/CMakeLists.txt @@ -1,12 +1,24 @@ add_executable(test_integrator test_integrator.cpp) add_executable(test_einstein_ring test_einstein_ring.cpp) +add_executable(test_rotating_boson_star_geodesic test_rotating_boson_star_geodesic.cpp) +add_executable(test_boson_star_geodesic test_boson_star_geodesic.cpp) +add_executable(test_boson_star_interpolators test_boson_star_interpolators.cpp) target_include_directories(test_integrator PUBLIC ${CMAKE_SOURCE_DIR}/FOORT/src) target_include_directories(test_einstein_ring PUBLIC ${CMAKE_SOURCE_DIR}/FOORT/src) +target_include_directories(test_rotating_boson_star_geodesic PUBLIC ${CMAKE_SOURCE_DIR}/FOORT/src) +target_include_directories(test_boson_star_geodesic PUBLIC ${CMAKE_SOURCE_DIR}/FOORT/src) +target_include_directories(test_boson_star_interpolators PUBLIC ${CMAKE_SOURCE_DIR}/FOORT/src) target_link_libraries(test_integrator PUBLIC config_lib metric_lib input_output_lib GTest::gtest_main) target_link_libraries(test_einstein_ring PUBLIC config_lib metric_lib input_output_lib GTest::gtest_main) +target_link_libraries(test_rotating_boson_star_geodesic PUBLIC config_lib metric_lib input_output_lib GTest::gtest_main) +target_link_libraries(test_boson_star_geodesic PUBLIC config_lib metric_lib input_output_lib GTest::gtest_main) +target_link_libraries(test_boson_star_interpolators PUBLIC config_lib metric_lib input_output_lib GTest::gtest_main) include(GoogleTest) gtest_discover_tests(test_integrator) gtest_discover_tests(test_einstein_ring) +gtest_discover_tests(test_rotating_boson_star_geodesic) +gtest_discover_tests(test_boson_star_geodesic) +gtest_discover_tests(test_boson_star_interpolators) diff --git a/FOORT/test/test_boson_star_geodesic.cpp b/FOORT/test/test_boson_star_geodesic.cpp new file mode 100644 index 0000000..8f7a51c --- /dev/null +++ b/FOORT/test/test_boson_star_geodesic.cpp @@ -0,0 +1,213 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "Config.h" +#include "Diagnostics.h" +#include "Geodesic.h" +#include "Geometry.h" +#include "InputOutput.h" +#include "Integrators.h" +#include "Metric.h" +#include "Terminations.h" + +namespace +{ + constexpr double kPhiInfinity = 1.376427; + constexpr int kNumLines = 10896; + constexpr real kStableLightRingRadius = 28.1; + constexpr real kUnstableLightRingRadius = 31.202133; + + struct GeodesicResult + { + int step_count; + real final_r; + real closest_r; + }; + + void EnsureBosonStarDataVisible() + { + namespace fs = std::filesystem; + + const fs::path phi_path = fs::path("BosonStar") / "Phi.dat"; + const fs::path m_path = fs::path("BosonStar") / "m.dat"; + + if (fs::exists(phi_path) && fs::exists(m_path)) + { + return; + } + + fs::path probe = fs::current_path(); + for (int i = 0; i < 8; ++i) + { + if (fs::exists(probe / phi_path) && fs::exists(probe / m_path)) + { + fs::current_path(probe); + return; + } + + const fs::path foort_root = probe / "FOORT"; + if (fs::exists(foort_root / phi_path) && fs::exists(foort_root / m_path)) + { + fs::current_path(foort_root); + return; + } + + if (!probe.has_parent_path()) + { + break; + } + probe = probe.parent_path(); + } + + throw std::runtime_error("Could not locate BosonStar/Phi.dat and BosonStar/m.dat."); + } + + OneIndex RaiseIndex(const TwoIndex &guu, const OneIndex &covector) + { + OneIndex vector = {0.0, 0.0, 0.0, 0.0}; + for (int i = 0; i < dimension; ++i) + { + for (int j = 0; j < dimension; ++j) + { + vector[i] += guu[i][j] * covector[j]; + } + } + return vector; + } + + real NullAngularMomentumAtRadius(const Metric *metric, real radius, real energy) + { + Point orbit_pos = {0.0, radius, pi / 2.0, 0.0}; + TwoIndex guu = metric->getMetric_uu(orbit_pos); + + const real gu_tt = guu[0][0]; + const real gu_tphi = guu[0][3]; + const real gu_phiphi = guu[3][3]; + + const real discriminant = gu_tphi * gu_tphi * energy * energy - gu_tt * gu_phiphi * energy * energy; + if (!std::isfinite(discriminant) || discriminant <= 0.0 || gu_phiphi == 0.0) + { + throw std::runtime_error("Failed to compute null angular momentum for BosonStarMetric."); + } + + return (gu_tphi * energy + std::sqrt(discriminant)) / gu_phiphi; + } + + GeodesicResult IntegrateBosonStarLightRing(real angular_momentum_radius, + real start_radius, + const std::string &output_filename) + { + EnsureBosonStarDataVisible(); + + const bool rLogScale = false; + std::unique_ptr theM = std::unique_ptr( + new BosonStarMetric(kPhiInfinity, kNumLines, rLogScale)); + std::unique_ptr theS = std::unique_ptr(new NoSource(theM.get())); + + GeodesicPositionDiagnostic::DiagOptions = + std::unique_ptr( + new GeodesicPositionOptions{0, UpdateFrequency{1, false, false}}); + + ClosestRadiusDiagnostic::DiagOptions = + std::unique_ptr( + new ClosestRadiusOptions{rLogScale, UpdateFrequency{1, false, false}}); + + const DiagBitflag all_diags = Diag_GeodesicPosition | Diag_ClosestRadius; + const DiagBitflag val_diag = Diag_ClosestRadius; + + const TermBitflag all_terms = Term_BoundarySphere | Term_TimeOut | Term_NaN; + + BoundarySphereTermination::TermOptions = + std::unique_ptr( + new BoundarySphereTermOptions{120.0, rLogScale, 1}); + TimeOutTermination::TermOptions = + std::unique_ptr( + new TimeOutTermOptions{3000, 1}); + NaNTermination::TermOptions = + std::unique_ptr( + new NaNTermOptions{true, 1}); + + GeodesicIntegratorFunc theIntegrator = Integrators::IntegrateGeodesicStep_RK4; + Integrators::IntegratorDescription = "RK4"; + Integrators::epsilon = 0.03; + + Geodesic theGeod(theM.get(), theS.get(), all_diags, val_diag, all_terms, theIntegrator); + + const real energy = 1.0; + const real angular_momentum = NullAngularMomentumAtRadius(theM.get(), angular_momentum_radius, energy); + const OneIndex initvel_d = {-energy, 0.0, 0.0, angular_momentum}; + + Point initpos = {0.0, start_radius, pi / 2.0, 0.0}; + TwoIndex guu_at_start = theM->getMetric_uu(initpos); + OneIndex initvel = RaiseIndex(guu_at_start, initvel_d); + + ScreenIndex scrindex = {0, 0}; + theGeod.Reset(scrindex, initpos, initvel); + + std::ofstream output(output_filename); + output << std::setprecision(15); + output << "# t r theta phi\n"; + + Point current_pos = initpos; + int step_count = 0; + + while (theGeod.getTermCondition() == Term::Continue) + { + theGeod.Update(); + current_pos = theGeod.getCurrentPos(); + + const real current_r = rLogScale ? exp(current_pos[1]) : current_pos[1]; + output << current_pos[0] << " " << current_r << " " + << current_pos[2] << " " << current_pos[3] << "\n"; + ++step_count; + } + + output.close(); + + const real final_r = rLogScale ? exp(current_pos[1]) : current_pos[1]; + const real closest_r = rLogScale ? exp(theGeod.getDiagnosticFinalValue()[0]) : theGeod.getDiagnosticFinalValue()[0]; + + return GeodesicResult{step_count, final_r, closest_r}; + } +} // namespace + +TEST(BosonStar, stable_light_ring) +{ + const GeodesicResult result = IntegrateBosonStarLightRing( + kStableLightRingRadius, + kStableLightRingRadius, + "boson_star_stable_light_ring.dat"); + + EXPECT_GT(result.step_count, 0) << "Geodesic should integrate for at least one step."; + EXPECT_TRUE(std::isfinite(result.final_r)) << "Final radius should be finite."; + EXPECT_TRUE(std::isfinite(result.closest_r)) << "Closest radius diagnostic should be finite."; + + EXPECT_NEAR(result.closest_r, kStableLightRingRadius, 0.3) + << "Stable light ring closest radius should stay near 28.1."; + EXPECT_NEAR(result.final_r, kStableLightRingRadius, 0.6) + << "Stable light ring final radius should remain near 28.1."; +} + +TEST(BosonStar, unstable_light_ring) +{ + constexpr real perturbation = 0.0000001; + const GeodesicResult result = IntegrateBosonStarLightRing( + kUnstableLightRingRadius, + kUnstableLightRingRadius + perturbation, + "boson_star_unstable_light_ring.dat"); + + EXPECT_GT(result.step_count, 0) << "Geodesic should integrate for at least one step."; + EXPECT_TRUE(std::isfinite(result.final_r)) << "Final radius should be finite."; + + EXPECT_GT(std::abs(result.final_r - kUnstableLightRingRadius), 0.05) + << "Unstable light ring trajectory should depart from areal radius 31.1."; + EXPECT_GT(std::abs(result.final_r - (kUnstableLightRingRadius + perturbation)), 0.01) + << "Trajectory should move away from the perturbed starting radius."; +} diff --git a/FOORT/test/test_boson_star_interpolators.cpp b/FOORT/test/test_boson_star_interpolators.cpp new file mode 100644 index 0000000..ce68497 --- /dev/null +++ b/FOORT/test/test_boson_star_interpolators.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include + +#include "Metric.h" // Metrics +#include "Geometry.h" + +/** + * @brief Test interpolators in RotatingBosonStar metric + * @details Tests f, l, g, and omega interpolators at various x and theta values + */ +TEST(RotatingBosonStarInterpolators, test_interpolators_at_various_points) +{ + // Setup metric with same parameters as main test + bool rLogScale = true; + std::string MetricFolder = "RotatingBosonStar/data_C38/"; + int NumX = 1500; + int NumTh = 200; + real L = 1.0; + + std::unique_ptr theM = std::unique_ptr( + new RotatingBosonStarMetric(rLogScale, MetricFolder, NumX, NumTh, L)); + + auto metric = static_cast(theM.get()); + + // Open output file to log results + std::ofstream outfile("interpolator_test_results.txt"); + outfile << std::scientific << std::setprecision(10); + outfile << "Testing RotatingBosonStarMetric interpolators\n"; + outfile << "MetricFolder: " << MetricFolder << "\n"; + outfile << "NumX: " << NumX << ", NumTh: " << NumTh << ", L: " << L << "\n"; + outfile << "rLogScale: " << (rLogScale ? "true" : "false") << "\n\n"; + + // Test points: various x and theta values + std::vector x_values = {0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.99}; + std::vector theta_values = {0.01, 0.1, M_PI / 4.0, M_PI / 2.0, 3.0 * M_PI / 4.0, M_PI - 0.1, M_PI - 0.01}; + + outfile << "x values to test: "; + for (auto x : x_values) + outfile << x << " "; + outfile << "\n"; + outfile << "theta values to test: "; + for (auto th : theta_values) + outfile << th << " "; + outfile << "\n\n"; + + outfile << std::left << std::setw(12) << "x" + << std::setw(12) << "theta" + << std::setw(15) << "f" + << std::setw(15) << "l" + << std::setw(15) << "g" + << std::setw(15) << "omega\n"; + outfile << std::string(69, '-') << "\n"; + + int nan_count = 0; + int total_count = 0; + + for (auto x : x_values) + { + for (auto theta : theta_values) + { + try + { + // Get interpolated values + real f_val = metric->m_fInterpolator->interpolate(x, theta); + real l_val = metric->m_lInterpolator->interpolate(x, theta); + real g_val = metric->m_gInterpolator->interpolate(x, theta); + real omega_val = metric->m_OmegaInterpolator->interpolate(x, theta); + + total_count++; + + // Check for NaN + bool has_nan = std::isnan(f_val) || std::isnan(l_val) || + std::isnan(g_val) || std::isnan(omega_val); + + if (has_nan) + nan_count++; + + // Output results + outfile << std::setw(12) << x + << std::setw(12) << theta + << std::setw(15) << f_val + << std::setw(15) << l_val + << std::setw(15) << g_val + << std::setw(15) << omega_val; + + if (has_nan) + outfile << " <- NAN"; + outfile << "\n"; + } + catch (const std::exception &e) + { + total_count++; + outfile << std::setw(12) << x + << std::setw(12) << theta + << "ERROR: " << e.what() << "\n"; + } + } + } + + outfile << "\n" + << std::string(69, '=') << "\n"; + outfile << "Summary:\n"; + outfile << "Total test points: " << total_count << "\n"; + outfile << "Points with NaN: " << nan_count << "\n"; + outfile << "Success rate: " << (total_count - nan_count) << "/" << total_count << "\n"; + + outfile.close(); + + // Special test for specific geodesic coordinates + double test_x = 0.998752; + double test_theta = 1.57; + + real f_test = metric->m_fInterpolator->interpolate(test_x, test_theta); + real l_test = metric->m_lInterpolator->interpolate(test_x, test_theta); + real g_test = metric->m_gInterpolator->interpolate(test_x, test_theta); + real omega_test = metric->m_OmegaInterpolator->interpolate(test_x, test_theta); + + // Basic assertions + ASSERT_GT(total_count, 0) << "Should have tested at least one point"; + ASSERT_EQ(nan_count, 0) << "Should not have NaN values in interpolations"; + ASSERT_TRUE(std::isfinite(f_test)); + ASSERT_TRUE(std::isfinite(l_test)); + ASSERT_TRUE(std::isfinite(g_test)); + ASSERT_TRUE(std::isfinite(omega_test)); +} diff --git a/FOORT/test/test_rotating_boson_star_geodesic.cpp b/FOORT/test/test_rotating_boson_star_geodesic.cpp new file mode 100644 index 0000000..84677ab --- /dev/null +++ b/FOORT/test/test_rotating_boson_star_geodesic.cpp @@ -0,0 +1,248 @@ +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "Config.h" +#include "Diagnostics.h" +#include "Geodesic.h" +#include "Geometry.h" +#include "InputOutput.h" +#include "Integrators.h" +#include "Metric.h" +#include "Terminations.h" + +namespace +{ + constexpr real kStableLightRingRadius = 15.47; + constexpr real kUnstableLightRingRadius = 56.507305; + + constexpr int kNumX = 1500; + constexpr int kNumTh = 200; + constexpr real kL = 1.0; + + struct GeodesicResult + { + int step_count; + real final_r; + real closest_r; + }; + + std::string ResolveRotatingBosonStarDataFolder() + { + namespace fs = std::filesystem; + + const fs::path relative_data_folder = fs::path("RotatingBosonStar") / "data_C38"; + if (fs::exists(relative_data_folder)) + { + return (relative_data_folder.string() + "/"); + } + + fs::path probe = fs::current_path(); + for (int i = 0; i < 8; ++i) + { + if (fs::exists(probe / relative_data_folder)) + { + fs::current_path(probe); + return (relative_data_folder.string() + "/"); + } + + const fs::path foort_root = probe / "FOORT"; + if (fs::exists(foort_root / relative_data_folder)) + { + fs::current_path(foort_root); + return (relative_data_folder.string() + "/"); + } + + if (!probe.has_parent_path()) + { + break; + } + probe = probe.parent_path(); + } + + throw std::runtime_error("Could not locate RotatingBosonStar/data_C38."); + } + + OneIndex RaiseIndex(const TwoIndex &guu, const OneIndex &covector) + { + OneIndex vector = {0.0, 0.0, 0.0, 0.0}; + for (int i = 0; i < dimension; ++i) + { + for (int j = 0; j < dimension; ++j) + { + vector[i] += guu[i][j] * covector[j]; + } + } + return vector; + } + + void CheckMetricInvertibility(const TwoIndex &gdd, const TwoIndex &guu) + { + for (int a = 0; a < dimension; ++a) + { + for (int d = 0; d < dimension; ++d) + { + real sum = 0.0; + for (int c = 0; c < dimension; ++c) + { + sum += guu[a][c] * gdd[c][d]; + } + + if (a == d) + { + EXPECT_NEAR(sum, 1.0, 1e-13) + << "Metric invertibility check failed at a = " << a << ", d = " << d; + } + else + { + EXPECT_NEAR(sum, 0.0, 1e-13) + << "Metric invertibility check failed at a = " << a << ", d = " << d; + } + } + } + } + + real NullAngularMomentumAtRadius(const Metric *metric, real radius, real energy) + { + Point orbit_pos = {0.0, radius, pi / 2.0, 0.0}; + const TwoIndex guu = metric->getMetric_uu(orbit_pos); + + const real gu_tt = guu[0][0]; + const real gu_tphi = guu[0][3]; + const real gu_phiphi = guu[3][3]; + + const real discriminant = gu_tphi * gu_tphi * energy * energy - gu_tt * gu_phiphi * energy * energy; + if (!std::isfinite(discriminant) || discriminant <= 0.0 || gu_phiphi == 0.0) + { + throw std::runtime_error("Failed to compute null angular momentum for RotatingBosonStarMetric."); + } + + return (gu_tphi * energy - std::sqrt(discriminant)) / gu_phiphi; + } + + GeodesicResult IntegrateRotatingBosonStarLightRing(real angular_momentum_radius, + real start_radius, + const std::string &output_filename) + { + const bool rLogScale = false; + const std::string metric_folder = ResolveRotatingBosonStarDataFolder(); + + std::unique_ptr theM = std::unique_ptr( + new RotatingBosonStarMetric(rLogScale, metric_folder, kNumX, + kNumTh, kL, true)); + std::unique_ptr theS = std::unique_ptr(new NoSource(theM.get())); + + GeodesicPositionDiagnostic::DiagOptions = + std::unique_ptr( + new GeodesicPositionOptions{0, UpdateFrequency{1, false, false}}); + + ClosestRadiusDiagnostic::DiagOptions = + std::unique_ptr( + new ClosestRadiusOptions{rLogScale, UpdateFrequency{1, false, false}}); + + const DiagBitflag all_diags = Diag_GeodesicPosition | Diag_ClosestRadius; + const DiagBitflag val_diag = Diag_ClosestRadius; + + const TermBitflag all_terms = Term_BoundarySphere | Term_TimeOut | Term_NaN; + BoundarySphereTermination::TermOptions = + std::unique_ptr( + new BoundarySphereTermOptions{120.0, rLogScale, 1}); + TimeOutTermination::TermOptions = + std::unique_ptr( + new TimeOutTermOptions{50000, 1}); + NaNTermination::TermOptions = + std::unique_ptr( + new NaNTermOptions{true, 1}); + + GeodesicIntegratorFunc theIntegrator = Integrators::IntegrateGeodesicStep_RK4; + Integrators::IntegratorDescription = "RK4"; + Integrators::epsilon = 0.01; + + Geodesic theGeod(theM.get(), theS.get(), all_diags, val_diag, all_terms, theIntegrator); + + const real energy = 1.0; + const real angular_momentum = NullAngularMomentumAtRadius(theM.get(), angular_momentum_radius, energy); + const OneIndex initvel_d = {-energy, 0.0, 0.0, angular_momentum}; + std::cout << "Integrating geodesic with angular momentum radius " << angular_momentum_radius + << " and corresponding null angular momentum " << angular_momentum << ".\n"; + + const Point initpos = {0.0, rLogScale ? log(start_radius) : start_radius, pi / 2.0, 0.0}; + const TwoIndex gdd_at_start = theM->getMetric_dd(initpos); + const TwoIndex guu_at_start = theM->getMetric_uu(initpos); + CheckMetricInvertibility(gdd_at_start, guu_at_start); + const OneIndex initvel = RaiseIndex(guu_at_start, initvel_d); + + ScreenIndex scrindex = {0, 0}; + theGeod.Reset(scrindex, initpos, initvel); + + std::ofstream output(output_filename); + output << std::setprecision(15); + output << "# t r theta phi\n"; + + Point current_pos = initpos; + int step_count = 0; + + while (theGeod.getTermCondition() == Term::Continue) + { + theGeod.Update(); + current_pos = theGeod.getCurrentPos(); + + const TwoIndex gdd = theM->getMetric_dd(current_pos); + const TwoIndex guu = theM->getMetric_uu(current_pos); + CheckMetricInvertibility(gdd, guu); + + const real current_r = rLogScale ? exp(current_pos[1]) : current_pos[1]; + output << current_pos[0] << " " << current_r << " " + << current_pos[2] << " " << current_pos[3] << "\n"; + ++step_count; + } + + output.close(); + + const real final_r = rLogScale ? exp(current_pos[1]) : current_pos[1]; + const real closest_r = rLogScale ? exp(theGeod.getDiagnosticFinalValue()[0]) : theGeod.getDiagnosticFinalValue()[0]; + + return GeodesicResult{step_count, final_r, closest_r}; + } +} // namespace + +TEST(RotatingBosonStar, stable_light_ring) +{ + const GeodesicResult result = IntegrateRotatingBosonStarLightRing( + kStableLightRingRadius, + kStableLightRingRadius, + "rotating_boson_star_stable_light_ring.dat"); + + EXPECT_GT(result.step_count, 0) << "Geodesic should integrate for at least one step."; + EXPECT_TRUE(std::isfinite(result.final_r)) << "Final radius should be finite."; + EXPECT_TRUE(std::isfinite(result.closest_r)) << "Closest radius diagnostic should be finite."; + + EXPECT_NEAR(result.closest_r, kStableLightRingRadius, 0.3) + << "Stable light ring closest radius should stay near 15.5."; + EXPECT_NEAR(result.final_r, kStableLightRingRadius, 0.6) + << "Stable light ring final radius should remain near 15.5."; +} + +TEST(RotatingBosonStar, unstable_light_ring) +{ + constexpr real perturbation = 1e-6; + const GeodesicResult result = IntegrateRotatingBosonStarLightRing( + kUnstableLightRingRadius, + kUnstableLightRingRadius + perturbation, + "rotating_boson_star_unstable_light_ring.dat"); + + EXPECT_GT(result.step_count, 0) << "Geodesic should integrate for at least one step."; + EXPECT_TRUE(std::isfinite(result.final_r)) << "Final radius should be finite."; + EXPECT_TRUE(std::isfinite(result.closest_r)) << "Closest radius diagnostic should be finite."; + + EXPECT_GT(std::abs(result.final_r - kUnstableLightRingRadius), 0.1) + << "Unstable light ring trajectory should depart from radius 56."; + EXPECT_GT(std::abs(result.final_r - (kUnstableLightRingRadius + perturbation)), 0.05) + << "Trajectory should move away from the perturbed starting radius."; +} diff --git a/PostProc/PhotonRingInterferometry.py b/PostProc/PhotonRingInterferometry.py index 250d58d..0598c2e 100644 --- a/PostProc/PhotonRingInterferometry.py +++ b/PostProc/PhotonRingInterferometry.py @@ -40,14 +40,22 @@ def RadonTransform( def RadonToComplexVis( - FOORTRadon: np.ndarray, PaddingFactor: float = 25, Verbose: bool = True -) -> np.ndarray: + FOORTRadon: np.ndarray, + PaddingFactor: float = 25, + Verbose: bool = True, + sample_spacing: float = 1.0, +) -> tuple[np.ndarray, np.ndarray]: """! @brief Calculates the complex visibility from a given radon transform. @param FOORTRadon: The radon transform to calculate the complex visibility from. @param PaddingFactor: The factor to pad the radon transform by (default = 25). @param Verbose: Whether to print out progress information (default = True). + @param sample_spacing: The spacing between samples in the radon transform (default = 1.0). @return complvis: The complex visibility of the radon transform. + @return freqs: The frequencies of the FFT of the radon transform. + @details The complex visibility is calculated by taking the FFT of the radon transform, shifting it to center the zero frequency, and selecting only the positive frequencies. + @note The padding factor is used to increase the resolution of the FFT. + @note The sample spacing is used to calculate the frequencies of the FFT. It corresponds to the pixel width in the original image. """ if Verbose: print("Calculating FFT...") @@ -56,16 +64,17 @@ def RadonToComplexVis( ) # 1D FFT of the projection radonshift = fftshift(radonff) # recenter FFT xfourier1 = fftshift( - fftfreq(PaddingFactor * FOORTRadon[0].shape[0], d=1) + fftfreq(PaddingFactor * FOORTRadon[0].shape[0], d=sample_spacing) ) # re centered frequencies indice1 = np.where((xfourier1 >= 0.0))[ 0 ] # select only the positive freqs the FFT is symmetrical anyway complvis = radonshift[:, indice1[0] : (indice1[-1] + 1)] + freqs = xfourier1[indice1[0] : (indice1[-1] + 1)] # frequencies of the visibilities if Verbose: print("Done calculating FFT.") - return complvis + return complvis, freqs def ComplexVisToNormVisAmp(ComplexVis: np.ndarray, Verbose: bool = True) -> np.ndarray: @@ -175,9 +184,15 @@ def FOORTToVisAmp( TruncateRange: tuple[float] = None, LimitRange: tuple[float] = (0.0, 100000.0), EquatPassesRange: tuple[int] = None, + EquatPassesSelection: list[int] = None, RadonPaddingFactor: float = 25, FileOutput: str = None, -) -> tuple[np.ndarray, str]: + LightRingRadius: float = None, + angular_size: float = 50.0 + * 1e-6 + * np.pi + / (180 * 60 * 60), # in radians (default = 50 microarcseconds) +) -> tuple[np.ndarray, np.ndarray, str]: """! @brief Converts FOORT output data to visibility amplitudes. @param FilePrefix: The prefix of the FOORT output files. @@ -189,19 +204,32 @@ def FOORTToVisAmp( @param TruncateRange: The range to truncate the data to (default = None). @param LimitRange: The range to limit the data to (default = (0.0, 100000.0)). @param EquatPassesRange: The range of equatorial passes to select (default = None). + @param EquatPassesSelection: The selection of equatorial passes to use (default = None). @param RadonPaddingFactor: The factor to pad the radon transform by (default = 25). @param FileOutput: The name of the file to save the visibility amplitudes to (default = None). + @param LightRingRadius: The radius of the light ring, if one is present (default = None). When used, this separets "inner" and "outer" photon rings. + @param angular_size: The angular size of the image in radians (default = 50 microarcseconds). @return visamps: The visibility amplitudes. + @return baselines: The baselines of the visibility amplitudes corresponding to the visamps. @return FirstLineInfo: The first line of the file, if applicable. """ - # Load in raw FOORT output data - FOORTData, FirstLineInfo = pyFOORT.LoadFOORTRawData( - FilePrefix, - "EquatorialEmission", - NrFiles=NrFiles, - FirstLineDescription=FirstLineDescription, - Verbose=Verbose, - ) + # If a light ring radius is given, we want to separate inner and outer photon rings. Inner photon rings get a minus sign, outer photon rings get a plus sign. + if LightRingRadius: + FOORTData, FirstLineInfo = pyFOORT.ModifiedEquatorialEmission( + FilePrefix, + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + LightRingRadius=LightRingRadius, + ) + else: + FOORTData, FirstLineInfo = pyFOORT.LoadFOORTRawData( + FilePrefix, + "EquatorialEmission", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) # Convert data to grid FOORTGrid = pyFOORT.DataToGrid( @@ -209,6 +237,7 @@ def FOORTToVisAmp( TruncateRange=TruncateRange, LimitRange=LimitRange, Diag2RangeSelect=EquatPassesRange, + Diag2AdvancedSelect=EquatPassesSelection, GridFraction=GridFraction, Verbose=Verbose, ) @@ -216,14 +245,20 @@ def FOORTToVisAmp( # Radon transform rad = RadonTransform(FOORTGrid, Angles, Verbose) # complex visibility - complvis = RadonToComplexVis(rad, PaddingFactor=RadonPaddingFactor, Verbose=Verbose) + complvis, freqs = RadonToComplexVis( + rad, + PaddingFactor=RadonPaddingFactor, + Verbose=Verbose, + sample_spacing=angular_size / FOORTGrid.shape[0], + ) + baselines = freqs / 1.0e9 # in Giga lambda # normalized visamp visamps = ComplexVisToNormVisAmp(complvis, Verbose=Verbose) if FileOutput: WriteVisAmpsToFile(visamps, FileOutput, FirstLineInfo=FirstLineInfo) - return visamps, FirstLineInfo + return visamps, baselines, FirstLineInfo # --- VISAMP ANALYSIS FUNCTIONS --- # diff --git a/PostProc/pyFOORT.py b/PostProc/pyFOORT.py index 84181cc..817bd3f 100644 --- a/PostProc/pyFOORT.py +++ b/PostProc/pyFOORT.py @@ -14,6 +14,7 @@ from matplotlib import cm # Color maps import matplotlib.colors as colors # Specific colors from scipy import ndimage # Used to map/distort background image +from scipy.interpolate import RectBivariateSpline # --- BASIC FUNCTIONS FOR ALL DIAGNOSTICS --- # @@ -124,6 +125,12 @@ def LoadFOORTRawData( df_list.append(new_data) data = pd.concat(df_list) + if data.isna().any().any(): + print( + "Warning: NaN values detected in data. These will be dropped before further processing." + ) + data.dropna(inplace=True) # drop any NaN values that may have appeared + # We are done loading! if Verbose: print("Done loading FOORT data.") @@ -138,6 +145,7 @@ def DataToGrid( TruncateRange: tuple[int] = None, LimitRange: tuple[float] = None, Diag2RangeSelect: tuple[int] = None, + Diag2AdvancedSelect: list[int] = None, Verbose: bool = True, ) -> np.ndarray: """! @@ -149,8 +157,10 @@ def DataToGrid( @param TruncateRange: Range to truncate data to (default None) @param LimitRange: Range to limit data to (default None) @param Diag2RangeSelect: Range of second diagnostic to select (default None) + @param Diag2AdvancedSelect: List of second diagnostic values to select (default None) @param Verbose: Whether to print progress information (default True) @return reshaped_data: Interpolated grid of specified size + @note: If both Diag2RangeSelect and Diag2AdvancedSelect are given, the advanced selection is used. """ if Verbose: print("Reshaping and interpolating grid...") @@ -183,14 +193,28 @@ def DataToGrid( RawData[RawData < TruncateRange[0]] = 0 # Limit the allowed range (e.g. for emission, if there is a minimum/maximum emission we want to allow) if LimitRange: + print("Limiting data to range " + str(LimitRange) + ".") RawData[RawData > LimitRange[1]] = LimitRange[1] RawData[RawData < LimitRange[0]] = LimitRange[0] + print("Max value after limiting: " + str(RawData.max())) # Only keep values of the first diagnostic (should be equatorial emission) # for pixels where the second diagnostic (=equatorial passes) is in a given range # note the abs(.) to take absolute value of the equatorial passes! - if Diag2RangeSelect and DiagToUse == 1: - RawDataDiag2 = np.abs(np.array(FOORTData["diag2"])) + if Diag2AdvancedSelect and Diag2RangeSelect: + print( + "Both a range and advanced selection are given for diagnostic 2. Continuing with the advanced selection, ignoring the range." + ) + if Diag2AdvancedSelect and DiagToUse == 1: + # Zero for all values that are not in the list, and all values in the list are set to 1 + RawDataDiag2 = np.array(FOORTData["diag2"]) + mask = np.zeros_like(RawDataDiag2) + for val in Diag2AdvancedSelect: + mask[RawDataDiag2 == val] = 1 + # Now we can simply multiply these two arrays to select the wanted pixels of the first diagnostic + RawData = np.multiply(RawData, mask) + elif Diag2RangeSelect and DiagToUse == 1: + RawDataDiag2 = np.array(FOORTData["diag2"]) # Zero out all values outside the range, and all set values inside the range to 1 RawDataDiag2[RawDataDiag2 < Diag2RangeSelect[0]] = 0 RawDataDiag2[RawDataDiag2 > Diag2RangeSelect[1]] = 0 @@ -212,6 +236,7 @@ def DisplayImage( ImageTitle: str = None, FileOutput: str = None, Verbose: bool = True, + Ax: plt.axes = None, ) -> None: """! @brief Display image from grid @@ -221,18 +246,23 @@ def DisplayImage( @param ImageTitle: Title of image (default None) @param FileOutput: File to save image to (default None) @param Verbose: Whether to print progress information (default True) + @param Ax: Axes to plot on (default None, creates new axes) """ if Verbose: print("Displaying image...") # plot the picture # leave room for title - if ImageTitle: - fig = plt.figure(figsize=(8, 10.3)) - plt.subplots_adjust(top=0.777) + if Ax == None: + if ImageTitle: + fig = plt.figure(figsize=(8, 10.3)) + plt.subplots_adjust(top=0.777) + else: + fig = plt.figure(figsize=(8, 8)) + if Ax: + ax = Ax else: - fig = plt.figure(figsize=(8, 8)) - ax = plt.axes() + ax = plt.axes() if ColorMinMax: ax.imshow( FOORTGrid, @@ -252,14 +282,16 @@ def DisplayImage( # Save the picture to file if applicable if FileOutput: - plt.savefig(FileOutput, format="pdf") + plt.tight_layout() + plt.savefig(FileOutput, format="png") if Verbose: print("Saved image to file " + FileOutput + ".") # Show the plot! - plt.show() - if Verbose: - print("Done displaying image.") + if Ax == None: + plt.show() + if Verbose: + print("Done displaying image.") # --- SPECIFIC GRID TO IMAGE AND COMBINATION FILE TO IMAGE FUNCTIONS PER DIAGNOSTIC --- # @@ -270,18 +302,27 @@ def GridToFourColorScreenImage( ImageTitle: str = None, FileOutput: str = None, Verbose: bool = True, + NoHorizon: bool = False, + Ax: plt.axes = None, ) -> None: """! @brief Convert grid to four-color screen image @param FOORTGrid: Grid data to display as image @param ImageTitle: Title of image (default None) @param FileOutput: File to save image to (default None) - @param Verbose: Whether to print progress information + @param Verbose: Whether to print progress information (default True) + @param NoHorizon: Set to True if the spacetime does not have a horizon, disabling the black color (default False) + @param Ax: Axes to plot on (default None, creates new axes) """ # Define our own color map for the four-color screen image - FourColorScreenColorMap = colors.ListedColormap( - ["black", "blue", "yellow", "red", "limegreen"] - ) + if NoHorizon: + FourColorScreenColorMap = colors.ListedColormap( + ["blue", "yellow", "red", "limegreen"] + ) + else: + FourColorScreenColorMap = colors.ListedColormap( + ["black", "blue", "yellow", "red", "limegreen"] + ) DisplayImage( FOORTGrid, @@ -289,6 +330,7 @@ def GridToFourColorScreenImage( ImageTitle=ImageTitle, FileOutput=FileOutput, Verbose=Verbose, + Ax=Ax ) @@ -300,6 +342,8 @@ def FOORTToFourColorScreenImage( Verbose: bool = False, GridFraction: float = 1, FileOutput: str = None, + NoHorizon: bool = False, + Ax: plt.axes = None, ) -> None: """! @brief Convert FOORT output to four-color screen image @@ -310,6 +354,8 @@ def FOORTToFourColorScreenImage( @param Verbose: Whether to print progress information (default False) @param GridFraction: Fraction of grid size to use (default 1) @param FileOutput: File to save image to (default None) + @param NoHorizon: Set to True if the spacetime does not have a horizon, disabling the black color (default False) + @param Ax: Axes to plot on (default None, creates new axes) """ # Load in raw FOORT output data FOORTData, FirstLineInfo = LoadFOORTRawData( @@ -325,7 +371,12 @@ def FOORTToFourColorScreenImage( FOORTGrid = DataToGrid(FOORTData, GridFraction=GridFraction, Verbose=Verbose) # Display image GridToFourColorScreenImage( - FOORTGrid, ImageTitle=FirstLineInfo, FileOutput=FileOutput, Verbose=Verbose + FOORTGrid, + ImageTitle=FirstLineInfo, + FileOutput=FileOutput, + Verbose=Verbose, + NoHorizon=NoHorizon, + Ax=Ax ) @@ -334,6 +385,7 @@ def GridToEquatorialPassesImage( ImageTitle: str = None, FileOutput: str = None, Verbose: bool = True, + Ax: plt.axes = None, ) -> None: """! @brief Convert grid to equatorial passes image @@ -341,6 +393,7 @@ def GridToEquatorialPassesImage( @param ImageTitle: Title of image (default None) @param FileOutput: File to save image to (default None) @param Verbose: Whether to print progress information (default True) + @param Ax: Axes to plot on (default None, creates new axes) """ # Color map for equatorial passes EquatorialPassesColorMap = cm.get_cmap( @@ -353,6 +406,7 @@ def GridToEquatorialPassesImage( ImageTitle=ImageTitle, FileOutput=FileOutput, Verbose=Verbose, + Ax=Ax ) @@ -368,6 +422,7 @@ def FOORTToEquatorialPassesImage( Verbose: bool = False, GridFraction: float = 1, FileOutput: str = None, + Ax: plt.axes = None, ) -> None: """! @brief Convert FOORT output to equatorial passes image @@ -382,6 +437,7 @@ def FOORTToEquatorialPassesImage( @param Verbose: Whether to print progress information (default False) @param GridFraction: Fraction of grid size to use (default 1) @param FileOutput: File to save image to (default None) + @param Ax: Axes to plot on (default None, creates new axes) """ # Load in raw FOORT output data FOORTData, FirstLineInfo = LoadFOORTRawData( @@ -421,9 +477,12 @@ def FOORTToEquatorialPassesImage( Verbose=Verbose, ) + if FOORTGrid.dtype != np.int64: + FOORTGrid = FOORTGrid.astype(np.int64) + # Display image GridToEquatorialPassesImage( - FOORTGrid, ImageTitle=FirstLineInfo, FileOutput=FileOutput, Verbose=Verbose + FOORTGrid, ImageTitle=FirstLineInfo, FileOutput=FileOutput, Verbose=Verbose, Ax=Ax ) @@ -432,6 +491,7 @@ def GridToEquatorialEmissionImage( ImageTitle: str = None, FileOutput: str = None, Verbose: bool = True, + Ax: plt.axes = None, ) -> None: """! @brief Convert grid to equatorial emission image @@ -439,6 +499,7 @@ def GridToEquatorialEmissionImage( @param ImageTitle: Title of image (default None) @param FileOutput: File to save image to (default None) @param Verbose: Whether to print progress information (default True) + @param Ax: Axes to plot on (default None, creates new axes) """ max_ring = np.max(FOORTGrid) min_ring = np.min(FOORTGrid) @@ -450,6 +511,7 @@ def GridToEquatorialEmissionImage( ImageTitle=ImageTitle, FileOutput=FileOutput, Verbose=Verbose, + Ax=Ax, ) @@ -463,7 +525,10 @@ def FOORTToEquatorialEmissionImage( TruncateRange: tuple[float] = None, LimitRange: tuple[float] = (0.0, 100000.0), EquatPassesRange: tuple[int] = None, + EquatPassesSelection: list[int] = None, FileOutput: str = None, + LightRingRadius: float = None, + Ax: plt.axes = None, ) -> None: """! @brief Convert FOORT output to equatorial emission image @@ -476,16 +541,28 @@ def FOORTToEquatorialEmissionImage( @param TruncateRange: Range to truncate data to (default None) @param LimitRange: Range to limit data to (default (0., 100000.)) @param EquatPassesRange: Range of equatorial passes to select (default None) + @param EquatPassesSelection: List of equatorial passes to select (default None) @param FileOutput: File to save image to (default None) + @param LightRingRadius: Radius of the light ring, if one is present (default None). When used, this separates "inner" and "outer" photon rings. + @param Ax: Axes to plot on (default None, creates new axes) """ # Load in raw FOORT output data - FOORTData, FirstLineInfo = LoadFOORTRawData( - FilePrefix, - "EquatorialEmission", - NrFiles=NrFiles, - FirstLineDescription=FirstLineDescription, - Verbose=Verbose, - ) + if LightRingRadius: + FOORTData, FirstLineInfo = ModifiedEquatorialEmission( + FilePrefix, + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + LightRingRadius=LightRingRadius, + ) + else: + FOORTData, FirstLineInfo = LoadFOORTRawData( + FilePrefix, + "EquatorialEmission", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) if DisplayImageTitle == False: FirstLineInfo = None @@ -495,13 +572,18 @@ def FOORTToEquatorialEmissionImage( TruncateRange=TruncateRange, LimitRange=LimitRange, Diag2RangeSelect=EquatPassesRange, + Diag2AdvancedSelect=EquatPassesSelection, GridFraction=GridFraction, Verbose=Verbose, ) # Display image GridToEquatorialEmissionImage( - FOORTGrid, ImageTitle=FirstLineInfo, FileOutput=FileOutput, Verbose=Verbose + FOORTGrid, + ImageTitle=FirstLineInfo, + FileOutput=FileOutput, + Verbose=Verbose, + Ax=Ax, ) @@ -644,3 +726,277 @@ def FOORTToDistortedBackground( FileOutput=FileOutput, Verbose=Verbose, ) + + +def ModifiedEquatorialEmission( + FilePrefix: str, + NrFiles: int = 1, + FirstLineDescription: bool = True, + Verbose: bool = False, + LightRingRadius: float = 3.0, +) -> tuple[pd.DataFrame, str]: + """! + @brief Load FOORT output data and modify equatorial emission data to include closest radius mask + @param FilePrefix: Prefix of the FOORT output files + @param NrFiles: Number of files to load (default 1) + @param FirstLineDescription: Whether the first line of the file contains information (default True). + @param Verbose: Whether to print progress information (default False) + @param LightRingRadius: Radius of the light ring (default 3.0, Schwarzschild value) + @return EquatorialEmission_df: DataFrame containing equatorial emission data with closest radius mask applied + @return FirstLineInfo: Information contained in first line of the file + @note: This function modifies the equatorial emission data to include a mask based on the closest radius encountered: if this is smaller than the light ring radius, the number of equatorial passes gets multiplied with a minus sign. + The mask is applied to the second diagnostic (diag2) of the equatorial emission data. + """ + EquatorialEmission_df, FirstLineInfo = LoadFOORTRawData( + FilePrefix, + "EquatorialEmission", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) + ClosestRadius_df, _ = LoadFOORTRawData( + FilePrefix, + "ClosestRadius", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) + ClosestRadius_mask = (ClosestRadius_df["diag1"] >= LightRingRadius) * 1 - ( + ClosestRadius_df["diag1"] < LightRingRadius + ) * 1 + EquatorialEmission_df["diag2"] = EquatorialEmission_df["diag2"] * ClosestRadius_mask + + return EquatorialEmission_df, FirstLineInfo + + +def FOORTToEquatorialEmissionLineout( + FilePrefix: str, + NrFiles: int = 1, + FirstLineDescription: bool = True, + DisplayImageTitle: bool = False, + Verbose: bool = False, + GridFraction: float = 1, + TruncateRange: tuple[float] = None, + LimitRange: tuple[float] = (0.0, 100000.0), + EquatPassesRange: tuple[int] = None, + EquatPassesSelection: list[int] = None, + FileOutput: str = None, + LightRingRadius: float = None, + Ax: plt.axes = None, + angle: float = 20.0, + ipl_points: int = 1000, + x_scale = 100, # fiducial scale in micro as, + return_lineout: bool = False, + **plot_kwargs, +) -> None | tuple[np.ndarray, np.ndarray]: + """! + @brief Convert FOORT output to equatorial emission image + @param FilePrefix: Prefix of the FOORT output files + @param NrFiles: Number of files to load (default 1) + @param FirstLineDescription: Whether the first line of the file contains information (default True) + @param DisplayImageTitle: Whether to display the image title (default False) + @param Verbose: Whether to print progress information (default False) + @param GridFraction: Fraction of grid size to use (default 1) + @param TruncateRange: Range to truncate data to (default None) + @param LimitRange: Range to limit data to (default (0., 100000.)) + @param EquatPassesRange: Range of equatorial passes to select (default None) + @param FileOutput: File to save image to (default None) + @param LightRingRadius: Radius of the light ring (default None) + @param Ax: Matplotlib axes to plot on (default None) + @param angle: Angle to plot the emission at (default 20.) + @param ipl_points: Number of points to interpolate (default 1000) + @param x_scale: Scale of the x-axis in microarcseconds (default 100) + @param return_lineout: Whether to return the lineout data as arrays (default False) + """ + + if not (min(abs(angle), abs(angle - 180)) <= 45): + raise ValueError( + "At the moment, only angles between -45 and 45 degrees are supported." + ) + + # Load in raw FOORT output data + if LightRingRadius: + FOORTData, FirstLineInfo = ModifiedEquatorialEmission( + FilePrefix, + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + LightRingRadius=LightRingRadius, + ) + else: + FOORTData, FirstLineInfo = LoadFOORTRawData( + FilePrefix, + "EquatorialEmission", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) + if DisplayImageTitle == False: + FirstLineInfo = None + + # Convert data to grid + FOORTGrid = DataToGrid( + FOORTData, + TruncateRange=TruncateRange, + LimitRange=LimitRange, + Diag2RangeSelect=EquatPassesRange, + Diag2AdvancedSelect=EquatPassesSelection, + GridFraction=GridFraction, + Verbose=Verbose, + ) + + x_dim = FOORTGrid.shape[0] + y_dim = FOORTGrid.shape[1] + intpl = RectBivariateSpline(range(x_dim), range(y_dim), FOORTGrid) + if min(abs(angle), abs(angle - 180)) <= 45: + x = np.linspace(0, x_dim, ipl_points) + y = np.tan(angle * np.pi / 180) * (x - x_dim // 2) + y_dim // 2 + + I = np.diag(intpl(x, y)) + + if Ax is None: + fig, ax = plt.subplots(figsize=(8, 8)) + else: + ax = Ax + + slice_length = x_scale / np.cos(angle * np.pi / 180) + + distance = x * slice_length/ x_dim + ax.plot(distance, I, linewidth=1.0, **plot_kwargs) + + # Save the picture to file if applicable + if FileOutput and (Ax is None): + fig.savefig(FileOutput, format="pdf") + if Verbose: + print("Saved image to file " + FileOutput + ".") + elif FileOutput: + print("Not saving figure as ax is given externally.") + + if return_lineout: + return distance, I + +def GridToClosestRadiusImage( + FOORTGrid: np.ndarray, + ImageTitle: str = None, + FileOutput: str = None, + Verbose: bool = True, + Ax: plt.axes = None, + logscale: bool = False, + contours: list[float] | None = None, + MaskFloorForContours: bool = True, +) -> None: + """! + @brief Convert grid to closest radius image + @param FOORTGrid: Grid data to display as image + @param ImageTitle: Title of image (default None) + @param FileOutput: File to save image to (default None) + @param Verbose: Whether to print progress information (default True) + @param Ax: Axes to plot on (default None, creates new axes) + """ + # Color map for closest radius: always scale from 0 to the max value + ClosestRadiusColorMap = "viridis" + if logscale: + FOORTGrid = np.log10(FOORTGrid + 1e-10) # add small value to avoid log(0) + min_radius = np.min(FOORTGrid) + max_radius = np.max(FOORTGrid) + print("Closest radius values range from " + str(min_radius) + " to " + str(max_radius) + ".") + + DisplayImage( + FOORTGrid, + ClosestRadiusColorMap, + ColorMinMax=(min_radius, max_radius), + ImageTitle=ImageTitle, + FileOutput=FileOutput, + Verbose=Verbose, + Ax=Ax + ) + + # Add contour lines for given radii + if contours is not None: + if Ax is None: + print("Cannot add contours if ax is None, as the figure and axes are created in the DisplayImage function. Not adding contours.") + return + + contour_grid = FOORTGrid + if MaskFloorForContours: + # Mask clipped floor values to avoid spurious contours around masked regions. + contour_grid = np.ma.masked_where( + np.isclose(FOORTGrid, min_radius, rtol=0.0, atol=1e-12), + FOORTGrid, + ) + + CS = Ax.contour(contour_grid, levels=contours, colors='white', linewidths=0.5) + Ax.clabel(CS, inline=True, fontsize=8, fmt='%.2f') + if FileOutput: + print("Not saving figure with contours as ax is given externally.") + + + + +def FOORTToClosestRadius( + FilePrefix: str, + NrFiles: int = 1, + FirstLineDescription: bool = True, + DisplayImageTitle: bool = False, + Verbose: bool = False, + GridFraction: float = 1, + TruncateRange: tuple[float] = None, + LimitRange: tuple[float] = None, + Diag2RangeSelect: tuple[int] = None, + Diag2AdvancedSelect: list[int] = None, + FileOutput: str = None, + Ax: plt.axes = None, + logscale: bool = False, + contours: list[float] | None = None, + MaskFloorForContours: bool = True, +) -> None: + """! + @brief Convert FOORT output to closest radius image + @param FilePrefix: Prefix of the FOORT output files + @param NrFiles: Number of files to load (default 1) + @param FirstLineDescription: Whether the first line of the file contains information (default True) + @param DisplayImageTitle: Whether to display the image title (default False) + @param Verbose: Whether to print progress information (default False) + @param GridFraction: Fraction of grid size to use (default 1) + @param TruncateRange: Range to truncate data to (default None) + @param LimitRange: Range to limit data to (default None) + @param Diag2RangeSelect: Range of second diagnostic to select (default None) + @param Diag2AdvancedSelect: List of second diagnostic values to select (default None) + @param FileOutput: File to save image to (default None) + @param Ax: Matplotlib axes to plot on (default None) + @param logscale: Whether to apply log scale to the closest radius values (default False) + """ + + # Load in raw FOORT output data + FOORTData, FirstLineInfo = LoadFOORTRawData( + FilePrefix, + "ClosestRadius", + NrFiles=NrFiles, + FirstLineDescription=FirstLineDescription, + Verbose=Verbose, + ) + if DisplayImageTitle == False: + FirstLineInfo = None + + # Convert data to grid + FOORTGrid = DataToGrid( + FOORTData, + TruncateRange=TruncateRange, + LimitRange=LimitRange, + Diag2RangeSelect=Diag2RangeSelect, + Diag2AdvancedSelect=Diag2AdvancedSelect, + GridFraction=GridFraction, + Verbose=Verbose, + ) + + # Display image + GridToClosestRadiusImage( + FOORTGrid, + ImageTitle=FirstLineInfo, + FileOutput=FileOutput, + Verbose=Verbose, + Ax=Ax, + logscale=logscale, + contours=contours, + MaskFloorForContours=MaskFloorForContours, + ) \ No newline at end of file diff --git a/README.md b/README.md index 9fd174a..50d8ce6 100644 --- a/README.md +++ b/README.md @@ -36,15 +36,23 @@ Recent versions of MacOS may encounter issues with CMake. The following seems to (Found in [this issue](https://gist.github.com/scivision/d69faebbc56da9714798087b56de925a)) ``` +export CC=/opt/homebrew/bin/gcc-14 export CXX=/opt/homebrew/bin/g++-14 +export FC=/opt/homebrew/bin/gfortran-14 export SDKROOT=/Library/Developer/CommandLineTools/SDKs/MacOSX14.sdk/ ``` -For convenience, this is stored in `alternative_builds/MacOS_workaround.sh`, which should be executed with `source \the\path\` in the terminal in case issues are encountered. You may need to `rm -rf build` in order to start clean. +For convenience, this is stored in `alternative_builds/MacOS_workaround.sh`, which should be executed with `source \the\path\` in the terminal in case issues are encountered. You may need to `rm -rf build` in order to start clean. In case this doesn't work (which is not unlikely - some of it is black magic), one can also revert to the `makefile_mac` as described below. #### Old MakeFiles -Alternatively, old Makefiles can be found in `alternative_builds/old_makefiles`. These can be system dependent, however, but should be adaptable to your needs. +Alternatively, old Makefiles can be found in `alternative_builds/old_makefiles`. These can be system dependent, however, but should be adaptable to your needs. If the user wants to use these, the relevant makefile should be copied into the `src` folder, in which + +``` +make -f +``` + +should be run, where `` should be replaced with the name of the relevant file. ### WINDOWS VISUAL STUDIO diff --git a/alternative_builds/MacOS_gcc15_workaround.sh b/alternative_builds/MacOS_gcc15_workaround.sh new file mode 100644 index 0000000..992ade7 --- /dev/null +++ b/alternative_builds/MacOS_gcc15_workaround.sh @@ -0,0 +1,4 @@ +export CC=/opt/homebrew/bin/gcc-15 +export CXX=/opt/homebrew/bin/g++-15 +export FC=/opt/homebrew/bin/gfortran-15 +export SDKROOT=/Library/Developer/CommandLineTools/SDKs/MacOSX15.sdk/