modified_equinoctial_derivs Subroutine

    subroutine modified_equinoctial_derivs(mu,evec,scn,evecd)

    implicit none

    real(wp),intent(in)               :: mu    !! central body gravitational parameter (\(km^3/s^2\))
    real(wp),dimension(6),intent(in)  :: evec  !! modified equinoctial element vector
    real(wp),dimension(3),intent(in)  :: scn   !! Perturbation (in the RSW frame)
    real(wp),dimension(6),intent(out) :: evecd !! derivative of `evec`

    real(wp) :: p,f,g,h,k,L,c,s,n,sqrtpm,sl,cl,s2no2w,s2,w

    p      = evec(1)
    f      = evec(2)
    g      = evec(3)
    h      = evec(4)
    k      = evec(5)
    L      = evec(6)
    s      = scn(1)
    c      = scn(2)
    n      = scn(3)
    sqrtpm = sqrt(p/mu)
    sl     = sin(L)
    cl     = cos(L)
    s2     = one + h*h + k*k
    w      = one + f*cl + g*sl
    s2no2w = s2*n / (two*w)

    evecd(1) = (two*p*c/w) * sqrtpm
    evecd(2) = sqrtpm * ( s*sl + ((w+one)*cl+f)*c/w - g*n*(h*sl-k*cl)/w )
    evecd(3) = sqrtpm * (-s*cl + ((w+one)*sl+g)*c/w + f*n*(h*sl-k*cl)/w )
    evecd(4) = sqrtpm * s2no2w * cl
    evecd(5) = sqrtpm * s2no2w * sl
    evecd(6) = sqrt(mu*p)*(w/p)**2 + sqrtpm * ((h*sl-k*cl)*n)/w

    end subroutine modified_equinoctial_derivs