Bug fix for the case when finite arithmetic errors causes the test fo…

…r the final value of the independent variable to fail.

README.md

@@ -28,9 +28,10 @@ FLINT is a modern object-oriented fortran library that provides four adaptive st
 + 4 Adaptive-step (fixed step-size also supported) ERK integrators: DOP54, DOP853, Verner98R, Verner65E
 + Any other ERK method can be implemented by just including their coefficients
 + Dense output with delayed interpolation (integrate once, interpolate as many times)
-+ Multiple event-detection as well as finding location of events using root-finding  (Brent's algorithm) with static and dynamic event masking
++ Multiple events can be defined for event detection during integration
++ Events location finding using root-finding  (Brent's algorithm) with static and dynamic event masking
 + Ability to set a maximum delay (referred to here as event step-size) after which events are guauanteed to be detected
-+ Ability to restart the integration or change solution on the detection of events
++ Ability to restart the integration or change solution on the detection of events for handling discontinuities in the differential equations
 + Stiffness detection
 
 
@@ -118,7 +119,7 @@ See the test program files, test.f90 and DiffEq.f90, in tests folder for a compa
     end subroutine SampleEventTB  
 ```
 
-2. Initialize the differential equation and ERK class objects for using Runge-Kutta intgerators.
+2. Initialize the differential equation and ERK class objects for using Runge-Kutta integrators.
 
 ```fortran
         use FLINT
@@ -156,7 +157,7 @@ See the test program files, test.f90 and DiffEq.f90, in tests folder for a compa
     end if    
 ```
 
-4. Call the Interpolate function for computing solution on the desired grid of x values. The last parameter must be specified as True if user wants FLINT to keep the internal storage for calling Interpolate again. Otherwise, the internal storage is deleted and the user must intgerate the equations again before calling Interpolate.
+4. Call the Interpolate function for computing solution on the desired grid of x values. The last parameter must be specified as True if user wants FLINT to keep the internal storage for calling Interpolate again. Otherwise, the internal storage is deleted and the user must integrate the equations again before calling Interpolate.
 
 ```fortran      
     real(WP), dimension(:), allocatable :: Xarr1, Xarr2

src/ERK.f90

@@ -45,7 +45,7 @@ module ERK
     !> Stiffness test default number of steps. See Hairer's DOP853 code.
     integer, parameter :: STIFFTEST_STEPS = 1000
 
-    !< All the supported ERK Methods
+    !> All the supported ERK Methods
     enum, bind(C)
         enumerator :: ERK_DOP853 = 17   !< Hairer's DOP853
         enumerator :: ERK_DOP54       !< Dormand-Prince 5(4)
@@ -61,7 +61,7 @@ module ERK
     public :: ERK_DOP853, ERK_DOP54, ERK_VERNER98R, ERK_VERNER65E  !, ERK_VERNER98E, ERK_VERNER87R, ERK_VERNER87E 
     !public :: ERK_VERNER76R, ERK_VERNER76E
 
-    !< The ERK class for all the ERK methods available to the user. It inherits FLINT_class.
+    !> The ERK class for all the ERK methods available to the user. It inherits FLINT_class.
     type, extends(FLINT_class), public :: ERK_class
 
         private

src/ERKIntegrate.f90

@@ -685,9 +685,9 @@ module subroutine erk_int(me, X0, Y0, Xf, Yf, StepSz, UseConstStepSz, IntStepsOn
                     EventStates = reshape(EventData, [n+2, int(size(EventData)/(n+2))])
                 end if
 
-                if (X == Xf .AND. status == FLINT_SUCCESS) then
+                if (LAST_STEP .AND. status == FLINT_SUCCESS) then
                     ! we finished like we should
-                elseif (EventsOn .AND. status == FLINT_EVENT_TERM) then
+                elseif (EventsOn .AND. LAST_STEP .AND. status == FLINT_EVENT_TERM) then
                     ! One of the terminal event has triggered, do nothing
                 else if (TotalSteps >= MaxSteps .AND. status == FLINT_SUCCESS) then
                     ! maximum steps reached

src/FLINT.f90

@@ -138,7 +138,7 @@
 !!                  end subroutine SampleEventTB  
 !!              \endcode
 !!
-!!              + Initialize the differential equation and ERK class objects for using Runge-Kutta intgerators
+!!              + Initialize the differential equation and ERK class objects for using Runge-Kutta integrators
 !!
 !!              \code{.f90}
 !!                      use FLINT
@@ -179,7 +179,7 @@
 !!              + Call the Interpolate function for computing solution on the desired grid of x values.
 !!              The last parameter must be specified as True if user wants FLINT to keep the internal storage
 !!              for calling Interpolate again. Otherwise, the internal storage is deleted and the user
-!!              must intgerate the equations again before calling Interpolate.
+!!              must integrate the equations again before calling Interpolate.
 !!
 !!              \code{.f90}
 !!                  real(WP), dimension(:), allocatable :: Xarr1, Xarr2

src/FLINT_base.f90

@@ -94,7 +94,7 @@ module FLINT_base
         !> Restart the integration
         enumerator :: FLINT_EVENTACTION_RESTARTINT  = 4
 
-        !> Mask the event for future detection. It can be masked with 
+        !> Mask the event for future detection. It can be combined with 
         !! any of the other event actions.
         enumerator :: FLINT_EVENTACTION_MASK        = 128   
     end enum
@@ -109,12 +109,12 @@ module FLINT_base
         !> Return the solution in EventStates at the beginning of the integrator's step 
         !! (with natural or event step size) just after which the event function 
         !! changed sign.
-        enumerator :: FLINT_EVENTOPTION_STEPBEGIN  = 1   !< Terminate the integration
+        enumerator :: FLINT_EVENTOPTION_STEPBEGIN  = 1
 
         !> Return the solution in EventStates at the end of the integrator's step 
         !! (with natural or event step size) just before which the event function
         !! changed sign.        
-        enumerator :: FLINT_EVENTOPTION_STEPEND    = 2   !< Change the solution value
+        enumerator :: FLINT_EVENTOPTION_STEPEND    = 2
     end enum
 
 
@@ -285,14 +285,14 @@ subroutine EventFunc(me, X, Y, EvalEvents, Value, Direction, LocEvent, LocEventA
             !! + Direction=0: all sign changes are detected
             integer, dimension(:), intent(out) :: Direction
 
-            !> FLINT will only call EventFunc with this option after an event is located.
-            !! In that case, LocEventIndex will contain the index between 1 and m of the 
+            !> FLINT will only call EventFunc with this parameter after an event is located.
+            !! In that case, it will contain the index between 1 and m of the 
             !! located event. User must check whether this parameter is "present" before
             !! using it.
             integer, intent(in), optional :: LocEvent
 
             !> This parameter will be only present if an event has been located and
-            !! LocEventIndex contains a valid index of the located event. In that case,
+            !! LocEvent contains a valid index of the located event. In that case,
             !! user can specify a one of following actions that the event handler
             !! will take:
             !! + FLINT_EVENTACTION_CONTINUE: Continue the integration (default)

src/tests/test.f90

@@ -69,6 +69,7 @@ program TestFLINT
 
     type(ERK_class) erkvar
     type(CR3BPSys) :: CR3BPDEObject
+    type(TBSys) :: TBDEObject
 
     real(WP) :: x0, xf, xfval, rnum
     real(WP), dimension(6) :: y0, yf, y0r
@@ -85,9 +86,54 @@ program TestFLINT
     integer :: method
 
     real(WP) :: t0, tf
-
-    ! Arenstorf orbit: Earth-moon mass-ratio
-    real(wp), parameter :: EMmu  = 0.012277471_WP
+
+    real(WP), parameter :: Moon_GM = 4902.800066163796_WP
+    real(WP), parameter :: EMmu  = 0.012277471_WP
+
+    ! Two-Body Orbit
+    real(WP), parameter :: TOF0 = -43139.98958152457_WP
+    real(WP), parameter :: TOF1 = TOF0 + 43200.0_WP
+    real(WP), parameter, dimension(6) :: y0tb = [-1.15997748e+04_WP, &
+                                                2.11325321e+04_WP,&
+                                                -1.53669890e+04_WP,&
+                                                1.08821230e-01_WP, &
+                                                -2.35407098e-01_WP,&
+                                                3.36994476e-01_WP]
+
+    TBDEObject%GM = Moon_GM    
+    TBDEObject%n = 6
+    TBDEObject%m = 0
+
+
+    stepsz0 = 0.0E-3_WP    ! let FLINT compute the initial step-size    
+    stifftestval = 1  ! check for stiffness and stop integration if stiff
+
+
+    ! Two-Body propagation
+    call erkvar%Init(TBDEObject, MAX_STEPS, Method=ERK_DOP853, ATol=[tol*1.0e-3], RTol=[tol],&
+        InterpOn=.FALSE., EventsOn=.FALSE. )
+    if (erkvar%status == FLINT_SUCCESS) then       
+        y0r = y0tb
+        stiffstatus = stifftestval
+        stepsz = stepsz0
+        xfval = TOF1
+        call erkvar%Integrate(TOF0, y0r, xfval, yf, StepSz=stepsz, UseConstStepSz=CONST_STEPSZ, &
+            IntStepsOn=.TRUE.,Xint = Xint, Yint = Yint, &
+            EventStates=EventStates, EventMask = EvMask,StiffTest=stiffstatus)
+
+        if (erkvar%status == FLINT_SUCCESS .OR. erkvar%status == FLINT_EVENT_TERM &
+                .OR. erkvar%status == FLINT_ERROR_MAXSTEPS) then 
+            call erkvar%Info(stiffstatus, errstr, nAccept=naccpt, nReject=nrejct, nFCalls=fcalls)
+
+            print *, xfval, yf, fcalls, naccpt, nrejct
+        else
+            call erkvar%Info(stiffstatus, errstr)
+            print *, mname//': Integrate failed: ', erkvar%status, ':', errstr
+        end if
+    end if
+
+
+    ! Arenstorf orbit: Earth-moon system
     x0    = 0.0
     xf    = 17.0652165601579625588917206249_WP
     y0    = [0.994_WP, 0.0_WP, 0.0_WP, 0.0_WP, -2.00158510637908252240537862224_WP, 0.0_WP]  !! initial `y` value
@@ -102,8 +148,6 @@ program TestFLINT
     Xarr = [(x0+ipdx*itr,itr=0,(nIp-1))]
     allocate(Yarr(6,size(Xarr)))
 
-    stepsz0 = 0.0E-3_WP    ! let FLINT compute the initial step-size    
-    stifftestval = 1  ! check for stiffness and stop integration if stiff
 
     ! create results file
     open(unit=17,file= fname, status = 'replace')