materials_for_mc 0.2.0

Material builder for Monte Carlo simulations
Documentation 
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

// ...existing code...
/// Utility functions for materials_for_mc

#[doc(hidden)]
/// Linear interpolation on a linear scale (hidden from docs)
///
/// Given arrays of x and y values, interpolate to find the y value at x_new.
/// If x_new is outside the range of x, returns the first or last y value.
pub fn interpolate_linear(x: &[f64], y: &[f64], x_new: f64) -> f64 {
    // Edge cases
    if x.is_empty() {
        return f64::NAN;
    }
    if x.len() == 1 {
        return y[0];
    }
    if x_new <= x[0] {
        return y[0];
    }
    if x_new >= x[x.len() - 1] {
        return y[y.len() - 1];
    }

    // Binary search for interval: find largest i with x[i] <= x_new
    let mut low = 0usize;
    let mut high = x.len() - 1; // invariant: target interval within (low, high]
    while high - low > 1 {
        // maintain at least two-point span
        let mid = (low + high) >> 1; // divide by 2
        if x[mid] <= x_new {
            low = mid;
        } else {
            high = mid;
        }
    }
    let idx = low; // x[idx] <= x_new < x[idx+1]
    let x1 = x[idx];
    let x2 = x[idx + 1];
    let y1 = y[idx];
    let y2 = y[idx + 1];
    y1 + (x_new - x1) * (y2 - y1) / (x2 - x1)
}

#[doc(hidden)]
/// Log-log interpolation (hidden from docs)
///
/// Given arrays of x and y values, interpolate on a log-log scale to find the y value at x_new.
/// If x_new is outside the range of x, returns the first or last y value.
/// All x and y values must be positive.
pub fn interpolate_log_log(x: &[f64], y: &[f64], x_new: f64) -> f64 {
    // Edge cases / validation
    if x.is_empty() {
        return f64::NAN;
    }
    if x.len() == 1 {
        return y[0];
    }
    if x_new <= x[0] {
        return y[0];
    }
    if x_new >= x[x.len() - 1] {
        return y[y.len() - 1];
    }

    // Binary search for interval
    let mut low = 0usize;
    let mut high = x.len() - 1;
    while high - low > 1 {
        let mid = (low + high) >> 1;
        if x[mid] <= x_new {
            low = mid;
        } else {
            high = mid;
        }
    }
    let idx = low;
    let x1 = x[idx];
    let x2 = x[idx + 1];
    let y1 = y[idx];
    let y2 = y[idx + 1];
    let log_x1 = x1.ln();
    let log_x2 = x2.ln();
    let log_y1 = y1.ln();
    let log_y2 = y2.ln();
    let log_x_new = x_new.ln();
    let log_y_new = log_y1 + (log_x_new - log_x1) * (log_y2 - log_y1) / (log_x2 - log_x1);
    log_y_new.exp()
}