crtbp_test Subroutine

public subroutine crtbp_test()

Uses

  • proc~~crtbp_test~~UsesGraph proc~crtbp_test crtbp_module::crtbp_test module~celestial_body_module celestial_body_module proc~crtbp_test->module~celestial_body_module module~base_class_module base_class_module module~celestial_body_module->module~base_class_module module~kind_module kind_module module~celestial_body_module->module~kind_module module~numbers_module numbers_module module~celestial_body_module->module~numbers_module iso_fortran_env iso_fortran_env module~kind_module->iso_fortran_env module~numbers_module->module~kind_module

Unit tests for CRTBP routines.

Arguments

None

Calls

proc~~crtbp_test~~CallsGraph proc~crtbp_test crtbp_module::crtbp_test proc~compute_crtpb_parameter crtbp_module::compute_crtpb_parameter proc~crtbp_test->proc~compute_crtpb_parameter proc~compute_jacobi_constant crtbp_module::compute_jacobi_constant proc~crtbp_test->proc~compute_jacobi_constant proc~compute_libration_points crtbp_module::compute_libration_points proc~crtbp_test->proc~compute_libration_points proc~compute_libration_points_v2 crtbp_module::compute_libration_points_v2 proc~crtbp_test->proc~compute_libration_points_v2 proc~crtbp_derivs crtbp_module::crtbp_derivs proc~crtbp_test->proc~crtbp_derivs proc~crtbp_derivs_with_stm crtbp_module::crtbp_derivs_with_stm proc~crtbp_test->proc~crtbp_derivs_with_stm proc~cube_root math_module::cube_root proc~compute_libration_points->proc~cube_root proc~compute_libration_points_v2->proc~compute_libration_points proc~compute_libration_points_v2->proc~cube_root

Source Code

    subroutine crtbp_test()

    use celestial_body_module

    implicit none

    real(wp),parameter :: mu_earth = body_earth%mu !! \( \mu_{Earth} ~ (\mathrm{km}^3/\mathrm{s}^2) \)
    real(wp),parameter :: mu_moon  = body_moon%mu  !! \( \mu_{Moon}  ~ (\mathrm{km}^3/\mathrm{s}^2) \)
    real(wp),parameter :: mu_sun   = body_sun%mu   !! \( \mu_{Sun}   ~ (\mathrm{km}^3/\mathrm{s}^2) \)

    !< sample state (normalized)
    !< see: [Celestial Mechanics Notes Set 4: The Circular Restricted
    !< Three Body Problem](http://cosweb1.fau.edu/~jmirelesjames/hw4Notes.pdf), p.40.
    real(wp),dimension(6),parameter :: x = [  0.30910452642073_wp, &
                                              0.07738174525518_wp, &
                                              0.0_wp,              &
                                             -0.72560796964234_wp, &
                                              1.55464233412773_wp, &
                                              0.0_wp               ]

    integer                 :: i       !! counter
    real(wp)                :: mu      !! CRTPB parameter
    real(wp)                :: mu1     !! primary body mu
    real(wp)                :: mu2     !! secondary body mu
    real(wp)                :: c       !! Jacobi constant
    real(wp),dimension(6)   :: xd      !! derivative vector: state
    real(wp),dimension(42)  :: x_phi   !! initial state + phi (identity)
    real(wp),dimension(42)  :: x_phi_d !! derivative vector: state + phi
    real(wp),dimension(6,6) :: eye     !! 6x6 identity matrix
    real(wp)                :: r1      !! L1 x coordinate (normalized)
    real(wp)                :: r2      !! L2 x coordinate (normalized)
    real(wp)                :: r3      !! L3 x coordinate (normalized)
    real(wp),dimension(2)   :: r4      !! L4 x coordinate (normalized)
    real(wp),dimension(2)   :: r5      !! L5 x coordinate (normalized)

    !create an identity matrix for stm initial condition:
    eye = zero
    do i = 1, 6
        eye(i,i) = one
    end do
    x_phi = [x, pack(eye,mask=.true.)]

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' crtbp_test'
    write(*,*) '---------------'
    write(*,*) ''

    do i=1,3

        select case (i)
        case(1)
            mu1 = mu_earth
            mu2 = mu_moon
        case(2)
            mu1 = mu_earth
            mu2 = mu_earth
        case(3)
            mu1 = mu_earth + mu_moon/four
            mu2 = mu_earth
        end select

        write(*,*) ''
        mu = compute_crtpb_parameter(mu1,mu2)
        c = compute_jacobi_constant(mu,x)
        call crtbp_derivs(mu,x,xd)
        call crtbp_derivs_with_stm(mu,x_phi,x_phi_d)
        call compute_libration_points(mu,r1,r2,r3,r4,r5)

        write(*,'(A,1X,*(F30.16,1X))')  'mu:         ', mu
        write(*,'(A,1X,*(F30.16,1X))' ) 'L1 x:       ', r1
        write(*,'(A,1X,*(F30.16,1X))' ) 'L2 x:       ', r2
        write(*,'(A,1X,*(F30.16,1X))' ) 'L3 x:       ', r3
        write(*,'(A,1X,*(F30.16,1X))' ) 'L4 x:       ', r4
        write(*,'(A,1X,*(F30.16,1X))' ) 'L5 x:       ', r5
        write(*,'(A,1X,*(F30.16,1X))' ) 'x:          ', x
        write(*,'(A,1X,*(F30.16,1X))' ) 'c:          ', c
        write(*,'(A,1X,*(F30.16,1X))' ) 'xd:         ', xd
        write(*,'(A,1X,*(F30.16,1X))' ) 'x+phi:      ', x_phi
        write(*,'(A,1X,*(F30.16,1X))' ) 'xd+phi_dot: ', x_phi_d
        write(*,*) ''

        call compute_libration_points_v2(mu,r1,r2,r3,r4,r5)
        write(*,*) ''
        write(*,*) 'alternate formulation:'
        write(*,'(A,1X,*(F30.16,1X))' ) 'L1 x:       ', r1
        write(*,'(A,1X,*(F30.16,1X))' ) 'L2 x:       ', r2
        write(*,'(A,1X,*(F30.16,1X))' ) 'L3 x:       ', r3
        write(*,'(A,1X,*(F30.16,1X))' ) 'L4 x:       ', r4
        write(*,'(A,1X,*(F30.16,1X))' ) 'L5 x:       ', r5
        write(*,*) ''

    end do

    end subroutine crtbp_test