Tweaked crater exterior function

MintonGroup · Dec 6, 2016 · 81bcac6 · 81bcac6
1 parent dca7b04
commit 81bcac6
Showing 1 changed file with 25 additions and 18 deletions.
diff --git a/src/crater/crater_form_exterior_rootfind.f90 b/src/crater/crater_form_exterior_rootfind.f90
@@ -43,34 +43,41 @@ subroutine crater_form_exterior_rootfind(user,surf,crater,domain,deltaMtot)
    real(DP) :: resultat,AA,BB,CC,DD,EE,FA,FB,FC,Tol1,PP,QQ,RR,SS,xm
    integer(I4B) :: i,j, ev,br
    integer(I4B),parameter :: NTRY = 50
-   integer(I4B),parameter :: NBRACKET = 5
-   real(DP),parameter :: FIRSTFACTOR = 1.6_DP
+   integer(I4B),parameter :: NBRACKET = 10
+   real(DP),parameter :: FIRSTFACTOR = 3.2_DP
    real(DP) :: factor
    real(DP) :: f1,f2
 
    ! Executable code
    if (abs(deltaMtot) < VSMALL) return
    if (deltaMtot > 0._DP) return
    lastloop = .false.
-   factor=FIRSTFACTOR
+   factor = FIRSTFACTOR
    startrd = RIMDROP
-   x1 = startrd
-   x2 = startrd / factor
+
 
    ! First bracket the root
-   f1 = crater_form_exterior_func(user,surf,crater,domain,x1,deltaMtot,lastloop)
-   f2 = crater_form_exterior_func(user,surf,crater,domain,x2,deltaMtot,lastloop)
-   do j = 1, NTRY
-      if (f1 * f2 < 0._DP) exit
-      if (abs(f1) < abs(f2)) then
-         x1 = x1 + factor * (x1 - x2)
-         f1 = crater_form_exterior_func(user,surf,crater,domain,x1,deltaMtot,lastloop)
-      else
-         x2 = x2 + factor * (x2 - x1)
-         f2 = crater_form_exterior_func(user,surf,crater,domain,x2,deltaMtot,lastloop)
-      end if
-
-   end do
+   bracket: do br=1,NBRACKET
+      factor = factor / 2
+      x1 = startrd
+      x2 = startrd / factor
+      f1 = crater_form_exterior_func(user,surf,crater,domain,x1,deltaMtot,lastloop)
+      f2 = crater_form_exterior_func(user,surf,crater,domain,x2,deltaMtot,lastloop)
+      do j = 1, NTRY
+         if (f1 * f2 < 0._DP) exit bracket
+         if (abs(f1) < abs(f2)) then
+            x1 = x1 + factor * (x1 - x2)
+            f1 = crater_form_exterior_func(user,surf,crater,domain,x1,deltaMtot,lastloop)
+         else
+            x2 = x2 + factor * (x2 - x1)
+            f2 = crater_form_exterior_func(user,surf,crater,domain,x2,deltaMtot,lastloop)
+         end if
+         if (abs(f1) > huge(f1).or.(abs(f2) > huge(f2))) cycle bracket
+         !write(*,*)
+         !write(*,*) j,factor,x1,x2,f1,f2,crater%fcrat
+         !read(*,*)
+      end do
+   end do bracket
    !!if (f1 * f2 >= 0._DP) then
    !   write(*,*) crater%fcrat
    !   write(*,*) deltaMtot