sciforge_lib/physics/quantum/
wavefunctions.rs

1use crate::maths::complex::Complex;
2
3pub fn plane_wave(x: f64, k: f64, omega: f64, t: f64) -> Complex {
4    Complex::from_polar(1.0, k * x - omega * t)
5}
6
7pub fn gaussian_packet(x: f64, x0: f64, sigma: f64, k0: f64) -> Complex {
8    let norm = 1.0 / (sigma * (2.0 * std::f64::consts::PI).sqrt()).sqrt();
9    let gauss = (-(x - x0).powi(2) / (4.0 * sigma * sigma)).exp();
10    Complex::from_polar(norm * gauss, k0 * x)
11}
12
13pub fn normalize(psi: &mut [Complex], dx: f64) {
14    let norm_sq: f64 = psi.iter().map(|c| c.norm() * c.norm()).sum::<f64>() * dx;
15    let inv_norm = 1.0 / norm_sq.sqrt();
16    for c in psi.iter_mut() {
17        *c = *c * Complex::new(inv_norm, 0.0);
18    }
19}
20
21pub fn probability_density(psi: &[Complex]) -> Vec<f64> {
22    psi.iter().map(|c| c.norm() * c.norm()).collect()
23}
24
25pub fn expectation_position(psi: &[Complex], x: &[f64], dx: f64) -> f64 {
26    psi.iter()
27        .zip(x.iter())
28        .map(|(p, &xi)| p.norm() * p.norm() * xi)
29        .sum::<f64>()
30        * dx
31}
32
33pub fn expectation_momentum(psi: &[Complex], dx: f64) -> f64 {
34    use crate::constants::HBAR;
35    let n = psi.len();
36    let mut sum = 0.0;
37    for i in 1..n - 1 {
38        let dpsi = (psi[i + 1] - psi[i - 1]) * Complex::new(1.0 / (2.0 * dx), 0.0);
39        let integrand = psi[i].conj() * dpsi * Complex::new(0.0, -HBAR);
40        sum += integrand.re;
41    }
42    sum * dx
43}
44
45pub fn time_evolve_split_step(
46    psi: &mut [Complex],
47    v: &[f64],
48    dx: f64,
49    dt: f64,
50    mass: f64,
51    steps: usize,
52) {
53    use crate::constants::HBAR;
54    let n = psi.len();
55    let dk = 2.0 * std::f64::consts::PI / (n as f64 * dx);
56
57    for _ in 0..steps {
58        for i in 0..n {
59            let phase = -v[i] * dt / (2.0 * HBAR);
60            psi[i] = psi[i] * Complex::from_polar(1.0, phase);
61        }
62
63        let mut ft = vec![Complex::new(0.0, 0.0); n];
64        for (k, ftk) in ft.iter_mut().enumerate() {
65            for (j, &psij) in psi.iter().enumerate() {
66                let angle = -2.0 * std::f64::consts::PI * (k as f64) * (j as f64) / (n as f64);
67                *ftk = *ftk + psij * Complex::from_polar(1.0, angle);
68            }
69        }
70
71        for (k, ftk) in ft.iter_mut().enumerate() {
72            let kk = if k <= n / 2 {
73                k as f64 * dk
74            } else {
75                (k as f64 - n as f64) * dk
76            };
77            let phase = -HBAR * kk * kk * dt / (2.0 * mass);
78            *ftk = *ftk * Complex::from_polar(1.0, phase);
79        }
80
81        for (j, psij) in psi.iter_mut().enumerate() {
82            *psij = Complex::new(0.0, 0.0);
83            for (k, &ftk) in ft.iter().enumerate() {
84                let angle = 2.0 * std::f64::consts::PI * (k as f64) * (j as f64) / (n as f64);
85                *psij = *psij + ftk * Complex::from_polar(1.0, angle);
86            }
87            *psij = *psij * Complex::new(1.0 / n as f64, 0.0);
88        }
89
90        for i in 0..n {
91            let phase = -v[i] * dt / (2.0 * HBAR);
92            psi[i] = psi[i] * Complex::from_polar(1.0, phase);
93        }
94    }
95}
96
97pub fn inner_product(psi: &[Complex], phi: &[Complex], dx: f64) -> Complex {
98    let mut sum = Complex::new(0.0, 0.0);
99    for i in 0..psi.len() {
100        sum = sum + psi[i].conj() * phi[i];
101    }
102    sum * Complex::new(dx, 0.0)
103}
104
105pub fn transition_probability(psi_initial: &[Complex], psi_final: &[Complex], dx: f64) -> f64 {
106    let ip = inner_product(psi_final, psi_initial, dx);
107    ip.norm() * ip.norm()
108}
sciforge_lib/physics/quantum/wavefunctions.rs

sciforge_lib/physics/quantum/
wavefunctions.rs
