rsspice 0.1.0

//
// GENERATED FILE
//

use super::*;
use crate::SpiceContext;
use f2rust_std::*;

struct SaveVars {
    TL: f64,
    G: StackArray<f64, 15>,
    REFPOS: StackArray<f64, 3>,
    REFVEL: StackArray<f64, 3>,
    DT: StackArray2D<f64, 45>,
    KQMAX1: i32,
    KQ: StackArray<i32, 3>,
    FC: StackArray<f64, 14>,
    SUM: f64,
    DELTA: f64,
    TP: f64,
    WC: StackArray<f64, 13>,
    W: StackArray<f64, 17>,
    MQ2: i32,
    KS1: i32,
    KS: i32,
    KQQ: i32,
    JX: i32,
}

impl SaveInit for SaveVars {
    fn new() -> Self {
        let mut TL: f64 = 0.0;
        let mut G = StackArray::<f64, 15>::new(1..=15);
        let mut REFPOS = StackArray::<f64, 3>::new(1..=3);
        let mut REFVEL = StackArray::<f64, 3>::new(1..=3);
        let mut DT = StackArray2D::<f64, 45>::new(1..=15, 1..=3);
        let mut KQMAX1: i32 = 0;
        let mut KQ = StackArray::<i32, 3>::new(1..=3);
        let mut FC = StackArray::<f64, 14>::new(1..=14);
        let mut SUM: f64 = 0.0;
        let mut DELTA: f64 = 0.0;
        let mut TP: f64 = 0.0;
        let mut WC = StackArray::<f64, 13>::new(1..=13);
        let mut W = StackArray::<f64, 17>::new(1..=17);
        let mut MQ2: i32 = 0;
        let mut KS1: i32 = 0;
        let mut KS: i32 = 0;
        let mut KQQ: i32 = 0;
        let mut JX: i32 = 0;

        {
            use f2rust_std::data::Val;

            let mut clist = [Val::D(1.0)].into_iter();
            FC[1] = clist.next().unwrap().into_f64();

            debug_assert!(clist.next().is_none(), "DATA not fully initialised");
        }

        Self {
            TL,
            G,
            REFPOS,
            REFVEL,
            DT,
            KQMAX1,
            KQ,
            FC,
            SUM,
            DELTA,
            TP,
            WC,
            W,
            MQ2,
            KS1,
            KS,
            KQQ,
            JX,
        }
    }
}

/// S/P Kernel, evaluate, type 1
///
/// Evaluate a single SPK data record from a segment of type 1
/// (Difference Lines).
///
/// # Required Reading
///
/// * [SPK](crate::required_reading::spk)
///
/// # Brief I/O
///
/// ```text
///  VARIABLE  I/O  DESCRIPTION
///  --------  ---  --------------------------------------------------
///  ET         I   Target epoch.
///  RECORD     I   Data record.
///  STATE      O   State (position and velocity).
/// ```
///
/// # Detailed Input
///
/// ```text
///  ET       is a target epoch, at which a state vector is to
///           be computed.
///
///  RECORD   is a data record which, when evaluated at epoch ET,
///           will give the state (position and velocity) of some
///           body, relative to some center, in some inertial
///           reference frame.
/// ```
///
/// # Detailed Output
///
/// ```text
///  STATE    is the state. Units are km and km/sec.
/// ```
///
/// # Particulars
///
/// ```text
///  The exact format and structure of type 1 (difference lines)
///  segments are described in the SPK Required Reading file.
///
///  Difference lines (DL's) are generated by JPL navigation
///  system programs P and PV. Each data record is equivalent
///  to the (slightly rearranged) 'P' portion of a NAVIO PV file
///  data record.
///
///  SPKE01 is a specialized version of Fred Krogh's subroutine DAINT.
///  Only the calling sequence has been changed.
///
///  Because the original version was undocumented, only Fred
///  knows how this really works.
/// ```
///
/// # Author and Institution
///
/// ```text
///  J. Diaz del Rio    (ODC Space)
///  F.T. Krogh         (JPL)
///  H.A. Neilan        (JPL)
///  W.L. Taber         (JPL)
///  I.M. Underwood     (JPL)
/// ```
///
/// # Version
///
/// ```text
/// -    SPICELIB Version 1.2.0, 14-APR-2021 (JDR)
///
///         Added IMPLICIT NONE statement.
///
///         Edited the header to comply with NAIF standard. Moved SPK
///         required reading from $Literature_References to
///         $Required_Reading section.
///
/// -    SPICELIB Version 1.1.0, 14-FEB-1997 (WLT)
///
///         The goto's were removed and loop and if structures
///         revealed. We still don't know exactly what's going
///         on, but at least the bones of this routine have been
///         cleaned off and are ready for assembly. (WLT)
///
/// -    SPICELIB Version 1.0.4, 30-OCT-1996 (WLT)
///
///         Removed redundant SAVE statements from the declaration
///         section. Thanks to Steve Schlaifer for finding this
///         error.
///
/// -    SPICELIB Version 1.0.3, 10-MAR-1992 (WLT)
///
///         Comment section for permuted index source lines was added
///         following the header.
///
/// -    SPICELIB Version 1.0.2, 23-AUG-1991 (HAN)
///
///         SPK01 was removed from the $Required_Reading section of the
///         header. The information in the SPK01 Required Reading file
///         is now part of the SPK Required Reading file.
///
/// -    SPICELIB Version 1.0.1, 22-MAR-1990 (HAN)
///
///         Literature references added to the header.
///
/// -    SPICELIB Version 1.0.0, 31-JAN-1990 (IMU) (FTK)
/// ```
pub fn spke01(ctx: &mut SpiceContext, et: f64, record: &[f64], state: &mut [f64; 6]) {
    SPKE01(et, record, state, ctx.raw_context());
}

