crseo 2.5.3

Cuda Engined Optics Rust Interface
Documentation
use crseo::analytic::*;
use skyangle::SkyAngle::*;

fn main() {
    let m1 = Conic::gmt_m1();
    let m2 = Conic::gmt_m2();

    println!("CHIEF RAY:");
    let mut ray = new_ray()
        .polar_direction_vector(Arcminute(10f64).to_radians(), 0f64)
        .build();
    println!("Init   : {}", ray);
    m1.reflect(&mut ray);
    println!("Reflect: {}", ray);
    ray.trace_to(&m2);
    println!("Trace  : {}", ray);
    m2.reflect(&mut ray);
    println!("Reflect: {}", ray);
    let z_xpupil = ray.solve_for_z(0f64, 0f64);
    println!("Exit pupil: {:.9}m", z_xpupil);

    println!("MARGINAL RAY:");
    let mut ray = new_ray()
        .point_of_origin(m1.height_at([10f64, 0f64, 0f64]))
        //        .polar_direction_vector(Arcminute(10f64).to_radians(), 0f64)
        .build();
    println!("Init   : {}", ray);
    m1.reflect(&mut ray);
    println!("Reflect: {}", ray);
    let z_gregorian = ray.solve_for_z(0f64, 0f64);
    println!("Gregorian focus: {:.9}m", z_gregorian);
    ray.trace_to(&m2);
    println!("Trace  : {}", ray);
    m2.reflect(&mut ray);
    println!("Reflect: {}", ray);
    let z_focal_plane = ray.solve_for_z(0f64, 0f64);
    println!("Focal plane: {:.9}m", z_focal_plane);

    println!("M1 walk");
    let gmt = Gmt::new();
    let m: Vec<Vector> = (1..11).map(|x| [x as f64, 0., 0.]).collect();
    let rays = gmt.focal_point(m.clone(), Arcminute(10f64).to_radians(), 0f64);
    rays.iter()
        .enumerate()
        .for_each(|x| println!("#{:3}: {}", x.0, x.1));

    println!("Field walk");
    for z in 1..11 {
        let rays = gmt.focal_point(vec![[1f64;3]], Arcsecond(z as f64).to_radians(), 0f64);
        println!("#{:3}\": {}", z, rays[0]);
    }

    println!("Field curvature");
    for z in 1..11 {
        let rays = gmt.focal_point(m.clone(), Arcminute(z as f64).to_radians(), 0f64);
        let z_mean = rays.iter().map(|r| r.p[2]).sum::<f64>() / rays.len() as f64;
        let z_rms = (rays.iter().map(|r| (r.p[2] - z_mean).powf(2.)).sum::<f64>()
            / rays.len() as f64)
            .sqrt();
        let mut p = rays[0].p;
        p[2] = z_mean - z_focal_plane;
        let r = -0.5 * p.norm_square() / p[2];
        println!(
            "z: {:2}': {:.9}m +/- {:3.0}nm ; curvature: {:.9}m",
            z,
            z_mean,
            z_rms * 1e9,
            r
        );
    }
}