frinZ 4.6.0

fringe search for Yamaguchi Interferometer and Japanese VLBI Network
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
89
90
91
92

use std::error::Error;
use std::fs;
use std::io::Cursor;
use std::path::Path;

use crate::args::Args;
use crate::header::parse_header;
use crate::input_support::read_input_bytes;
use crate::npy_output::{npz_sidecar_path, write_complex_2d, NpyMeta};
use crate::plot;
use crate::read::read_visibility_data;
use num_complex::Complex;
type C32 = Complex<f32>;

/// Executes the raw visibility plotting.
pub fn run_raw_visibility_plot(args: &Args) -> Result<(), Box<dyn Error>> {
    //println!("# Starting raw visibility plotting...");

    let input_path = args.input.as_ref().unwrap();

    // --- Create Output Directory ---
    let parent_dir = input_path.parent().unwrap_or_else(|| Path::new(""));
    let output_dir = parent_dir.join("frinZ").join("rawvis");
    fs::create_dir_all(&output_dir)?;
    let base_filename = input_path.file_stem().unwrap().to_str().unwrap();

    let buffer = read_input_bytes(input_path)?;
    let mut cursor = Cursor::new(buffer.as_slice());

    let header = parse_header(&mut cursor)?;

    let mut all_spectra: Vec<Vec<C32>> = Vec::new();
    for l1 in 0..header.number_of_sector {
        let (complex_vec, _, _) = match read_visibility_data(
            &mut cursor,
            &header,
            1,  // length in sectors
            0,  // skip in sectors
            l1, // loop_idx, which acts as sector index here
            false,
            &[], // Add empty pp_flag_ranges
        ) {
            Ok(data) => data,
            Err(_) => break, // Stop if we can't read more data
        };
        if complex_vec.is_empty() {
            eprintln!("Warning: Empty sector {} found, stopping read.", l1);
            break;
        }
        all_spectra.push(complex_vec);
    }

    if all_spectra.is_empty() {
        eprintln!("No visibility data found in the file.");
        return Ok(());
    }

    let amp_heatmap_filepath = output_dir.join(format!("{}_heatmap_amp.png", base_filename));
    let phase_heatmap_filepath = output_dir.join(format!("{}_heatmap_phase.png", base_filename));
    for legacy_filename in [
        format!("{}_heatmap_amp_phase.png", base_filename),
        format!("{}_scatter_real_imag.png", base_filename),
        format!("{}_scatter_amp_phase.png", base_filename),
    ] {
        let _ = fs::remove_file(output_dir.join(legacy_filename));
    }

    // Use a default sigma of 0.0 for blurring, as in the original frinZrawvis.
    plot::plot_spectrum_amplitude_heatmap(&amp_heatmap_filepath, &all_spectra, 0.0)?;
    plot::plot_spectrum_phase_heatmap(&phase_heatmap_filepath, &all_spectra, 0.0)?;
    let rows = all_spectra.len();
    let cols = all_spectra.first().map_or(0, Vec::len);
    if args.npz && rows > 0 && cols > 0 && all_spectra.iter().all(|row| row.len() == cols) {
        let time_axis: Vec<f64> = (0..rows).map(|index| index as f64).collect();
        let channel_axis: Vec<f64> = (0..cols).map(|index| index as f64).collect();
        write_complex_2d(
            &npz_sidecar_path(&amp_heatmap_filepath, "raw_visibility"),
            NpyMeta::new(
                "raw_visibility",
                header.fft_point as u32,
                header.number_of_sector as u32,
            )
            .axes("sector", "", "channel", ""),
            (rows, cols),
            all_spectra.iter().flatten().copied(),
            &time_axis,
            &channel_axis,
        )?;
    }

    Ok(())
}