quantrs2_device/photonic/
squeezed_states.rs

1//! Squeezed State Generation and Manipulation
2//!
3//! This module implements advanced squeezed state operations for photonic quantum computing,
4//! including generation, characterization, and application of squeezed light.
5
6use serde::{Deserialize, Serialize};
7use std::collections::HashMap;
8use std::f64::consts::{E, PI};
9use thiserror::Error;
10
11use super::continuous_variable::{CVError, CVResult, Complex, GaussianState};
12use super::{PhotonicMode, PhotonicSystemType};
13use crate::DeviceResult;
14use scirs2_core::random::prelude::*;
15
16/// Squeezed state operation errors
17#[derive(Error, Debug)]
18pub enum SqueezedStateError {
19    #[error("Invalid squeezing parameters: {0}")]
20    InvalidParameters(String),
21    #[error("Squeezing generation failed: {0}")]
22    GenerationFailed(String),
23    #[error("State characterization failed: {0}")]
24    CharacterizationFailed(String),
25    #[error("Squeezing measurement error: {0}")]
26    MeasurementError(String),
27}
28
29/// Types of squeezed states
30#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
31#[allow(non_snake_case)]
32pub enum SqueezedStateType {
33    /// Single-mode squeezed vacuum
34    SingleMode { squeezing_dB: f64, angle: f64 },
35    /// Two-mode squeezed state
36    TwoMode { squeezing_dB: f64, phase: f64 },
37    /// Multimode squeezed state
38    Multimode { squeezing_matrix: Vec<Vec<f64>> },
39    /// Spin squeezed state
40    SpinSqueezed {
41        collective_spin: f64,
42        squeezing_parameter: f64,
43    },
44    /// Amplitude squeezed state
45    AmplitudeSqueezed { alpha: Complex, squeezing_dB: f64 },
46    /// Phase squeezed state
47    PhaseSqueezed { alpha: Complex, squeezing_dB: f64 },
48}
49
50/// Squeezing generation method
51#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
52pub enum SqueezingMethod {
53    /// Parametric down-conversion
54    ParametricDownConversion {
55        pump_power: f64,
56        crystal_length: f64,
57        phase_matching: Phasematching,
58    },
59    /// Four-wave mixing
60    FourWaveMixing {
61        pump_powers: Vec<f64>,
62        fiber_length: f64,
63        nonlinearity: f64,
64    },
65    /// Kerr squeezing
66    KerrSqueezing {
67        kerr_coefficient: f64,
68        interaction_time: f64,
69    },
70    /// Cavity squeezing
71    CavitySqueezing {
72        cavity_finesse: f64,
73        pump_detuning: f64,
74    },
75}
76
77/// Phase matching conditions
78#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
79pub enum Phasematching {
80    /// Type I phase matching
81    TypeI { crystal_angle: f64 },
82    /// Type II phase matching
83    TypeII {
84        crystal_angle: f64,
85        temperature: f64,
86    },
87    /// Quasi-phase matching
88    QuasiPhaseMatching { poling_period: f64 },
89}
90
91/// Squeezed state characterization
92#[derive(Debug, Clone, Serialize, Deserialize)]
93#[allow(non_snake_case)]
94pub struct SqueezingCharacterization {
95    /// Measured squeezing level (dB)
96    pub measured_squeezing_dB: f64,
97    /// Squeezing angle (radians)
98    pub squeezing_angle: f64,
99    /// Anti-squeezing level (dB)
100    pub anti_squeezing_dB: f64,
101    /// Noise figure
102    pub noise_figure: f64,
103    /// Squeezing bandwidth
104    pub bandwidth_hz: f64,
105    /// Purity measure
106    pub purity: f64,
107    /// Loss estimation
108    pub estimated_loss: f64,
109}
110
111/// Squeezed state generator
112pub struct SqueezedStateGenerator {
113    /// Generation method
114    pub method: SqueezingMethod,
115    /// Operating parameters
116    pub parameters: SqueezingParameters,
117    /// Performance metrics
118    pub performance: GeneratorPerformance,
119    /// Calibration data
120    pub calibration: GeneratorCalibration,
121}
122
123/// Squeezing generation parameters
124#[derive(Debug, Clone, Serialize, Deserialize)]
125#[allow(non_snake_case)]
126pub struct SqueezingParameters {
127    /// Target squeezing level (dB)
128    pub target_squeezing_dB: f64,
129    /// Operating wavelength (nm)
130    pub wavelength_nm: f64,
131    /// Pump power (mW)
132    pub pump_power_mw: f64,
133    /// Temperature (K)
134    pub temperature_k: f64,
135    /// Phase stabilization enabled
136    pub phase_stabilization: bool,
137    /// Feedback control enabled
138    pub feedback_control: bool,
139}
140
141impl Default for SqueezingParameters {
142    fn default() -> Self {
143        Self {
144            target_squeezing_dB: 10.0,
145            wavelength_nm: 1550.0,
146            pump_power_mw: 100.0,
147            temperature_k: 293.0,
148            phase_stabilization: true,
149            feedback_control: true,
150        }
151    }
152}
153
154/// Generator performance metrics
155#[derive(Debug, Clone, Serialize, Deserialize)]
156#[allow(non_snake_case)]
157pub struct GeneratorPerformance {
158    /// Achieved squeezing level (dB)
159    pub achieved_squeezing_dB: f64,
160    /// Squeezing stability (%)
161    pub stability_percent: f64,
162    /// Generation efficiency (%)
163    pub efficiency_percent: f64,
164    /// Output power (μW)
165    pub output_power_uw: f64,
166    /// Noise floor (dB below shot noise)
167    pub noise_floor_dB: f64,
168}
169
170impl Default for GeneratorPerformance {
171    fn default() -> Self {
172        Self {
173            achieved_squeezing_dB: 8.5,
174            stability_percent: 95.0,
175            efficiency_percent: 85.0,
176            output_power_uw: 50.0,
177            noise_floor_dB: -15.0,
178        }
179    }
180}
181
182/// Generator calibration data
183#[derive(Debug, Clone, Serialize, Deserialize)]
184pub struct GeneratorCalibration {
185    /// Power calibration curve
186    pub power_curve: Vec<(f64, f64)>, // (input_power, output_squeezing)
187    /// Temperature dependence
188    pub temperature_curve: Vec<(f64, f64)>, // (temperature, squeezing)
189    /// Phase calibration
190    pub phase_calibration: Vec<(f64, f64)>, // (phase_setting, actual_phase)
191    /// Last calibration timestamp (Unix timestamp)
192    #[serde(with = "instant_serde")]
193    pub last_calibration: std::time::Instant,
194}
195
196mod instant_serde {
197    use serde::{Deserialize, Deserializer, Serialize, Serializer};
198    use std::time::{Duration, Instant, SystemTime, UNIX_EPOCH};
199
200    pub fn serialize<S>(instant: &Instant, serializer: S) -> Result<S::Ok, S::Error>
201    where
202        S: Serializer,
203    {
204        let duration_since_epoch = SystemTime::now()
205            .duration_since(UNIX_EPOCH)
206            .expect("System time should be after UNIX epoch");
207        duration_since_epoch.as_secs().serialize(serializer)
208    }
209
210    pub fn deserialize<'de, D>(deserializer: D) -> Result<Instant, D::Error>
211    where
212        D: Deserializer<'de>,
213    {
214        let secs = u64::deserialize(deserializer)?;
215        Ok(Instant::now()) // Simplified - in practice would need proper epoch handling
216    }
217}
218
219impl Default for GeneratorCalibration {
220    fn default() -> Self {
221        Self {
222            power_curve: vec![(0.0, 0.0), (50.0, 5.0), (100.0, 10.0), (200.0, 12.0)],
223            temperature_curve: vec![(273.0, 8.0), (293.0, 10.0), (313.0, 9.0)],
224            phase_calibration: vec![(0.0, 0.0), (PI / 2.0, PI / 2.0), (PI, PI)],
225            last_calibration: std::time::Instant::now(),
226        }
227    }
228}
229
230#[allow(non_snake_case)]
231impl SqueezedStateGenerator {
232    pub fn new(method: SqueezingMethod) -> Self {
233        Self {
234            method,
235            parameters: SqueezingParameters::default(),
236            performance: GeneratorPerformance::default(),
237            calibration: GeneratorCalibration::default(),
238        }
239    }
240
241    /// Generate a squeezed state
242    pub fn generate_squeezed_state(
243        &mut self,
244        state_type: SqueezedStateType,
245        num_modes: usize,
246    ) -> Result<GaussianState, SqueezedStateError> {
247        match state_type {
248            SqueezedStateType::SingleMode {
249                squeezing_dB,
250                angle,
251            } => self.generate_single_mode_squeezed(squeezing_dB, angle, num_modes),
252            SqueezedStateType::TwoMode {
253                squeezing_dB,
254                phase,
255            } => self.generate_two_mode_squeezed(squeezing_dB, phase, num_modes),
256            SqueezedStateType::AmplitudeSqueezed {
257                alpha,
258                squeezing_dB,
259            } => self.generate_amplitude_squeezed(alpha, squeezing_dB, num_modes),
260            SqueezedStateType::PhaseSqueezed {
261                alpha,
262                squeezing_dB,
263            } => self.generate_phase_squeezed(alpha, squeezing_dB, num_modes),
264            _ => Err(SqueezedStateError::GenerationFailed(format!(
265                "State type {state_type:?} not yet implemented"
266            ))),
267        }
268    }
269
270    /// Generate single-mode squeezed vacuum
271    fn generate_single_mode_squeezed(
272        &mut self,
273        squeezing_dB: f64,
274        angle: f64,
275        num_modes: usize,
276    ) -> Result<GaussianState, SqueezedStateError> {
277        if squeezing_dB < 0.0 {
278            return Err(SqueezedStateError::InvalidParameters(
279                "Squeezing must be non-negative".to_string(),
280            ));
281        }
282
283        // Convert dB to linear squeezing parameter
284        let squeezing_r = squeezing_dB * (10.0_f64.ln() / 20.0);
285
286        // Account for realistic limitations
287        let realistic_r = self.apply_realistic_limitations(squeezing_r)?;
288
289        // Create squeezed vacuum state
290        match GaussianState::squeezed_vacuum(realistic_r, angle, 0, num_modes) {
291            Ok(state) => {
292                self.update_performance_metrics(squeezing_dB, realistic_r);
293                Ok(state)
294            }
295            Err(e) => Err(SqueezedStateError::GenerationFailed(format!(
296                "Failed to generate squeezed state: {e:?}"
297            ))),
298        }
299    }
300
301    /// Generate two-mode squeezed state
302    fn generate_two_mode_squeezed(
303        &mut self,
304        squeezing_dB: f64,
305        phase: f64,
306        num_modes: usize,
307    ) -> Result<GaussianState, SqueezedStateError> {
308        if num_modes < 2 {
309            return Err(SqueezedStateError::InvalidParameters(
310                "Two-mode squeezing requires at least 2 modes".to_string(),
311            ));
312        }
313
314        let squeezing_r = squeezing_dB * (10.0_f64.ln() / 20.0);
315        let realistic_r = self.apply_realistic_limitations(squeezing_r)?;
316
317        let mut state = GaussianState::vacuum(num_modes);
318
319        match state.two_mode_squeeze(realistic_r, phase, 0, 1) {
320            Ok(()) => {
321                self.update_performance_metrics(squeezing_dB, realistic_r);
322                Ok(state)
323            }
324            Err(e) => Err(SqueezedStateError::GenerationFailed(format!(
325                "Failed to generate two-mode squeezed state: {e:?}"
326            ))),
327        }
328    }
329
330    /// Generate amplitude squeezed state
331    fn generate_amplitude_squeezed(
332        &mut self,
333        alpha: Complex,
334        squeezing_dB: f64,
335        num_modes: usize,
336    ) -> Result<GaussianState, SqueezedStateError> {
337        let squeezing_r = squeezing_dB * (10.0_f64.ln() / 20.0);
338        let realistic_r = self.apply_realistic_limitations(squeezing_r)?;
339
340        // Create coherent state
341        let mut state = match GaussianState::coherent(alpha, 0, num_modes) {
342            Ok(s) => s,
343            Err(e) => {
344                return Err(SqueezedStateError::GenerationFailed(format!(
345                    "Failed to create coherent state: {e:?}"
346                )))
347            }
348        };
349
350        // Apply amplitude squeezing (squeezing in phase with displacement)
351        let squeezing_angle = alpha.phase();
352        match state.squeeze(realistic_r, squeezing_angle, 0) {
353            Ok(()) => {
354                self.update_performance_metrics(squeezing_dB, realistic_r);
355                Ok(state)
356            }
357            Err(e) => Err(SqueezedStateError::GenerationFailed(format!(
358                "Failed to apply amplitude squeezing: {e:?}"
359            ))),
360        }
361    }
362
363    /// Generate phase squeezed state
364    fn generate_phase_squeezed(
365        &mut self,
366        alpha: Complex,
367        squeezing_dB: f64,
368        num_modes: usize,
369    ) -> Result<GaussianState, SqueezedStateError> {
370        let squeezing_r = squeezing_dB * (10.0_f64.ln() / 20.0);
371        let realistic_r = self.apply_realistic_limitations(squeezing_r)?;
372
373        // Create coherent state
374        let mut state = match GaussianState::coherent(alpha, 0, num_modes) {
375            Ok(s) => s,
376            Err(e) => {
377                return Err(SqueezedStateError::GenerationFailed(format!(
378                    "Failed to create coherent state: {e:?}"
379                )))
380            }
381        };
382
383        // Apply phase squeezing (squeezing perpendicular to displacement)
384        let squeezing_angle = alpha.phase() + PI / 2.0;
385        match state.squeeze(realistic_r, squeezing_angle, 0) {
386            Ok(()) => {
387                self.update_performance_metrics(squeezing_dB, realistic_r);
388                Ok(state)
389            }
390            Err(e) => Err(SqueezedStateError::GenerationFailed(format!(
391                "Failed to apply phase squeezing: {e:?}"
392            ))),
393        }
394    }
395
396    /// Apply realistic limitations to squeezing
397    fn apply_realistic_limitations(&self, ideal_r: f64) -> Result<f64, SqueezedStateError> {
398        // Limit maximum squeezing based on method
399        let max_r = match &self.method {
400            SqueezingMethod::ParametricDownConversion { .. } => 2.3, // ~20 dB typical maximum
401            SqueezingMethod::FourWaveMixing { .. } => 1.15,          // ~10 dB typical maximum
402            SqueezingMethod::KerrSqueezing { .. } => 0.46,           // ~4 dB typical maximum
403            SqueezingMethod::CavitySqueezing { .. } => 3.45,         // ~30 dB possible in cavities
404        };
405
406        if ideal_r > max_r {
407            return Err(SqueezedStateError::InvalidParameters(format!(
408                "Requested squeezing {} exceeds maximum {} for method {:?}",
409                ideal_r, max_r, self.method
410            )));
411        }
412
413        // Apply efficiency and loss factors
414        let efficiency = self.performance.efficiency_percent / 100.0;
415        let realistic_r = ideal_r * efficiency;
416
417        // Add noise and imperfections
418        let noise_factor = 1.0 - (self.performance.noise_floor_dB / 20.0);
419        let final_r = realistic_r * noise_factor;
420
421        Ok(final_r.max(0.0))
422    }
423
424    /// Update performance metrics based on generation
425    fn update_performance_metrics(&mut self, target_dB: f64, achieved_r: f64) {
426        let achieved_dB = achieved_r * (20.0 / 10.0_f64.ln());
427        self.performance.achieved_squeezing_dB = achieved_dB;
428
429        // Update efficiency estimate
430        let efficiency = achieved_dB / target_dB;
431        self.performance.efficiency_percent = (efficiency * 100.0).min(100.0);
432
433        // Simulate stability fluctuations
434        let stability_noise = (thread_rng().gen::<f64>() - 0.5) * 0.1;
435        self.performance.stability_percent = (95.0 + stability_noise).clamp(90.0, 99.0);
436    }
437
438    /// Characterize generated squeezed state
439    pub fn characterize_state(
440        &self,
441        state: &GaussianState,
442        mode: usize,
443    ) -> Result<SqueezingCharacterization, SqueezedStateError> {
444        if mode >= state.num_modes {
445            return Err(SqueezedStateError::CharacterizationFailed(format!(
446                "Mode {mode} does not exist"
447            )));
448        }
449
450        // Extract squeezing information from covariance matrix
451        let x_idx = 2 * mode;
452        let p_idx = 2 * mode + 1;
453
454        let var_x = state.covariance[x_idx][x_idx];
455        let var_p = state.covariance[p_idx][p_idx];
456        let cov_xp = state.covariance[x_idx][p_idx];
457
458        // Calculate eigenvalues of covariance matrix
459        let trace = var_x + var_p;
460        let det = var_x.mul_add(var_p, -(cov_xp * cov_xp));
461        let discriminant = trace.mul_add(trace, -(4.0 * det));
462
463        if discriminant < 0.0 {
464            return Err(SqueezedStateError::CharacterizationFailed(
465                "Invalid covariance matrix".to_string(),
466            ));
467        }
468
469        let sqrt_discriminant = discriminant.sqrt();
470        let eigenval1 = f64::midpoint(trace, sqrt_discriminant);
471        let eigenval2 = (trace - sqrt_discriminant) / 2.0;
472
473        let min_eigenval = eigenval1.min(eigenval2);
474        let max_eigenval = eigenval1.max(eigenval2);
475
476        // Calculate squeezing and anti-squeezing in dB
477        let squeezing_dB = -10.0 * (min_eigenval / 0.5).log10();
478        let anti_squeezing_dB = 10.0 * (max_eigenval / 0.5).log10();
479
480        // Calculate squeezing angle
481        let squeezing_angle = if cov_xp.abs() > 1e-10 {
482            0.5 * (2.0 * cov_xp / (var_x - var_p)).atan()
483        } else if var_x < var_p {
484            0.0
485        } else {
486            PI / 2.0
487        };
488
489        // Estimate purity
490        let purity = 1.0 / (4.0 * det).sqrt();
491
492        // Estimate loss (simplified)
493        let theoretical_min =
494            0.5 * (-self.performance.achieved_squeezing_dB * 10.0_f64.ln() / 10.0).exp();
495        let estimated_loss = (min_eigenval - theoretical_min) / (0.5 - theoretical_min);
496
497        Ok(SqueezingCharacterization {
498            measured_squeezing_dB: squeezing_dB,
499            squeezing_angle,
500            anti_squeezing_dB,
501            noise_figure: anti_squeezing_dB - squeezing_dB,
502            bandwidth_hz: 1e6, // Placeholder - would depend on specific measurement
503            purity: purity.min(1.0),
504            estimated_loss: estimated_loss.clamp(0.0, 1.0),
505        })
506    }
507
508    /// Optimize generator parameters
509    pub fn optimize_parameters(
510        &mut self,
511        target_squeezing_dB: f64,
512    ) -> Result<SqueezingParameters, SqueezedStateError> {
513        // Simple optimization based on calibration curves
514        let mut optimized = self.parameters.clone();
515        optimized.target_squeezing_dB = target_squeezing_dB;
516
517        // Find optimal power from calibration curve
518        let optimal_power = self.find_optimal_power(target_squeezing_dB);
519        optimized.pump_power_mw = optimal_power;
520
521        // Adjust temperature for stability
522        let optimal_temp = self.find_optimal_temperature(target_squeezing_dB);
523        optimized.temperature_k = optimal_temp;
524
525        // Enable feedback for high squeezing levels
526        optimized.feedback_control = target_squeezing_dB >= 15.0;
527        optimized.phase_stabilization = target_squeezing_dB >= 10.0;
528
529        self.parameters = optimized.clone();
530        Ok(optimized)
531    }
532
533    /// Find optimal pump power from calibration
534    fn find_optimal_power(&self, target_squeezing_dB: f64) -> f64 {
535        for (power, squeezing) in &self.calibration.power_curve {
536            if *squeezing >= target_squeezing_dB {
537                return *power;
538            }
539        }
540        // Extrapolate if beyond calibration range
541        self.calibration
542            .power_curve
543            .last()
544            .map_or(100.0, |(p, _)| *p * 1.5)
545    }
546
547    /// Find optimal temperature from calibration
548    fn find_optimal_temperature(&self, _target_squeezing_dB: f64) -> f64 {
549        // Find temperature with maximum squeezing
550        self.calibration
551            .temperature_curve
552            .iter()
553            .max_by(|(_, s1), (_, s2)| s1.partial_cmp(s2).unwrap_or(std::cmp::Ordering::Equal))
554            .map_or(293.0, |(t, _)| *t)
555    }
556
557    /// Create multimode squeezed state
558    pub fn create_multimode_squeezed(
559        &mut self,
560        num_modes: usize,
561        squeezing_matrix: &[Vec<f64>],
562    ) -> Result<GaussianState, SqueezedStateError> {
563        if squeezing_matrix.len() != num_modes
564            || squeezing_matrix.iter().any(|row| row.len() != num_modes)
565        {
566            return Err(SqueezedStateError::InvalidParameters(
567                "Squeezing matrix dimensions don't match number of modes".to_string(),
568            ));
569        }
570
571        let mut state = GaussianState::vacuum(num_modes);
572
573        // Apply squeezing transformations based on matrix
574        for i in 0..num_modes {
575            for j in i + 1..num_modes {
576                let squeezing_strength = squeezing_matrix[i][j];
577                if squeezing_strength.abs() > 1e-10 {
578                    let realistic_r = self.apply_realistic_limitations(squeezing_strength)?;
579                    if let Err(e) = state.two_mode_squeeze(realistic_r, 0.0, i, j) {
580                        return Err(SqueezedStateError::GenerationFailed(format!(
581                            "Multimode squeezing failed: {e:?}"
582                        )));
583                    }
584                }
585            }
586        }
587
588        Ok(state)
589    }
590}
591
592/// Squeezed state measurement apparatus
593pub struct SqueezedStateMeasurement {
594    /// Homodyne detection efficiency
595    pub homodyne_efficiency: f64,
596    /// Electronic noise level
597    pub electronic_noise: f64,
598    /// Dark count rate
599    pub dark_count_rate: f64,
600    /// Measurement bandwidth
601    pub bandwidth_hz: f64,
602}
603
604#[allow(non_snake_case)]
605impl SqueezedStateMeasurement {
606    pub const fn new() -> Self {
607        Self {
608            homodyne_efficiency: 0.95,
609            electronic_noise: 0.01,
610            dark_count_rate: 1000.0, // Hz
611            bandwidth_hz: 1e6,
612        }
613    }
614
615    /// Measure squeezing via homodyne detection
616    pub fn measure_squeezing_homodyne(
617        &self,
618        state: &GaussianState,
619        mode: usize,
620        local_oscillator_phase: f64,
621        measurement_time: f64,
622    ) -> Result<Vec<f64>, SqueezedStateError> {
623        if mode >= state.num_modes {
624            return Err(SqueezedStateError::MeasurementError(format!(
625                "Mode {mode} does not exist"
626            )));
627        }
628
629        let x_idx = 2 * mode;
630        let p_idx = 2 * mode + 1;
631
632        // Calculate quadrature being measured
633        let cos_phi = local_oscillator_phase.cos();
634        let sin_phi = local_oscillator_phase.sin();
635
636        let mean_quad = cos_phi.mul_add(state.mean[x_idx], sin_phi * state.mean[p_idx]);
637        let var_quad = (2.0 * cos_phi * sin_phi).mul_add(
638            state.covariance[x_idx][p_idx],
639            (cos_phi * cos_phi).mul_add(
640                state.covariance[x_idx][x_idx],
641                sin_phi * sin_phi * state.covariance[p_idx][p_idx],
642            ),
643        );
644
645        // Account for detection efficiency
646        let effective_variance = self
647            .homodyne_efficiency
648            .mul_add(var_quad, (1.0 - self.homodyne_efficiency) * 0.5)
649            + self.electronic_noise;
650
651        // Generate measurement samples
652        let num_samples = (measurement_time * self.bandwidth_hz) as usize;
653        let mut measurements = Vec::with_capacity(num_samples);
654
655        for _ in 0..num_samples {
656            // Add dark counts
657            let dark_noise = if thread_rng().gen::<f64>() < self.dark_count_rate / self.bandwidth_hz
658            {
659                (thread_rng().gen::<f64>() - 0.5) * 0.1
660            } else {
661                0.0
662            };
663
664            // Generate Gaussian measurement result
665            let measurement =
666                self.generate_gaussian_sample(mean_quad, effective_variance) + dark_noise;
667            measurements.push(measurement);
668        }
669
670        Ok(measurements)
671    }
672
673    /// Generate Gaussian random sample
674    fn generate_gaussian_sample(&self, mean: f64, variance: f64) -> f64 {
675        // Box-Muller transform
676        let u1 = thread_rng().gen::<f64>();
677        let u2 = thread_rng().gen::<f64>();
678        let z = (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos();
679        variance.sqrt().mul_add(z, mean)
680    }
681
682    /// Analyze measurement data for squeezing
683    pub fn analyze_squeezing(
684        &self,
685        measurements: &[f64],
686        reference_measurements: &[f64],
687    ) -> Result<f64, SqueezedStateError> {
688        if measurements.is_empty() || reference_measurements.is_empty() {
689            return Err(SqueezedStateError::MeasurementError(
690                "Insufficient measurement data".to_string(),
691            ));
692        }
693
694        // Calculate variances
695        let squeezed_variance = self.calculate_variance(measurements);
696        let reference_variance = self.calculate_variance(reference_measurements);
697
698        // Squeezing in dB relative to reference (shot noise)
699        let squeezing_dB = 10.0 * (squeezed_variance / reference_variance).log10();
700
701        Ok(-squeezing_dB) // Negative because squeezing reduces noise
702    }
703
704    /// Calculate sample variance
705    fn calculate_variance(&self, data: &[f64]) -> f64 {
706        let mean = data.iter().sum::<f64>() / data.len() as f64;
707        let variance =
708            data.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / (data.len() - 1) as f64;
709        variance
710    }
711}
712
713impl Default for SqueezedStateMeasurement {
714    fn default() -> Self {
715        Self::new()
716    }
717}
718
719#[cfg(test)]
720mod tests {
721    use super::*;
722
723    #[test]
724    fn test_squeezed_state_generator_creation() {
725        let method = SqueezingMethod::ParametricDownConversion {
726            pump_power: 100.0,
727            crystal_length: 10.0,
728            phase_matching: Phasematching::TypeI {
729                crystal_angle: 45.0,
730            },
731        };
732
733        let generator = SqueezedStateGenerator::new(method);
734        assert_eq!(generator.parameters.target_squeezing_dB, 10.0);
735    }
736
737    #[test]
738    fn test_single_mode_squeezed_generation() {
739        let method = SqueezingMethod::ParametricDownConversion {
740            pump_power: 100.0,
741            crystal_length: 10.0,
742            phase_matching: Phasematching::TypeI {
743                crystal_angle: 45.0,
744            },
745        };
746
747        let mut generator = SqueezedStateGenerator::new(method);
748
749        let state_type = SqueezedStateType::SingleMode {
750            squeezing_dB: 10.0,
751            angle: 0.0,
752        };
753
754        let result = generator.generate_squeezed_state(state_type, 1);
755        assert!(result.is_ok());
756
757        let state = result.expect("Single-mode squeezed state generation should succeed");
758        assert_eq!(state.num_modes, 1);
759
760        // Check that squeezing was applied (X quadrature should be squeezed)
761        assert!(state.covariance[0][0] < 0.5); // Squeezed quadrature
762        assert!(state.covariance[1][1] > 0.5); // Anti-squeezed quadrature
763    }
764
765    #[test]
766    fn test_two_mode_squeezed_generation() {
767        let method = SqueezingMethod::FourWaveMixing {
768            pump_powers: vec![50.0, 50.0],
769            fiber_length: 100.0,
770            nonlinearity: 0.001,
771        };
772
773        let mut generator = SqueezedStateGenerator::new(method);
774
775        let state_type = SqueezedStateType::TwoMode {
776            squeezing_dB: 8.0,
777            phase: 0.0,
778        };
779
780        let result = generator.generate_squeezed_state(state_type, 2);
781        assert!(result.is_ok());
782
783        let state = result.expect("Two-mode squeezed state generation should succeed");
784        assert_eq!(state.num_modes, 2);
785    }
786
787    #[test]
788    fn test_squeezing_characterization() {
789        let method = SqueezingMethod::ParametricDownConversion {
790            pump_power: 100.0,
791            crystal_length: 10.0,
792            phase_matching: Phasematching::TypeI {
793                crystal_angle: 45.0,
794            },
795        };
796
797        let mut generator = SqueezedStateGenerator::new(method);
798
799        let state_type = SqueezedStateType::SingleMode {
800            squeezing_dB: 10.0,
801            angle: 0.0,
802        };
803
804        let state = generator
805            .generate_squeezed_state(state_type, 1)
806            .expect("Squeezed state generation should succeed");
807        let characterization = generator
808            .characterize_state(&state, 0)
809            .expect("State characterization should succeed");
810
811        assert!(characterization.measured_squeezing_dB > 0.0);
812        assert!(characterization.anti_squeezing_dB > 0.0);
813        assert!(characterization.purity > 0.0 && characterization.purity <= 1.0);
814    }
815
816    #[test]
817    fn test_measurement_apparatus() {
818        let measurement = SqueezedStateMeasurement::new();
819
820        // Create a simple squeezed state for testing
821        let state = GaussianState::squeezed_vacuum(1.0, 0.0, 0, 1)
822            .expect("Squeezed vacuum state creation should succeed");
823
824        let measurements = measurement
825            .measure_squeezing_homodyne(&state, 0, 0.0, 0.1)
826            .expect("Homodyne measurement should succeed");
827
828        assert!(!measurements.is_empty());
829        assert!(measurements.len() > 10000); // Should have many samples
830    }
831
832    #[test]
833    fn test_parameter_optimization() {
834        let method = SqueezingMethod::ParametricDownConversion {
835            pump_power: 100.0,
836            crystal_length: 10.0,
837            phase_matching: Phasematching::TypeI {
838                crystal_angle: 45.0,
839            },
840        };
841
842        let mut generator = SqueezedStateGenerator::new(method);
843
844        let optimized = generator
845            .optimize_parameters(15.0)
846            .expect("Parameter optimization should succeed");
847        assert_eq!(optimized.target_squeezing_dB, 15.0);
848        assert!(optimized.feedback_control); // Should be enabled for high squeezing
849        assert!(optimized.phase_stabilization);
850    }
851}