diff --git a/src/binary_model.rs b/src/binary_model.rs index e4f4bb9..5391bab 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; @@ -15,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; @@ -79,19 +80,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, @@ -215,31 +216,31 @@ impl BinaryModel { /// "disc", /// "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(()) + // } /// /// Computes a model light curve for an array of times and exposure times @@ -409,11 +410,11 @@ impl BinaryModel { }; 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(), @@ -578,13 +579,13 @@ fn build_grids( model: &Model, ) -> Result< ( - Vec, - Vec, - Vec, - Vec, - Vec, - Vec, - Vec, + Grid, + Grid, + Grid, + Grid, + Grid, + Grid, + Grid, Ginterp, f64, bool, @@ -595,8 +596,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)?; @@ -627,9 +628,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![]); + 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 { @@ -668,7 +669,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, @@ -682,7 +683,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, @@ -696,7 +697,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, @@ -710,7 +711,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, @@ -768,10 +769,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)?; @@ -1019,14 +1020,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/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 new file mode 100644 index 0000000..44d1613 --- /dev/null +++ b/src/grid.rs @@ -0,0 +1,202 @@ +use roche::{self, Point, Vec3}; +use std::f64::consts::TAU; +use numpy::{IntoPyArray, PyArray1}; +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, +} + +impl Grid { + + pub fn new(points: Vec) -> Self { + Self { + points, + } + } + + 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(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, iangle: f64, phase: Option) -> Vec { + + let direction: Vec = match phase { + Some(phase) => { + let mut direction: Vec = vec![]; + 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); + } + } + direction + }, + None => { + let mut direction: Vec = vec![]; + for point in &self.points { + direction.push(point.direction); + } + direction + } + }; + + direction + } + + 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(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 + } + + 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(iangle, phase); + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + area.push(point.area); + } + } + area + }, + None => { + let mut area: Vec = vec![]; + for point in &self.points { + area.push(point.area); + } + area + } + }; + + area + } + + 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(iangle, phase); + for point in &self.points { + if earth.dot(&point.direction) > 0.0 && point.is_visible(phase) { + flux.push(point.flux); + } + } + flux + }, + None => { + let mut flux: Vec = vec![]; + for point in &self.points { + flux.push(point.flux); + } + flux + } + }; + + flux + } + +} + +#[pymethods] +impl Grid { + + #[pyo3(name="area", signature = (iangle, phase=None))] + pub fn python_area(&self, py: Python, iangle: f64, phase: Option) -> Py> { + + let area: Vec = self.area(iangle, phase); + area.into_pyarray(py).unbind() + } + + #[pyo3(name="flux", signature = (iangle, phase=None))] + pub fn python_flux(&self, py: Python, iangle: f64, phase: Option) -> Py> { + + let flux: Vec = self.flux(iangle, 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. + /// + /// 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 = (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(q/(1.0+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 862113a..24640a4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -5,6 +5,7 @@ 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; @@ -29,6 +30,9 @@ mod lcurve { #[pymodule_export] use crate::ldc::LDCType; + + #[pymodule_export] + use crate::grid::Grid; #[allow(non_upper_case_globals)] #[pymodule_export] diff --git a/src/set_bright_spot_grid.rs b/src/set_bright_spot_grid.rs index b28855a..0274847 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)) } 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..02502f7 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)) } /// @@ -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)) } 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 e2bff2e..605e4ea 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,9 @@ 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) + ); let acc: f64 = model.delta_phase / 10.0; @@ -308,7 +311,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, @@ -488,7 +491,7 @@ pub fn add_faces( for band in bands { - star_grid.extend(band); + star_grid.points.extend(band); } }