From ba7812c243390e85409d89d409fd3d28793a3fd5 Mon Sep 17 00:00:00 2001 From: Alex-J-Brown Date: Thu, 20 Aug 2026 00:10:45 +0200 Subject: [PATCH 1/4] first go at integrating Grid type --- src/binary_model.rs | 157 ++++++++++++++++++++++---------------------- src/grid.rs | 42 ++++++++++++ src/lib.rs | 4 +- 3 files changed, 124 insertions(+), 79 deletions(-) create mode 100644 src/grid.rs diff --git a/src/binary_model.rs b/src/binary_model.rs index cb48ffe..e35e981 100644 --- a/src/binary_model.rs +++ b/src/binary_model.rs @@ -2,6 +2,7 @@ use crate::comp_gravity::{comp_gravity1, comp_gravity2}; use crate::comp_light::{comp_bright_spot, comp_disc, comp_disc_edge, comp_star1, comp_star2}; use crate::comp_radius::comp_radius; use crate::ginterp::Ginterp; +use crate::grid::Grid; use crate::ldc::LDC; use crate::model::{Entry, Model, ModelUpdate}; use crate::set_bright_spot_grid::set_bright_spot_grid; @@ -78,19 +79,19 @@ pub struct LightCurve { #[pyclass] pub struct BinaryModel { #[pyo3(get)] - star1_coarse_grid: Vec, + star1_coarse_grid: Grid, #[pyo3(get)] - star2_coarse_grid: Vec, + star2_coarse_grid: Grid, #[pyo3(get)] - star1_fine_grid: Vec, + star1_fine_grid: Grid, #[pyo3(get)] - star2_fine_grid: Vec, + star2_fine_grid: Grid, #[pyo3(get)] - disc_grid: Vec, + disc_grid: Grid, #[pyo3(get)] - disc_edge_grid: Vec, + disc_edge_grid: Grid, #[pyo3(get)] - bright_spot_grid: Vec, + bright_spot_grid: Grid, gint: Ginterp, rlens1: f64, model_beaming1: bool, @@ -202,30 +203,30 @@ impl BinaryModel { /// "disc_edge", /// "bright_spot" /// - pub fn set_grid_fluxes(&mut self, grid: &str, fluxes: Vec) -> Result<(), RocheError> { - let chosen_grid = match grid { - "star1_fine" => &mut self.star1_fine_grid, - "star1_coarse" => &mut self.star1_coarse_grid, - "star2_fine" => &mut self.star2_fine_grid, - "star2_coarse" => &mut self.star2_coarse_grid, - "disc" => &mut self.disc_grid, - "disc_edge" => &mut self.disc_edge_grid, - "bright_spot" => &mut self.bright_spot_grid, - _ => return Err(RocheError::ParameterError("Not a valid grid.".to_string())), - }; - - apply_fluxes(chosen_grid, fluxes)?; - self.gint = set_ginterp( - &self.model, - self.rlens1, - &self.star1_coarse_grid, - &self.star2_coarse_grid, - &self.star1_fine_grid, - &self.star2_fine_grid, - )?; - - Ok(()) - } + // pub fn set_grid_fluxes(&mut self, grid: &str, fluxes: Vec) -> Result<(), RocheError> { + // let chosen_grid = match grid { + // "star1_fine" => &mut self.star1_fine_grid, + // "star1_coarse" => &mut self.star1_coarse_grid, + // "star2_fine" => &mut self.star2_fine_grid, + // "star2_coarse" => &mut self.star2_coarse_grid, + // "disc" => &mut self.disc_grid, + // "disc_edge" => &mut self.disc_edge_grid, + // "bright_spot" => &mut self.bright_spot_grid, + // _ => return Err(RocheError::ParameterError("Not a valid grid.".to_string())), + // }; + + // apply_fluxes(chosen_grid, fluxes)?; + // self.gint = set_ginterp( + // &self.model, + // self.rlens1, + // &self.star1_coarse_grid, + // &self.star2_coarse_grid, + // &self.star1_fine_grid, + // &self.star2_fine_grid, + // )?; + + // Ok(()) + // } #[pyo3(signature = ( time, @@ -364,19 +365,19 @@ impl BinaryModel { }; let (logg1, logg2) = if self.model.velocity_scale.defined { - let logg1 = comp_gravity1(&self.model, &self.star1_fine_grid)?; - let logg2 = comp_gravity2(&self.model, &self.star2_fine_grid)?; + let logg1 = comp_gravity1(&self.model, &self.star1_fine_grid.points)?; + let logg2 = comp_gravity2(&self.model, &self.star2_fine_grid.points)?; (Some(logg1), Some(logg2)) } else { (None, None) }; let rva1: f64 = if self.model.roche1 { - comp_radius(&self.star1_coarse_grid, Star::Primary) + comp_radius(&self.star1_coarse_grid.points, Star::Primary) } else { self.model.r1.value }; - let rva2: f64 = comp_radius(&self.star2_coarse_grid, Star::Secondary); + let rva2: f64 = comp_radius(&self.star2_coarse_grid.points, Star::Secondary); Ok(LightCurve { star1: star1.into_pyarray(py).unbind(), @@ -410,8 +411,8 @@ impl BinaryModel { self.model.velocity_scale.value, self.model_beaming1, &self.gint, - &self.star1_fine_grid, - &self.star1_coarse_grid, + &self.star1_fine_grid.points, + &self.star1_coarse_grid.points, ) } @@ -429,8 +430,8 @@ impl BinaryModel { self.model.glens1, self.rlens1, &self.gint, - &self.star2_fine_grid, - &self.star2_coarse_grid, + &self.star2_fine_grid.points, + &self.star2_coarse_grid.points, ) } @@ -442,7 +443,7 @@ impl BinaryModel { phase, expose, n_div, - &self.disc_grid, + &self.disc_grid.points, ) } @@ -454,7 +455,7 @@ impl BinaryModel { phase, expose, n_div, - &self.disc_edge_grid, + &self.disc_edge_grid.points, ) } @@ -464,7 +465,7 @@ impl BinaryModel { phase, expose, n_div, - &self.bright_spot_grid, + &self.bright_spot_grid.points, ) } @@ -484,22 +485,22 @@ impl BinaryModel { set_star_continuum( &self.model, - &mut self.star1_fine_grid, - &mut self.star2_fine_grid, + &mut self.star1_fine_grid.points, + &mut self.star2_fine_grid.points, )?; set_star_continuum( &self.model, - &mut self.star1_coarse_grid, - &mut self.star2_coarse_grid, + &mut self.star1_coarse_grid.points, + &mut self.star2_coarse_grid.points, )?; self.gint = set_ginterp( &self.model, self.rlens1, - &self.star1_coarse_grid, - &self.star2_coarse_grid, - &self.star1_fine_grid, - &self.star2_fine_grid, + &self.star1_coarse_grid.points, + &self.star2_coarse_grid.points, + &self.star1_fine_grid.points, + &self.star2_fine_grid.points, )?; if self.model.add_disc { @@ -515,7 +516,7 @@ impl BinaryModel { self.model.temp_disc.value, self.model.texp_disc.value, self.model.wavelength, - &mut self.disc_grid, + &mut self.disc_grid.points, ); // Set the surface brightness of outer edge, accounting for @@ -526,12 +527,12 @@ impl BinaryModel { self.model.t2.value.abs(), self.model.absorb_edge.value, self.model.wavelength, - &mut self.disc_edge_grid, + &mut self.disc_edge_grid.points, ); } if self.model.add_spot { - self.bright_spot_grid = set_bright_spot_grid(&self.model)?; + self.bright_spot_grid.points = set_bright_spot_grid(&self.model)?; } Ok(()) } @@ -541,13 +542,13 @@ fn build_grids( model: &Model, ) -> Result< ( - Vec, - Vec, - Vec, - Vec, - Vec, - Vec, - Vec, + Grid, + Grid, + Grid, + Grid, + Grid, + Grid, + Grid, Ginterp, f64, bool, @@ -714,13 +715,13 @@ fn build_grids( } Ok(( - star1_coarse_grid, - star2_coarse_grid, - star1_fine_grid, - star2_fine_grid, - disc_grid, - disc_edge_grid, - bright_spot_grid, + Grid{points: star1_coarse_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: star2_coarse_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: star1_fine_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: star2_fine_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: disc_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: disc_edge_grid, q: model.q.value, iangle: model.iangle.value}, + Grid{points: bright_spot_grid, q: model.q.value, iangle: model.iangle.value}, gint, rlens1, model_beaming1, @@ -982,14 +983,14 @@ pub fn chisq_log_prob( (chisq_sum, log_prob) } -fn apply_fluxes(points: &mut Vec, fluxes: Vec) -> Result<(), RocheError> { - if points.len() != fluxes.len() { - return Err(RocheError::ParameterError( - "Selected grid and flux array have mismatched lengths.".to_string(), - )); - } - for (point, flux) in points.iter_mut().zip(fluxes) { - point.set_flux(flux); - } - Ok(()) -} +// fn apply_fluxes(points: &mut Vec, fluxes: Vec) -> Result<(), RocheError> { +// if points.len() != fluxes.len() { +// return Err(RocheError::ParameterError( +// "Selected grid and flux array have mismatched lengths.".to_string(), +// )); +// } +// for (point, flux) in points.iter_mut().zip(fluxes) { +// point.set_flux(flux); +// } +// Ok(()) +// } diff --git a/src/grid.rs b/src/grid.rs new file mode 100644 index 0000000..70f5b75 --- /dev/null +++ b/src/grid.rs @@ -0,0 +1,42 @@ +use roche::{self, Point, Vec3}; +use std::f64::consts::TAU; +use pyo3::prelude::*; + +#[pyclass(skip_from_py_object)] +#[derive(Clone, Debug)] +pub struct Grid { + #[pyo3(get)] + pub points: Vec, + #[pyo3(get)] + pub q: f64, + #[pyo3(get)] + pub iangle: f64, +} + +#[pymethods] +impl Grid { + #[pyo3(signature = (phase, iangle=None))] + pub fn project_2d(&self, phase: f64, iangle: Option) -> (Vec, Vec) { + let iangle = match iangle { + Some(iangle) => iangle, + None => self.iangle, + }; + let mut x_arr: Vec = vec![]; + let mut y_arr: Vec = vec![]; + let cofm = Vec3::new(self.q/(1.0+self.q), 0.0, 0.0); + let (sinp, cosp) = (TAU*phase).sin_cos(); + let earth = roche::set_earth_iangle(iangle, phase); + let xsky = Vec3::new(sinp, cosp, 0.0); + let ysky = earth.cross(&xsky); + + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + let r = point.position - cofm; + x_arr.push(r.dot(&xsky)); + y_arr.push(r.dot(&ysky)); + } + } + (x_arr, y_arr) + } + +} diff --git a/src/lib.rs b/src/lib.rs index 1a0bd6c..d967687 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,12 +1,13 @@ use pyo3::prelude::*; -use crate::ldc::LDCType; +use crate::{grid::Grid, ldc::LDCType}; pub mod binary_model; pub mod comp_gravity; pub mod comp_light; pub mod comp_radius; pub mod ginterp; +pub mod grid; pub mod ldc; pub mod model; pub mod numface; @@ -23,5 +24,6 @@ fn lcurve(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> { m.add_class::()?; m.add_class::()?; m.add_class::()?; + m.add_class::()?; Ok(()) } From bfd0d2f5dadfae6cbc8bdb1e9ca7ddfae2adef69 Mon Sep 17 00:00:00 2001 From: Alex-J-Brown Date: Thu, 20 Aug 2026 14:56:36 +0200 Subject: [PATCH 2/4] Clean up implementation of Grid struct --- src/binary_model.rs | 82 ++++++++++++++++++------------------- src/comp_gravity.rs | 11 ++--- src/comp_light.rs | 47 ++++++++++----------- src/grid.rs | 39 ++++++++++++++---- src/set_bright_spot_grid.rs | 5 ++- src/set_disc_continuum.rs | 12 +++--- src/set_disc_grid.rs | 11 ++--- src/set_star_continuum.rs | 12 +++--- src/set_star_grid.rs | 15 ++++--- 9 files changed, 133 insertions(+), 101 deletions(-) diff --git a/src/binary_model.rs b/src/binary_model.rs index e35e981..4c915f9 100644 --- a/src/binary_model.rs +++ b/src/binary_model.rs @@ -16,7 +16,7 @@ use pyo3::types::{PyDict, PyDictMethods}; use rayon::prelude::*; use roche::constants::{C, DAY}; use roche::errors::RocheError; -use roche::{self, Etype, Point, Star, disc_eclipse}; +use roche::{self, Etype, Star, disc_eclipse}; use serde_pyobject::from_pyobject; use std::collections::HashMap; use std::f64::consts::TAU; @@ -365,8 +365,8 @@ impl BinaryModel { }; let (logg1, logg2) = if self.model.velocity_scale.defined { - let logg1 = comp_gravity1(&self.model, &self.star1_fine_grid.points)?; - let logg2 = comp_gravity2(&self.model, &self.star2_fine_grid.points)?; + let logg1 = comp_gravity1(&self.model, &self.star1_fine_grid)?; + let logg2 = comp_gravity2(&self.model, &self.star2_fine_grid)?; (Some(logg1), Some(logg2)) } else { (None, None) @@ -411,8 +411,8 @@ impl BinaryModel { self.model.velocity_scale.value, self.model_beaming1, &self.gint, - &self.star1_fine_grid.points, - &self.star1_coarse_grid.points, + &self.star1_fine_grid, + &self.star1_coarse_grid, ) } @@ -430,8 +430,8 @@ impl BinaryModel { self.model.glens1, self.rlens1, &self.gint, - &self.star2_fine_grid.points, - &self.star2_coarse_grid.points, + &self.star2_fine_grid, + &self.star2_coarse_grid, ) } @@ -443,7 +443,7 @@ impl BinaryModel { phase, expose, n_div, - &self.disc_grid.points, + &self.disc_grid, ) } @@ -455,7 +455,7 @@ impl BinaryModel { phase, expose, n_div, - &self.disc_edge_grid.points, + &self.disc_edge_grid, ) } @@ -465,7 +465,7 @@ impl BinaryModel { phase, expose, n_div, - &self.bright_spot_grid.points, + &self.bright_spot_grid, ) } @@ -485,22 +485,22 @@ impl BinaryModel { set_star_continuum( &self.model, - &mut self.star1_fine_grid.points, - &mut self.star2_fine_grid.points, + &mut self.star1_fine_grid, + &mut self.star2_fine_grid, )?; set_star_continuum( &self.model, - &mut self.star1_coarse_grid.points, - &mut self.star2_coarse_grid.points, + &mut self.star1_coarse_grid, + &mut self.star2_coarse_grid, )?; self.gint = set_ginterp( &self.model, self.rlens1, - &self.star1_coarse_grid.points, - &self.star2_coarse_grid.points, - &self.star1_fine_grid.points, - &self.star2_fine_grid.points, + &self.star1_coarse_grid, + &self.star2_coarse_grid, + &self.star1_fine_grid, + &self.star2_fine_grid, )?; if self.model.add_disc { @@ -516,7 +516,7 @@ impl BinaryModel { self.model.temp_disc.value, self.model.texp_disc.value, self.model.wavelength, - &mut self.disc_grid.points, + &mut self.disc_grid, ); // Set the surface brightness of outer edge, accounting for @@ -527,12 +527,12 @@ impl BinaryModel { self.model.t2.value.abs(), self.model.absorb_edge.value, self.model.wavelength, - &mut self.disc_edge_grid.points, + &mut self.disc_edge_grid, ); } if self.model.add_spot { - self.bright_spot_grid.points = set_bright_spot_grid(&self.model)?; + self.bright_spot_grid = set_bright_spot_grid(&self.model)?; } Ok(()) } @@ -559,8 +559,8 @@ fn build_grids( model.validate()?; let mut star1_fine_grid = set_star_grid(model, Star::Primary, true)?; let mut star2_fine_grid = set_star_grid(model, Star::Secondary, true)?; - let mut star1_coarse_grid: Vec; - let mut star2_coarse_grid: Vec; + let mut star1_coarse_grid: Grid; + let mut star2_coarse_grid: Grid; let (r1, mut r2) = model.get_r1r2(); let rl2: f64 = 1.0 - roche::x_l1_2(model.q.value, model.spin2.value)?; @@ -591,9 +591,9 @@ fn build_grids( set_star_continuum(model, &mut star1_coarse_grid, &mut star2_coarse_grid)?; } - let mut disc_grid: Vec = vec![]; - let mut disc_edge_grid: Vec = vec![]; - let mut bright_spot_grid: Vec = vec![]; + let mut disc_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); + let mut disc_edge_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); + let mut bright_spot_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); let mut rlens1 = 0.0; if model.glens1 { @@ -632,7 +632,7 @@ fn build_grids( let mut eclipses: Etype; if model.opaque { - for point in &mut star1_fine_grid { + for point in &mut star1_fine_grid.points { eclipses = disc_eclipse( model.iangle.value, rdisc1, @@ -646,7 +646,7 @@ fn build_grids( } } - for point in &mut star1_coarse_grid { + for point in &mut star1_coarse_grid.points { eclipses = disc_eclipse( model.iangle.value, rdisc1, @@ -660,7 +660,7 @@ fn build_grids( } } - for point in &mut star2_fine_grid { + for point in &mut star2_fine_grid.points { eclipses = disc_eclipse( model.iangle.value, rdisc1, @@ -674,7 +674,7 @@ fn build_grids( } } - for point in &mut star2_coarse_grid { + for point in &mut star2_coarse_grid.points { eclipses = disc_eclipse( model.iangle.value, rdisc1, @@ -715,13 +715,13 @@ fn build_grids( } Ok(( - Grid{points: star1_coarse_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: star2_coarse_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: star1_fine_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: star2_fine_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: disc_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: disc_edge_grid, q: model.q.value, iangle: model.iangle.value}, - Grid{points: bright_spot_grid, q: model.q.value, iangle: model.iangle.value}, + star1_coarse_grid, + star2_coarse_grid, + star1_fine_grid, + star2_fine_grid, + disc_grid, + disc_edge_grid, + bright_spot_grid, gint, rlens1, model_beaming1, @@ -732,10 +732,10 @@ fn build_grids( pub fn set_ginterp( model: &Model, rlens1: f64, - star1c: &Vec, - star2c: &Vec, - star1f: &Vec, - star2f: &Vec, + star1c: &Grid, + star2c: &Grid, + star1f: &Grid, + star2f: &Grid, ) -> Result { let (r1, mut r2) = model.get_r1r2(); let rl2: f64 = 1.0 - roche::x_l1_2(model.q.value, model.spin2.value)?; diff --git a/src/comp_gravity.rs b/src/comp_gravity.rs index 4107037..1bf5e25 100644 --- a/src/comp_gravity.rs +++ b/src/comp_gravity.rs @@ -1,7 +1,8 @@ +use crate::grid::Grid; use crate::model::Model; use roche::constants::DAY; use roche::errors::RocheError; -use roche::{Point, RocheContext, Star, Vec3}; +use roche::{RocheContext, Star, Vec3}; use std::f64::consts::TAU; // @@ -14,7 +15,7 @@ use std::f64::consts::TAU; // \return the value of logg // -pub fn comp_gravity1(model: &Model, star1_fine_grid: &Vec) -> Result { +pub fn comp_gravity1(model: &Model, star1_fine_grid: &Grid) -> Result { // Calculate the unit scaling factor to get CGS gravity let gm1m2: f64 = (1000.0 * model.velocity_scale.value).powi(3) * model.tperiod * DAY / TAU; let a: f64 = (gm1m2 / (TAU / DAY / model.tperiod).powi(2)).powf(1.0 / 3.0); @@ -44,7 +45,7 @@ pub fn comp_gravity1(model: &Model, star1_fine_grid: &Vec) -> Result) -> Result) -> Result { +pub fn comp_gravity2(model: &Model, star2_fine_grid: &Grid) -> Result { // Calculate the unit scaling factor to get CGS gravity let gm1m2: f64 = (1000.0 * model.velocity_scale.value).powi(3) * model.tperiod * DAY / TAU; let a: f64 = (gm1m2 / (TAU / DAY / model.tperiod).powi(2)).powf(1.0 / 3.0); @@ -99,7 +100,7 @@ pub fn comp_gravity2(model: &Model, star2_fine_grid: &Vec) -> Result, - star2f: &Vec, - star1c: &Vec, - star2c: &Vec, + star1f: &Grid, + star2f: &Grid, + star1c: &Grid, + star2c: &Grid, ) -> f64 { let x_cofm: f64 = q / (1.0 + q); let (sini, cosi) = iangle.to_radians().sin_cos(); @@ -98,14 +99,14 @@ pub fn comp_light( earth = roche::set_earth(cosi, sini, phi); ptype = gint.interp_type(phi); - let star1: &Vec = if ptype == 1 { star1f } else { star1c }; + let star1: &Grid = if ptype == 1 { star1f } else { star1c }; - let star2: &Vec = if ptype == 3 { star2f } else { star2c }; + let star2: &Grid = if ptype == 3 { star2f } else { star2c }; ssum = 0.0; // star 1 - for point in star1 { + for point in &star1.points { if point.is_visible(phi) { mu = earth.dot(&point.direction); if ldc1.see(mu) { @@ -128,7 +129,7 @@ pub fn comp_light( // star 2 ssum2 = 0.0; - for point in star2 { + for point in &star2.points { if point.is_visible(phi) { mu = earth.dot(&point.direction); @@ -201,8 +202,8 @@ pub fn comp_star1( vscale: f64, model_beaming1: bool, gint: &Ginterp, - star1f: &Vec, - star1c: &Vec, + star1f: &Grid, + star1c: &Grid, ) -> f64 { let x_cofm: f64 = q / (1.0 + q); let (sini, cosi) = iangle.to_radians().sin_cos(); @@ -242,12 +243,12 @@ pub fn comp_star1( // Define the grid to use ptype = gint.interp_type(phi); - let star1: &Vec = if ptype == 1 { star1f } else { star1c }; + let star1: &Grid = if ptype == 1 { star1f } else { star1c }; ssum = 0.0; let phi_normed: f64 = phi - phi.floor(); // star 1 - for point in star1 { + for point in &star1.points { if point.is_visible_phase_normed(phi_normed) { mu = earth.dot(&point.direction); if ldc1.see(mu) { @@ -285,8 +286,8 @@ pub fn comp_star2( glens1: bool, rlens1: f64, gint: &Ginterp, - star2f: &Vec, - star2c: &Vec, + star2f: &Grid, + star2c: &Grid, ) -> f64 { let x_cofm: f64 = q / (1.0 + q); let (sini, cosi) = iangle.to_radians().sin_cos(); @@ -333,12 +334,12 @@ pub fn comp_star2( // Define the grid to use ptype = gint.interp_type(phi); - let star2: &Vec = if ptype == 3 { star2f } else { star2c }; + let star2: &Grid = if ptype == 3 { star2f } else { star2c }; ssum = 0.0; let phi_normed: f64 = phi - phi.floor(); // star 2 - for point in star2 { + for point in &star2.points { if point.is_visible_phase_normed(phi_normed) { mu = earth.dot(&point.direction); if ldc2.see(mu) { @@ -404,7 +405,7 @@ pub fn comp_disc( phase: f64, expose: f64, n_div: i32, - disc_grid: &Vec, + disc_grid: &Grid, ) -> f64 { let ri = iangle.to_radians(); let (sini, cosi) = ri.sin_cos(); @@ -435,7 +436,7 @@ pub fn comp_disc( ssum = 0.0; let phi_normed: f64 = phi - phi.floor(); // Disc - for point in disc_grid { + for point in &disc_grid.points { mu = earth.dot(&point.direction); if mu > 0.0 && point.is_visible_phase_normed(phi_normed) { ssum += mu @@ -457,7 +458,7 @@ pub fn comp_disc_edge( phase: f64, expose: f64, n_div: i32, - disc_edge_grid: &Vec, + disc_edge_grid: &Grid, ) -> f64 { let ri: f64 = iangle.to_radians(); let (sini, cosi) = ri.sin_cos(); @@ -488,7 +489,7 @@ pub fn comp_disc_edge( ssum = 0.0; let phi_normed: f64 = phi - phi.floor(); // Disc edge - for point in disc_edge_grid { + for point in &disc_edge_grid.points { mu = earth.dot(&point.direction); if mu > 0.0 && point.is_visible_phase_normed(phi_normed) { ssum += mu @@ -508,7 +509,7 @@ pub fn comp_bright_spot( phase: f64, expose: f64, n_div: i32, - bright_spot_grid: &Vec, + bright_spot_grid: &Grid, ) -> f64 { let ri = iangle.to_radians(); let (sini, cosi) = ri.sin_cos(); @@ -538,7 +539,7 @@ pub fn comp_bright_spot( ssum = 0.0; // Bright spot - for point in bright_spot_grid { + for point in &bright_spot_grid.points { mu = earth.dot(&point.direction); if mu > 0.0 && point.is_visible(phi) { ssum += mu * (point.flux as f64); diff --git a/src/grid.rs b/src/grid.rs index 70f5b75..1d66280 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -2,30 +2,51 @@ use roche::{self, Point, Vec3}; use std::f64::consts::TAU; use pyo3::prelude::*; +/// +/// Grid is a struct to hold the vector of `Vec` defining a +/// component grid along with methods to act on this. +/// #[pyclass(skip_from_py_object)] #[derive(Clone, Debug)] pub struct Grid { #[pyo3(get)] pub points: Vec, - #[pyo3(get)] pub q: f64, - #[pyo3(get)] pub iangle: f64, } +impl Grid { + + pub fn new(points: Vec, q: f64, iangle: f64) -> Self { + Self { + points, + q, + iangle, + } + } +} + #[pymethods] impl Grid { - #[pyo3(signature = (phase, iangle=None))] - pub fn project_2d(&self, phase: f64, iangle: Option) -> (Vec, Vec) { - let iangle = match iangle { - Some(iangle) => iangle, - None => self.iangle, - }; + + /// + /// Projects the grid onto a 2D plane as seen at the model inclination + /// at the supplied phase with the binary centre of mass as the origin. + /// + /// Arguments + /// * `phase` - Orbital phase at which to project the grid + /// + /// Returns + /// (x, y) - Arrays of projected grid point positions + /// + #[pyo3(signature = (phase))] + pub fn project_2d(&self, phase: f64) -> (Vec, Vec) { + let mut x_arr: Vec = vec![]; let mut y_arr: Vec = vec![]; let cofm = Vec3::new(self.q/(1.0+self.q), 0.0, 0.0); let (sinp, cosp) = (TAU*phase).sin_cos(); - let earth = roche::set_earth_iangle(iangle, phase); + let earth = roche::set_earth_iangle(self.iangle, phase); let xsky = Vec3::new(sinp, cosp, 0.0); let ysky = earth.cross(&xsky); diff --git a/src/set_bright_spot_grid.rs b/src/set_bright_spot_grid.rs index b28855a..8b5d1c4 100644 --- a/src/set_bright_spot_grid.rs +++ b/src/set_bright_spot_grid.rs @@ -1,3 +1,4 @@ +use crate::grid::Grid; use crate::model::Model; use crate::set_star_grid::star_eclipse; use roche::{self, Etype, Point, RocheContext, Star, Vec3, errors::RocheError}; @@ -22,7 +23,7 @@ use std::panic; // \exception Exceptions are thrown if the specified radii over-fill the Roche lobes. // -pub fn set_bright_spot_grid(model: &Model) -> Result, RocheError> { +pub fn set_bright_spot_grid(model: &Model) -> Result { let (mut r1, mut r2) = model.get_r1r2(); let roche_context1 = RocheContext::new(model.q.value, Star::Primary, model.spin1.value)?; @@ -152,5 +153,5 @@ pub fn set_bright_spot_grid(model: &Model) -> Result, RocheError> { flux: flux_parallel, }; } - Ok(bright_spot_grid) + Ok(Grid::new(bright_spot_grid, model.q.value, model.iangle.value)) } diff --git a/src/set_disc_continuum.rs b/src/set_disc_continuum.rs index 2d0813e..8aac361 100644 --- a/src/set_disc_continuum.rs +++ b/src/set_disc_continuum.rs @@ -1,6 +1,8 @@ -use roche::{self, Point, Vec3}; +use roche::{self, Vec3}; use std::f64::consts::{FRAC_PI_2, PI}; +use crate::grid::Grid; + // set_disc_continuum computes the face-on brightness of each element of the // disc assuming a power law with radius. // @@ -13,11 +15,11 @@ use std::f64::consts::{FRAC_PI_2, PI}; // \param wave wavelength of interest, nm // \param disc grid of elements over disc -pub fn set_disc_continuum(rdisc: f64, tdisc: f64, texp: f64, wave: f64, disc: &mut Vec) { +pub fn set_disc_continuum(rdisc: f64, tdisc: f64, texp: f64, wave: f64, disc: &mut Grid) { // Reference surface brightness let bright: f64 = roche::planck(wave, tdisc); - for point in disc { + for point in &mut disc.points { let r: f64 = point.position.length(); point.flux = (bright * (r / rdisc).powf(texp) * point.area as f64) as f32; } @@ -43,7 +45,7 @@ pub fn set_edge_continuum( t2: f64, absorb: f64, wave: f64, - edge: &mut Vec, + edge: &mut Grid, ) { let mut vec: Vec3; let cofm2: Vec3 = Vec3::cofm2(); @@ -52,7 +54,7 @@ pub fn set_edge_continuum( let mut mu: f64; let mut r: f64; - for point in edge { + for point in &mut edge.points { vec = cofm2 - point.position; r = vec.length(); mu = vec.dot(&point.direction) / r; diff --git a/src/set_disc_grid.rs b/src/set_disc_grid.rs index 1137736..fbd3e57 100644 --- a/src/set_disc_grid.rs +++ b/src/set_disc_grid.rs @@ -1,3 +1,4 @@ +use crate::grid::Grid; use crate::model::Model; use crate::set_star_grid::star_eclipse; use rayon::iter::{IntoParallelIterator, ParallelIterator}; @@ -25,7 +26,7 @@ use std::panic; /// \exception Exceptions are thrown if the specified radii over-fill the /// Roche lobes. /// -pub fn set_disc_grid(model: &Model) -> Result, RocheError> { +pub fn set_disc_grid(model: &Model) -> Result { const EFAC: f64 = 1.0000001; let (mut r1, mut r2) = model.get_r1r2(); @@ -83,7 +84,7 @@ pub fn set_disc_grid(model: &Model) -> Result, RocheError> { let drad: f64 = (rdisc2 - rdisc1) / model.nrad as f64; let drrad: f64 = rdisc2 / model.nrad as f64; - let disc_grid: Vec = (0..model.nrad) + let disc_grid_points: Vec = (0..model.nrad) .into_par_iter() .flat_map_iter(|i| { let rad: f64 = rdisc1 + (rdisc2 - rdisc1) * (i as f64 + 0.5) / model.nrad as f64; @@ -150,7 +151,7 @@ pub fn set_disc_grid(model: &Model) -> Result, RocheError> { }) }) .collect(); - Ok(disc_grid) + Ok(Grid::new(disc_grid_points, model.q.value, model.iangle.value)) } /// @@ -177,7 +178,7 @@ pub fn set_disc_edge_grid( model: &Model, outer: bool, visual: bool, -) -> Result, RocheError> { +) -> Result { const EFAC: f64 = 1.0000001; let (mut r1, mut r2) = model.get_r1r2(); @@ -374,5 +375,5 @@ pub fn set_disc_edge_grid( } } - Ok(edge_grid) + Ok(Grid::new(edge_grid, model.q.value, model.iangle.value)) } diff --git a/src/set_star_continuum.rs b/src/set_star_continuum.rs index 4b3688f..b7b3799 100644 --- a/src/set_star_continuum.rs +++ b/src/set_star_continuum.rs @@ -1,11 +1,11 @@ -use crate::model::Model; -use roche::{self, Point, Vec3, constants::EFAC, errors::RocheError}; +use crate::{grid::Grid, model::Model}; +use roche::{self, Vec3, constants::EFAC, errors::RocheError}; use std::f64::consts::PI; pub fn set_star_continuum( model: &Model, - star1: &mut Vec, - star2: &mut Vec, + star1: &mut Grid, + star2: &mut Grid, ) -> Result<(), RocheError> { let (mut r1, mut r2) = model.get_r1r2(); @@ -90,7 +90,7 @@ pub fn set_star_continuum( Vec3::new(0.0, 0.0, 0.0) }; - for point in star1 { + for point in &mut star1.points { let vec: Vec3 = cofm2 - point.position; let r: f64 = vec.length(); let mu: f64 = point.direction.dot(&vec) / r; @@ -226,7 +226,7 @@ pub fn set_star_continuum( Vec3::new(0.0, 0.0, 0.0) }; - for point in star2 { + for point in &mut star2.points { let vec: Vec3 = cofm1 - point.position; let r: f64 = vec.length(); let mu: f64 = point.direction.dot(&vec) / r; diff --git a/src/set_star_grid.rs b/src/set_star_grid.rs index 1a765cd..db2cef4 100644 --- a/src/set_star_grid.rs +++ b/src/set_star_grid.rs @@ -1,3 +1,4 @@ +use crate::grid::Grid; use crate::model::Model; use crate::numface::numface; use rayon::prelude::*; @@ -22,7 +23,7 @@ pub fn envelope(rangle: f64, lambda: f64, r1: f64) -> Xy { } } -pub fn set_star_grid(model: &Model, star: Star, fine: bool) -> Result, RocheError> { +pub fn set_star_grid(model: &Model, star: Star, fine: bool) -> Result { let (mut r1, mut r2) = model.get_r1r2(); let eclipse: bool = match star { @@ -166,7 +167,11 @@ pub fn set_star_grid(model: &Model, star: Star, fine: bool) -> Result let nface: u32 = numface(nlat, infill, thelo, thehi, nlatfill, nlngfill); // Generate arrays over the star's face - let mut star_grid: Vec = Vec::with_capacity(nface as usize); + let mut star_grid: Grid = Grid::new( + Vec::with_capacity(nface as usize), + model.q.value, + model.iangle.value + ); let acc: f64 = model.delta_phase / 10.0; @@ -308,7 +313,7 @@ pub fn set_star_grid(model: &Model, star: Star, fine: bool) -> Result } pub fn add_faces( - star_grid: &mut Vec, + star_grid: &mut Grid, tlo: f64, thi: f64, dtheta: f64, @@ -486,10 +491,10 @@ pub fn add_faces( }) .collect(); - star_grid.clear(); + star_grid.points.clear(); for band in bands { - star_grid.extend(band); + star_grid.points.extend(band); } } From 3e253cc3c7361c48c1c56e61accb5d9749e907a5 Mon Sep 17 00:00:00 2001 From: Alex-J-Brown Date: Fri, 21 Aug 2026 13:46:21 +0200 Subject: [PATCH 3/4] add more access methods to Grid --- src/grid.rs | 99 ++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 94 insertions(+), 5 deletions(-) diff --git a/src/grid.rs b/src/grid.rs index a6d3cd6..046ce05 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -1,5 +1,6 @@ use roche::{self, Point, Vec3}; use std::f64::consts::TAU; +use numpy::{IntoPyArray, PyArray1}; use pyo3::prelude::*; /// @@ -24,12 +25,82 @@ impl Grid { iangle, } } -} -#[pymethods] -impl Grid { + pub fn position(&self, phase: Option) -> Vec { + + let position: Vec = match phase { + Some(phase) => { + let mut position: Vec = vec![]; + let earth = roche::set_earth_iangle(self.iangle, phase); + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + position.push(point.position); + } + } + position + }, + None => { + let mut position: Vec = vec![]; + for point in &self.points { + position.push(point.position); + } + position + } + }; + + position + } + + pub fn direction(&self, phase: Option) -> Vec { + + let direction: Vec = match phase { + Some(phase) => { + let mut direction: Vec = vec![]; + let earth = roche::set_earth_iangle(self.iangle, phase); + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + direction.push(point.direction); + } + } + direction + }, + None => { + let mut direction: Vec = vec![]; + for point in &self.points { + direction.push(point.direction); + } + direction + } + }; + + direction + } + + pub fn gravity(&self, phase: Option) -> Vec { + + let gravity: Vec = match phase { + Some(phase) => { + let mut gravity: Vec = vec![]; + let earth = roche::set_earth_iangle(self.iangle, phase); + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + gravity.push(point.gravity); + } + } + gravity + }, + None => { + let mut gravity: Vec = vec![]; + for point in &self.points { + gravity.push(point.gravity); + } + gravity + } + }; + + gravity + } - #[pyo3(signature = (phase=None))] pub fn area(&self, phase: Option) -> Vec { let area: Vec = match phase { @@ -55,7 +126,6 @@ impl Grid { area } - #[pyo3(signature = (phase=None))] pub fn flux(&self, phase: Option) -> Vec { let flux: Vec = match phase { @@ -81,6 +151,25 @@ impl Grid { flux } +} + +#[pymethods] +impl Grid { + + #[pyo3(name="area", signature = (phase=None))] + pub fn python_area(&self, py: Python, phase: Option) -> Py> { + + let area: Vec = self.area(phase); + area.into_pyarray(py).unbind() + } + + #[pyo3(name="flux", signature = (phase=None))] + pub fn python_flux(&self, py: Python, phase: Option) -> Py> { + + let flux: Vec = self.flux(phase); + flux.into_pyarray(py).unbind() + } + /// /// Projects the grid onto a 2D plane as seen at the model inclination /// at the supplied phase with the binary centre of mass as the origin. From e13b23c9c502c1545eb5c44dd4b39ccb296393b5 Mon Sep 17 00:00:00 2001 From: Alex-J-Brown Date: Fri, 4 Sep 2026 16:59:29 +0200 Subject: [PATCH 4/4] remove q and iangle as Grid attributes --- src/binary_model.rs | 6 ++--- src/grid.rs | 48 ++++++++++++++++++------------------- src/set_bright_spot_grid.rs | 2 +- src/set_disc_grid.rs | 4 ++-- src/set_star_grid.rs | 4 +--- 5 files changed, 30 insertions(+), 34 deletions(-) diff --git a/src/binary_model.rs b/src/binary_model.rs index a9e1afe..5391bab 100644 --- a/src/binary_model.rs +++ b/src/binary_model.rs @@ -628,9 +628,9 @@ fn build_grids( set_star_continuum(model, &mut star1_coarse_grid, &mut star2_coarse_grid)?; } - let mut disc_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); - let mut disc_edge_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); - let mut bright_spot_grid: Grid = Grid::new(vec![], model.q.value, model.iangle.value); + let mut disc_grid: Grid = Grid::new(vec![]); + let mut disc_edge_grid: Grid = Grid::new(vec![]); + let mut bright_spot_grid: Grid = Grid::new(vec![]); let mut rlens1 = 0.0; if model.glens1 { diff --git a/src/grid.rs b/src/grid.rs index 046ce05..44d1613 100644 --- a/src/grid.rs +++ b/src/grid.rs @@ -12,26 +12,22 @@ use pyo3::prelude::*; pub struct Grid { #[pyo3(get)] pub points: Vec, - pub q: f64, - pub iangle: f64, } impl Grid { - pub fn new(points: Vec, q: f64, iangle: f64) -> Self { + pub fn new(points: Vec) -> Self { Self { points, - q, - iangle, } } - pub fn position(&self, phase: Option) -> Vec { + pub fn position(&self, iangle: f64, phase: Option) -> Vec { let position: Vec = match phase { Some(phase) => { let mut position: Vec = vec![]; - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); for point in &self.points { if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { position.push(point.position); @@ -51,12 +47,12 @@ impl Grid { position } - pub fn direction(&self, phase: Option) -> Vec { + pub fn direction(&self, iangle: f64, phase: Option) -> Vec { let direction: Vec = match phase { Some(phase) => { let mut direction: Vec = vec![]; - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); for point in &self.points { if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { direction.push(point.direction); @@ -76,12 +72,12 @@ impl Grid { direction } - pub fn gravity(&self, phase: Option) -> Vec { + pub fn gravity(&self, iangle: f64, phase: Option) -> Vec { let gravity: Vec = match phase { Some(phase) => { let mut gravity: Vec = vec![]; - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); for point in &self.points { if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { gravity.push(point.gravity); @@ -101,12 +97,12 @@ impl Grid { gravity } - pub fn area(&self, phase: Option) -> Vec { + pub fn area(&self, iangle: f64, phase: Option) -> Vec { let area: Vec = match phase { Some(phase) => { let mut area: Vec = vec![]; - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); for point in &self.points { if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { area.push(point.area); @@ -126,12 +122,12 @@ impl Grid { area } - pub fn flux(&self, phase: Option) -> Vec { + pub fn flux(&self, iangle: f64, phase: Option) -> Vec { let flux: Vec = match phase { Some(phase) => { let mut flux: Vec = vec![]; - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); for point in &self.points { if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { flux.push(point.flux); @@ -156,17 +152,17 @@ impl Grid { #[pymethods] impl Grid { - #[pyo3(name="area", signature = (phase=None))] - pub fn python_area(&self, py: Python, phase: Option) -> Py> { + #[pyo3(name="area", signature = (iangle, phase=None))] + pub fn python_area(&self, py: Python, iangle: f64, phase: Option) -> Py> { - let area: Vec = self.area(phase); + let area: Vec = self.area(iangle, phase); area.into_pyarray(py).unbind() } - #[pyo3(name="flux", signature = (phase=None))] - pub fn python_flux(&self, py: Python, phase: Option) -> Py> { + #[pyo3(name="flux", signature = (iangle, phase=None))] + pub fn python_flux(&self, py: Python, iangle: f64, phase: Option) -> Py> { - let flux: Vec = self.flux(phase); + let flux: Vec = self.flux(iangle, phase); flux.into_pyarray(py).unbind() } @@ -175,19 +171,21 @@ impl Grid { /// at the supplied phase with the binary centre of mass as the origin. /// /// Arguments + /// * `q` - Binary mass ratio M2/M1. + /// * `iangle` - Orbital inclination at which to project the grid. /// * `phase` - Orbital phase at which to project the grid /// /// Returns /// (x, y) - Arrays of projected grid point positions /// - #[pyo3(signature = (phase))] - pub fn project_2d(&self, phase: f64) -> (Vec, Vec) { + #[pyo3(signature = (q, iangle, phase))] + pub fn project_2d(&self, q: f64, iangle: f64, phase: f64) -> (Vec, Vec) { let mut x_arr: Vec = vec![]; let mut y_arr: Vec = vec![]; - let cofm = Vec3::new(self.q/(1.0+self.q), 0.0, 0.0); + let cofm = Vec3::new(q/(1.0+q), 0.0, 0.0); let (sinp, cosp) = (TAU*phase).sin_cos(); - let earth = roche::set_earth_iangle(self.iangle, phase); + let earth = roche::set_earth_iangle(iangle, phase); let xsky = Vec3::new(sinp, cosp, 0.0); let ysky = earth.cross(&xsky); diff --git a/src/set_bright_spot_grid.rs b/src/set_bright_spot_grid.rs index 8b5d1c4..0274847 100644 --- a/src/set_bright_spot_grid.rs +++ b/src/set_bright_spot_grid.rs @@ -153,5 +153,5 @@ pub fn set_bright_spot_grid(model: &Model) -> Result { flux: flux_parallel, }; } - Ok(Grid::new(bright_spot_grid, model.q.value, model.iangle.value)) + Ok(Grid::new(bright_spot_grid)) } diff --git a/src/set_disc_grid.rs b/src/set_disc_grid.rs index fbd3e57..02502f7 100644 --- a/src/set_disc_grid.rs +++ b/src/set_disc_grid.rs @@ -151,7 +151,7 @@ pub fn set_disc_grid(model: &Model) -> Result { }) }) .collect(); - Ok(Grid::new(disc_grid_points, model.q.value, model.iangle.value)) + Ok(Grid::new(disc_grid_points)) } /// @@ -375,5 +375,5 @@ pub fn set_disc_edge_grid( } } - Ok(Grid::new(edge_grid, model.q.value, model.iangle.value)) + Ok(Grid::new(edge_grid)) } diff --git a/src/set_star_grid.rs b/src/set_star_grid.rs index e11d0a1..605e4ea 100644 --- a/src/set_star_grid.rs +++ b/src/set_star_grid.rs @@ -168,9 +168,7 @@ pub fn set_star_grid(model: &Model, star: Star, fine: bool) -> Result