accrete/structs/
dust.rs

1use crate::consts::*;
2use crate::utils::*;
3
4use serde::{Deserialize, Serialize};
5
6pub type DustBands = Vec<DustBand>;
7#[derive(Debug, Copy, Clone, Serialize, Deserialize, PartialEq)]
8pub struct DustBand {
9    pub outer_edge: f64,
10    pub inner_edge: f64,
11    pub dust_present: bool,
12    pub gas_present: bool,
13}
14
15impl DustBand {
16    pub fn new(outer_edge: f64, inner_edge: f64, dust_present: bool, gas_present: bool) -> Self {
17        Self {
18            outer_edge,
19            inner_edge,
20            dust_present,
21            gas_present,
22        }
23    }
24}
25
26pub fn dust_availible(dust_bands: &[DustBand], inside_range: &f64, outside_range: &f64) -> bool {
27    dust_bands.iter().rev().fold(false, |mut acc, band| {
28        if band.dust_present && &band.outer_edge > inside_range && &band.inner_edge < outside_range
29        {
30            acc = true;
31        }
32        acc
33    })
34}
35
36/// "The center of mass is occupied by a star with a mass of one unit (one solar mass). All particles in the cloud are moving on elliptical orbits, with the center of mass at one focus. The density of dust (p1) within the cloud depends on a function of the form p1 = A exp (-arl/n). The overall density of gas and dust (p2) within the cloud equals Kpl, where r is distance from the center of mass (in astronomical units, a.u.) and A. a. n. and K (the vas/dust ratio) are constants."
37/// "There is a spherically symmetrical cloud of dust and gas with a constant ratio of gas to dust, the density decreasing with distance from the center."
38pub fn accrete_dust(
39    mass: &mut f64,
40    a: &f64,
41    e: &f64,
42    crit_mass: &f64,
43    dust_bands: &mut [DustBand],
44    cloud_eccentricity: &f64,
45    dust_density: &f64,
46    k: &f64,
47) {
48    let mut new_mass = *mass;
49    loop {
50        *mass = new_mass;
51        new_mass = 0.0;
52
53        for d in dust_bands.iter_mut() {
54            new_mass += collect_dust(
55                mass,
56                a,
57                e,
58                crit_mass,
59                cloud_eccentricity,
60                dust_density,
61                k,
62                d,
63            );
64        }
65
66        if !((new_mass - *mass) >= (0.0001 * *mass)) {
67            break;
68        }
69    }
70    *mass = new_mass;
71}
72
73pub fn collect_dust(
74    mass: &f64,
75    a: &f64,
76    e: &f64,
77    crit_mass: &f64,
78    cloud_eccentricity: &f64,
79    dust_density: &f64,
80    k: &f64,
81    band: &mut DustBand,
82) -> f64 {
83    let mut r_inner = inner_swept_limit(a, e, mass, cloud_eccentricity);
84    let r_outer = outer_swept_limit(a, e, mass, cloud_eccentricity);
85
86    if r_inner < 0.0 {
87        r_inner = 0.0;
88    }
89
90    if band.outer_edge <= r_inner || band.inner_edge >= r_outer || !band.dust_present {
91        return 0.0;
92    };
93
94    let density = match !band.gas_present || mass < crit_mass {
95        true => *dust_density,
96        false => get_mass_density(k, dust_density, crit_mass, mass),
97    };
98    let bandwidth = r_outer - r_inner;
99    let temp1 = match r_outer - band.outer_edge > 0.0 {
100        true => r_outer - band.outer_edge,
101        false => 0.0,
102    };
103    let temp2 = match band.inner_edge - r_inner > 0.0 {
104        true => band.inner_edge - r_inner,
105        false => 0.0,
106    };
107    let width = bandwidth - temp1 - temp2;
108    let term1 = 4.0 * PI * a.powf(2.0);
109    let term2 = 1.0 - e * (temp1 - temp2) / bandwidth;
110    let volume = term1 * reduced_mass(mass) * width * term2;
111
112    volume * density
113}
114
115pub fn update_dust_lanes(
116    dust_bands: &mut Vec<DustBand>,
117    min: f64,
118    max: f64,
119    mass: &f64,
120    crit_mass: &f64,
121) {
122    let mut gas = true;
123    if mass > crit_mass {
124        gas = false;
125    }
126    *dust_bands = dust_bands.iter_mut().fold(Vec::new(), |mut acc, band| {
127        let new_gas = band.gas_present && gas;
128        if band.inner_edge < min && band.outer_edge > max {
129            let mut inner = *band;
130            inner.outer_edge = min;
131            let middle = DustBand::new(max, min, false, new_gas);
132            let outer = DustBand::new(band.outer_edge, max, band.dust_present, band.gas_present);
133            acc.push(inner);
134            acc.push(middle);
135            acc.push(outer);
136        } else if band.inner_edge < max && band.outer_edge > max {
137            let outer = DustBand::new(band.outer_edge, max, band.dust_present, band.gas_present);
138            let inner = DustBand::new(max, band.inner_edge, false, new_gas);
139            acc.push(inner);
140            acc.push(outer);
141        } else if band.inner_edge < min && band.outer_edge > min {
142            let outer = DustBand::new(band.outer_edge, min, false, new_gas);
143            let inner = DustBand::new(min, band.inner_edge, band.dust_present, band.gas_present);
144            acc.push(inner);
145            acc.push(outer);
146        } else if band.inner_edge >= min && band.outer_edge <= max {
147            let dust_band = DustBand::new(band.outer_edge, band.inner_edge, false, new_gas);
148            acc.push(dust_band)
149        } else if band.outer_edge < min || band.inner_edge > max {
150            acc.push(*band);
151        }
152        acc
153    });
154}
155
156pub fn compress_dust_lanes(dust_bands: &mut Vec<DustBand>) {
157    *dust_bands = dust_bands
158        .iter()
159        .enumerate()
160        .fold(Vec::new(), |mut acc, (i, band)| {
161            match dust_bands.get(i + 1) {
162                Some(next_band) => {
163                    if band.dust_present == next_band.dust_present
164                        && band.gas_present == next_band.gas_present
165                    {
166                        let mut band = *band;
167                        band.outer_edge = next_band.outer_edge;
168                        acc.push(band);
169                    } else {
170                        acc.push(*band);
171                    }
172                }
173                None => acc.push(*band),
174            }
175            acc
176        });
177}
178
179pub fn get_mass_density(k: &f64, dust_density: &f64, critical_mass: &f64, mass: &f64) -> f64 {
180    k * dust_density / (1.0 + (critical_mass / mass).sqrt() * (k - 1.0))
181}
182
183pub fn dust_density(dust_density_coeff: &f64, stellar_mass: &f64, oribital_radius: &f64) -> f64 {
184    dust_density_coeff * stellar_mass.sqrt() * (-ALPHA * oribital_radius.powf(1.0 / N)).exp()
185}
accrete/structs/dust.rs

accrete/structs/
dust.rs