//$Procedure SPKE01 ( S/P Kernel, evaluate, type 1 )
pub fn SPKE01(ET: f64, RECORD: &[f64], STATE: &mut [f64], ctx: &mut Context) {
    let save = ctx.get_vars::<SaveVars>();
    let save = &mut *save.borrow_mut();

    let RECORD = DummyArray::new(RECORD, 1..);
    let mut STATE = DummyArrayMut::new(STATE, 1..=6);

    //
    // SPICELIB functions
    //

    //
    // Local variables
    //
    // The names below are original to the routine. They correspond
    // roughly to the original memos written by Fred Krogh to explain
    // how all this stuff really works.
    //

    //
    // Save everything between calls.
    //

    //
    // If the RETURN function is set, don't even bother with this.
    //
    if RETURN(ctx) {
        return;
    }

    //
    // Unpack the contents of the MDA array.
    //
    //    Name    Dimension  Description
    //    ------  ---------  -------------------------------
    //    TL              1  Final epoch of record
    //    G              15  Stepsize function vector
    //    REFPOS          3  Reference position vector
    //    REFVEL          3  Reference velocity vector
    //    DT         15,NTE  Modified divided difference arrays
    //    KQMAX1          1  Maximum integration order plus 1
    //    KQ            NTE  Integration order array
    //
    // For our purposes, NTE is always 3.
    //
    MOVED(RECORD.subarray(1), 1, std::slice::from_mut(&mut save.TL));
    MOVED(RECORD.subarray(2), 15, save.G.as_slice_mut());
    //
    // Collect the reference position and velocity.
    //
    save.REFPOS[1] = RECORD[17];
    save.REFVEL[1] = RECORD[18];

    save.REFPOS[2] = RECORD[19];
    save.REFVEL[2] = RECORD[20];

    save.REFPOS[3] = RECORD[21];
    save.REFVEL[3] = RECORD[22];

    MOVED(RECORD.subarray(23), 45, save.DT.as_slice_mut());

    save.KQMAX1 = (RECORD[68] as i32);
    save.KQ[1] = (RECORD[69] as i32);
    save.KQ[2] = (RECORD[70] as i32);
    save.KQ[3] = (RECORD[71] as i32);

    //
    // Next we set up for the computation of the various differences
    //
    save.DELTA = (ET - save.TL);
    save.TP = save.DELTA;
    save.MQ2 = (save.KQMAX1 - 2);
    save.KS = (save.KQMAX1 - 1);

    //
    // This is clearly collecting some kind of coefficients.
    // The problem is that we have no idea what they are...
    //
    // The G coefficients are supposed to be some kind of step size
    // vector.
    //
    // TP starts out as the delta t between the request time
    // and the time for which we last had a state in the MDL file.
    // We then change it from DELTA  by the components of the stepsize
    // vector G.
    //
    for J in 1..=save.MQ2 {
        save.FC[(J + 1)] = (save.TP / save.G[J]);
        save.WC[J] = (save.DELTA / save.G[J]);
        save.TP = (save.DELTA + save.G[J]);
    }

    //
    // Collect KQMAX1 reciprocals.
    //
    for J in 1..=save.KQMAX1 {
        save.W[J] = (1.0 / (J as f64));
    }

    //
    // Compute the W(K) terms needed for the position interpolation
    // (Note,  it is assumed throughout this routine that KS, which
    // starts out as KQMAX1-1 (the ``maximum integration'')
    // is at least 2.
    //
    save.JX = 0;
    save.KS1 = (save.KS - 1);

    while (save.KS >= 2) {
        save.JX = (save.JX + 1);

        for J in 1..=save.JX {
            save.W[(J + save.KS)] = ((save.FC[(J + 1)] * save.W[(J + save.KS1)])
                - (save.WC[J] * save.W[(J + save.KS)]));
        }

        save.KS = save.KS1;
        save.KS1 = (save.KS1 - 1);
    }

    //
    // Perform position interpolation: (Note that KS = 1 right now.
    // We don't know much more than that.)
    //
    for I in 1..=3 {
        save.KQQ = save.KQ[I];
        save.SUM = 0.0;

        for J in intrinsics::range(save.KQQ, 1, -1) {
            save.SUM = (save.SUM + (save.DT[[J, I]] * save.W[(J + save.KS)]));
        }

        STATE[I] = (save.REFPOS[I] + (save.DELTA * (save.REFVEL[I] + (save.DELTA * save.SUM))));
    }

    //
    // Again we need to compute the W(K) coefficients that are
    // going to be used in the velocity interpolation.
    // (Note, at this point, KS = 1, KS1 = 0.)
    //
    for J in 1..=save.JX {
        save.W[(J + save.KS)] =
            ((save.FC[(J + 1)] * save.W[(J + save.KS1)]) - (save.WC[J] * save.W[(J + save.KS)]));
    }

    save.KS = (save.KS - 1);

    //
    // Perform velocity interpolation:
    //
    for I in 1..=3 {
        save.KQQ = save.KQ[I];
        save.SUM = 0.0;

        for J in intrinsics::range(save.KQQ, 1, -1) {
            save.SUM = (save.SUM + (save.DT[[J, I]] * save.W[(J + save.KS)]));
        }

        STATE[(I + 3)] = (save.REFVEL[I] + (save.DELTA * save.SUM));
    }

    //
    // That's all folks.  We don't know why we did anything, but
    // at least we can tell structurally what we did.
    //
}