gravity_j2_j3_j4 Subroutine

    subroutine gravity_j2_j3_j4(r,mu,req,j2,j3,j4,acc)

    implicit none

    real(wp),dimension(3),intent(in)  :: r   !! satellite position vector [km]
    real(wp),intent(in)               :: mu  !! central body gravitational parameter [km^3/s^2]
    real(wp),intent(in)               :: req !! body equatorial radius [km]
    real(wp),intent(in)               :: j2  !! j2 coefficient
    real(wp),intent(in)               :: j3  !! j3 coefficient
    real(wp),intent(in)               :: j4  !! j4 coefficient
    real(wp),dimension(3),intent(out) :: acc !! gravity acceleration vector [km/s^2]

    real(wp) :: mmor3,reqor,reqor2,reqor3,reqor4,&
                rmag,rmag2,rmag3,rzor,rzor2,rzor3,rzor4,c,d

    rmag2 = dot_product(r,r)
    rmag  = sqrt(rmag2)

    if (rmag==zero) then

        write(output_unit,'(A)') 'Error in gravity_j2_j3_j4: spacecraft at center of body.'
        acc = zero

    else

        rmag3  = rmag*rmag2
        mmor3  = -mu/rmag3
        reqor  = req/rmag
        reqor2 = reqor*reqor
        reqor3 = reqor2*reqor
        reqor4 = reqor3*reqor
        rzor   = r(3)/rmag
        rzor2  = rzor*rzor
        rzor3  = rzor2*rzor
        rzor4  = rzor3*rzor

        c = mmor3 * (1.0_wp - 1.5_wp*J2*reqor2*(5.0_wp*rzor2-1.0_wp) + &
                     2.5_wp*J3*reqor3*(3.0_wp*rzor-7.0_wp*rzor3) - &
                     0.625_wp*J4*reqor4*(3.0_wp-42.0_wp*rzor2+63.0_wp*rzor4))
        d = mmor3 * (r(3) + 1.5_wp*J2*reqor2*(3.0_wp-5.0_wp*rzor2)*r(3) + &
                   0.5_wp*J3*reqor3*(30.0_wp*rzor*r(3)-35.0_wp*rzor3*r(3)-3.0_wp*rmag) - &
                   0.625_wp*J4*reqor4*(15.0_wp-70.0_wp*rzor2+63.0_wp*rzor4)*r(3))

        acc(1) = c * r(1)
        acc(2) = c * r(2)
        acc(3) = d

    end if

    end subroutine gravity_J2_J3_J4