lambert_test Subroutine

    subroutine lambert_test()

    use gooding_module,    only: pv3els
    use random_module,     only: get_random_number

    implicit none

    real(wp),parameter :: tol = 1.0e-11_wp

    real(wp),dimension(:,:),allocatable :: izzo_v1,izzo_v2
    real(wp),dimension(:,:),allocatable :: gooding_v1,gooding_v2
    real(wp),dimension(6)   :: rv1,rv2,pv,e
    integer                 :: imeth,icases,i,multi_revs,iway,n_cases
    real(wp)                :: tof,err1,err2
    logical                 :: long_way,status_ok
    character(len=100)      :: str
    real                    :: tstart, tend
    real(wp),dimension(3)   :: arora_v1,arora_v2

    real(wp),dimension(6),parameter :: rv1_base = &
                                        [1e6_wp,1e6_wp,1e6_wp,10.0_wp,10.0_wp,10.0_wp]
    real(wp),dimension(6),parameter :: rv2_base = &
                                        [1e6_wp,1e6_wp,1e6_wp,10.0_wp,10.0_wp,10.0_wp]
    real(wp),parameter :: tof_base = 86400.0_wp  !sec
    real(wp),parameter :: mu = 3.986004362330e+05_wp

    !real(wp),parameter :: auora_tol      = 1.0e-14_wp   !  122946 cases/sec
    real(wp),parameter :: auora_tol       = 1.0e-9_wp    !! 203025 cases/sec
    integer,parameter  :: max_iterations = 100
    logical,parameter  :: shortperiod    = .true.   !! "solution 1" for mult-rev case.
    real(wp),parameter :: lu             = 1.0e5_wp

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' lambert_test'
    write(*,*) '---------------'
    write(*,*) ''

    write(*,*) ''
    write(*,*) '   Test 1'
    write(*,*) ''

    n_cases = 10
    multi_revs = 0  !1

    !reseed the random number generator:
    call random_seed()

    do icases=1,n_cases

        write(*,*) ''
        do i=1,6
            rv1(i) = get_random_number(-rv1_base(i),rv1_base(i))
        end do
        do i=1,6
            rv2(i) = get_random_number(-rv2_base(i),rv2_base(i))
        end do
        tof = get_random_number(1000.0_wp, tof_base)

        do iway=1,2    !short and long way

            long_way = (iway==1)

            do imeth=1,3    ![izzo, gooding, auora]

                !if (icases==1 .and. imeth==1 .and. iway==1) &
                !        write(*,'(*(A30,1X))') &
                !        'case',&
                !        'x1','y1','z1','vx1','vy1','vz1',&
                !        'x2','y2','z2','vx2','vy2','vz2','tof'
                !if (imeth==1) write(*,'(I30,1X,*(F30.6,1X))') icases, rv1, rv2, tof

                select case (imeth)
                case(1)
                    call solve_lambert_izzo(   rv1(1:3),rv2(1:3),tof,mu,long_way,&
                                                multi_revs,izzo_v1,izzo_v2,&
                                                status_ok)
                case(2)
                    call solve_lambert_gooding(rv1(1:3),rv2(1:3),tof,mu,long_way,&
                                                multi_revs,gooding_v1,gooding_v2,&
                                                status_ok)
                case(3)
                    call solve_lambert_arorarussell(rv1(1:3),rv2(1:3),tof,mu,lu,multi_revs,long_way,&
                                                            shortperiod,auora_tol,&
                                                            max_iterations,arora_v1,arora_v2)
                end select

            end do

            !results:
            if (size(izzo_v1,2)==size(gooding_v1,2)) then

                do i = 1, size(izzo_v1,2)

                    !orb. elements of transfer orbit:
                    pv = [rv1(1:3),gooding_v1(:,i)]
                    call pv3els(mu, pv, e)

                    err1 = norm2( izzo_v1(:,i) - gooding_v1(:,i) )
                    err2 = norm2( izzo_v2(:,i) - gooding_v2(:,i) )

                    if (err1>tol) then
                        str = '*****IZZO-GOODING ERROR*****'
                    else
                        str = '     IZZO-GOODING v1'
                    end if
                    write(*,'(I5,1X,E25.16,1X,E25.16,1X,E25.16,1X,A)')&
                                 icases, e(1:2), err1, trim(str)

                    if (err2>tol) then
                        str = '*****IZZO-GOODING ERROR*****'
                    else
                        str = '     IZZO-GOODING v2'
                    end if
                    write(*,'(I5,1X,E25.16,1X,E25.16,1X,E25.16,1X,A)')&
                                 icases, e(1:2), err2, trim(str)

                     err1 = norm2( arora_v1 - gooding_v1(:,i) )
                     err2 = norm2( arora_v2 - gooding_v2(:,i) )

                     if (err1>tol) then
                         str = '*****ARORA-GOODING ERROR*****'
                     else
                         str = '     ARORA-GOODING v1'
                     end if
                     write(*,'(I5,1X,E25.16,1X,E25.16,1X,E25.16,1X,A)')&
                                  icases, e(1:2), err1, trim(str)

                     if (err2>tol) then
                         str = '*****ARORA-GOODING ERROR*****'
                     else
                         str = '     ARORA-GOODING v2'
                     end if
                     write(*,'(I5,1X,E25.16,1X,E25.16,1X,E25.16,1X,A)')&
                                  icases, e(1:2), err2, trim(str)

                end do

            else
                write(*,*) icases, 'Error: arrays sizes are not the same'
                stop
            end if

        end do

    end do

    write(*,*) ''
    write(*,*) '   Test 2: Speed test'
    write(*,*) ''

    n_cases = 1000000

    do imeth=1,3    ![izzo, gooding, auora]

        !reseed the random number generator:
        call random_seed()

        call cpu_time(tstart)

        do icases=1,n_cases

            do i=1,6
                rv1(i) = get_random_number(-rv1_base(i),rv1_base(i))
            end do
            do i=1,6
                rv2(i) = get_random_number(-rv2_base(i),rv2_base(i))
            end do
            tof = get_random_number(1000.0_wp, tof_base)

            do iway=1,2    !short and long way

                long_way = (iway==1)

                select case (imeth)
                case(1)
                    call solve_lambert_izzo(   rv1(1:3),rv2(1:3),tof,mu,long_way,&
                                                multi_revs,izzo_v1,izzo_v2,&
                                                status_ok)
                case(2)
                    call solve_lambert_gooding(rv1(1:3),rv2(1:3),tof,mu,long_way,&
                                                multi_revs,gooding_v1,gooding_v2,&
                                                status_ok)
                case(3)
                    call solve_lambert_arorarussell(rv1(1:3),rv2(1:3),tof,mu,lu,multi_revs,long_way,&
                                                            shortperiod,auora_tol,&
                                                            max_iterations,arora_v1,arora_v2)
                end select

            end do

        end do

        call cpu_time(tend)
        select case (imeth)
        case(1)
            write(*,*) 'run time for Izzo   : ', tend-tstart, ' sec.  ', dble(n_cases)/(tend-tstart),' cases/sec'
        case(2)
            write(*,*) 'run time for Gooding: ', tend-tstart, ' sec.  ', dble(n_cases)/(tend-tstart),' cases/sec'
        case(3)
            write(*,*) 'run time for Arora  : ', tend-tstart, ' sec.  ', dble(n_cases)/(tend-tstart),' cases/sec'
        end select
        write(*,*) ''

    end do

    end subroutine lambert_test