Using DIVA to simultaneously integrate the variational equations two-body problem scenario

[OrbitalMechanic]

I'm experimenting with the DIVA numerical integrator from MATH77 to simultaneously integrate the variational equations for a two-body problem scenario.

Now as I understand it the 3x12 matrix that represents the state elements for the variational equations is given below.

The questions that come to mind are:

1. What is the storage order for the partial derivatives values stored in the array containing the state (position and velocity) estimates?

2. How should the vector equation of motion and the variational equation values be stored in the array that contains the right-hand-sides and error tolerances?

3. How should the control flags, KORD, IOPT be set up to simultaneously numerically integrate the vector equation of motion with the variational equations?

4. Does anyone have a sample code using DIVA to simultaneously numerically integrate the vector equation of motion with the variational equations with output?

Any suggestions?

[jacobwilliams]

In general, if you're integrating a matrix, you just need to pack it into the vector like so:

    real(wp),dimension(6)   :: rv     
    real(wp),dimension(42)  :: x     
    real(wp),dimension(6,6) :: phi

    x(1:6)  = rv
    x(7:42) = pack (phi, mask=.true.)