subroutine get_commutator_coefs(qtype, nnodes, dt, coefs)
integer, intent(in) :: qtype, nnodes
real(pfdp), intent(in) :: dt
real(pfdp), intent(inout) :: coefs(:,:,:)
! coefs has the structure coefs(coefs, magnus_order, node)
! for a given node, pass subroutines the coefs for a magnus order, then
! loop over coefs
if (qtype == 1) then
! we're talking Lobatto nodes, where nnodes=3 includes, t0, t_1/2, tn
! need some way to differentiate whether you want full collocation or not
! coefs(1:3, 1, 1) = dt**2 * [real(8)::11/480., -1/480., 1/480.]
! coefs(1:3, 1, 2) = dt**2 * [real(8)::1/15., 1/60., 1/15.]
coefs(1, 1, 1) = -1/48.d0 * dt**2
coefs(2, 1, 2) = -1/12.d0 * dt**2
elseif (qtype == 5) then
coefs(1:3, 1, 1) = 1.d-3 * [real(8)::-0.708256232441739d0, 0.201427439334681d0, -0.002608155816283d0]
coefs(1:3, 1, 2) = [real(8)::-0.035291589565775d0, 0.004482619613666d0, -0.000569367343553d0]
coefs(1:3, 1, 3) = [real(8)::-0.078891497044705d0, -0.018131905893999d0, -0.035152700676886d0]
coefs(1:3, 1, 4) = [real(8)::-0.071721913818656d0, -0.035860956909328d0, -0.071721913818656d0]
coefs(:,1,:) = dt**2 * coefs(:,1,:)
coefs(:, 2, 1) = &
[real(8)::1.466782892818107d-6, -2.546845448743404d-6, 7.18855795894042d-7, &
-3.065370250683271d-7, 6.962336322868984d-7, -1.96845581200288d-7, &
-2.262216360714434d-8, -2.72797194008496d-9, 8.54843541920492d-10]
coefs(:, 2, 2) = &
[real(8) ::0.001040114336531742d0, -0.001714330280871491d0, 0.0001980882752518163d0, &
-0.00006910549596945875d0, 0.0002905401601450182d0, -0.00003465884693947625d0, &
0.0000924518848932026d0, 0.0000125950571649574d0, -2.4709074423913880d-6]
coefs(:, 2, 3) = &
[real(8)::0.004148295975360902d0, -0.006387421893168941d0, -0.003594231910817328d0, &
0.000997378110327084d0, 0.0001241530237557625d0, -0.0003805975423160699d0, &
0.003718384934573079d0, 0.001693514295056844d0, -0.001060408584538103d0]
coefs(:, 2, 4) = &
[real(8)::0.003453850676072909d0, -0.005584950029394391d0, -0.007128159905937654d0, &
0.001653439153439147d0, 0.0d0, -0.001653439153439143d0, &
0.007128159905937675d0, 0.005584950029394475d0, -0.003453850676072897d0]
coefs(:,2,:) = dt**3 * coefs(:,2,:)
coefs(1, 3, 4) = dt**4 / 60.d0
else
stop 'oh no! unsupported qtype'
endif
end subroutine get_commutator_coefs