1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
use super::{binning, shift_and_add};
use crate::{Field, Intensity, Observer, Observing, ObservingModes};
use indicatif::ProgressBar;
use rand_distr::{Distribution, Poisson};
impl<T, Mode> Field<T, Mode>
where
T: Observer + Sync + Send,
Mode: ObservingModes + Send,
Observing<Mode>: Intensity,
{
/// Computes field-of-view intensity map
pub fn intensity(&mut self, bar: Option<ProgressBar>) -> Vec<f64> {
// Telescope Nyquist-Shannon sampling criteria
// let nyquist = 0.5 * self.photometry.wavelength / self.observer.diameter();
// Image resolution to sampling criteria ratio
let b = self
.pixel_scale
.to_nyquist_clamped_ratio(&self.observer, &self.photometry);
// Intensity sampling (oversampled wrt. image by factor b>=1)
let intensity_sampling = (b * self.field_of_view.to_pixelscale_ratio(self)).ceil() as usize;
// Pupil size according to intensity angular resolution
let pupil_size = b * self.photometry.wavelength / self.resolution();
// FFT sampling based on pupil spatial resolution
let mut n_dft = (pupil_size / self.observer.resolution()).ceil() as usize;
// Match parity of FFT and intensity sampling if the latter is larger
if intensity_sampling > n_dft && intensity_sampling % 2 != n_dft % 2 {
n_dft += 1;
}
log::debug!(
r"
. Image sampling: {intensity_sampling}:{b}
. Pupil size : {pupil_size:.3}m
. DFT sampling : {n_dft}
"
);
// Zero-padding discrete Fourier transform
self.observing_mode
.init_fft(n_dft, self.observer.resolution());
// star image stacking buffer
let mut buffer = vec![0f64; intensity_sampling.pow(2)];
let n = intensity_sampling as i32;
let alpha = self.resolution() / b;
let n_threads = match std::env::var("N_THREAD") {
Ok(n) => n.parse::<usize>().unwrap(),
Err(__) => num_cpus::get(),
};
log::info!("Computing intensity in parallel using {n_threads} threads");
for star_chunk in self.objects.chunks(n_threads) {
let intensities: Vec<_> = std::thread::scope(|s| {
log::info!("starting intensity batch");
let mut intensity_scope = vec![];
for star in star_chunk {
bar.as_ref().map(|b| b.inc(1));
// todo: check if star is within FOV (rejection criteria?)
if !star.inside_box(self.field_of_view() + self.resolution() * 2.) {
continue;
}
let n_photon = self.flux.unwrap_or(
self.photometry.n_photon(star.magnitude)
* self.exposure
* self.observer.resolution().powi(2), // * self.observer.area() ,
);
// star coordinates
let (x, y) = star.coordinates;
// integer part
let x0 = -(y / alpha).round();
let y0 = (x / alpha).round();
// fractional part
let fr_x0 = -y.to_radians() - x0 * alpha;
let fr_y0 = x.to_radians() - y0 * alpha;
// image fractional translation by Fourier interpolation
let shift = if intensity_sampling % 2 == 0 {
Some((
0.5 / pupil_size + fr_x0 / self.photometry.wavelength,
0.5 / pupil_size + fr_y0 / self.photometry.wavelength,
))
} else {
Some((
fr_x0 / self.photometry.wavelength,
fr_y0 / self.photometry.wavelength,
))
};
let mut observing_mode = self.observing_mode.clone();
let pupil_resolution = self.observer.resolution();
let mut pupil = self.observer.pupil(shift);
pupil.iter_mut().for_each(|p| *p *= n_photon.sqrt());
let poisson_noise = self.poisson_noise;
let intensity_thread = s.spawn(move || {
let mut rng = rand::thread_rng();
// star intensity map
// Zero-padding discrete Fourier transform
observing_mode.init_fft(n_dft, pupil_resolution);
let mut intensity = observing_mode
.intensity(pupil, intensity_sampling, star)
.unwrap();
// intensity set to # of photon & Poisson noise
// log::debug!("Image flux: {n_photon}");
if poisson_noise {
intensity.iter_mut().for_each(|i| {
if *i == 0f64 {
*i = 0f64;
} else {
let poi = Poisson::new(*i).unwrap();
*i = poi.sample(&mut rng)
}
})
};
(x0, y0, intensity)
});
intensity_scope.push(intensity_thread);
}
intensity_scope
.into_iter()
.map(|intensity_thread| intensity_thread.join().unwrap())
.collect()
});
log::info!("star images shift & add");
for (x0, y0, intensity) in intensities.into_iter() {
// shift and add star images
shift_and_add(buffer.as_mut_slice(), x0, y0, n, intensity);
}
}
bar.as_ref().map(|b| b.finish());
let m = b as usize;
if m == 1 {
return buffer;
}
// binning
binning(intensity_sampling, m, buffer)
}
}