From d01d3084ffc235a7efff8c41e23b97d148bc9ba4 Mon Sep 17 00:00:00 2001 From: Austin Michael Blevins Date: Wed, 25 Jan 2023 14:36:14 -0500 Subject: [PATCH] Fixed the first quartic bug; now working on one when ri=0 --- src/regolith/regolith_quartic_func.f90 | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/regolith/regolith_quartic_func.f90 b/src/regolith/regolith_quartic_func.f90 index a91be50e..ae3ead39 100755 --- a/src/regolith/regolith_quartic_func.f90 +++ b/src/regolith/regolith_quartic_func.f90 @@ -34,27 +34,29 @@ function regolith_quartic_func(rpj,r) result(z) real(DP) :: e,delta0,delta1,delta,S,T,x12,x34,costheta real(DP),dimension(4) :: x integer(I4B) :: i + costheta = 0.0 -e = (rpj/r)**2 - 1.0 -delta0 = 12.0*(rpj/r)**2 -delta1 = 108.0*(rpj/r)**2 +e = (rpj/r)**2 - 1.0_DP +delta0 = 12.0_DP*(rpj/r)**2 +delta1 = 108.0_DP*(rpj/r)**2 delta = sqrt(delta1**2 - 4.0*delta0**3) -T = ((delta1 + delta)/2.0)**(1.0/3.0) -S = 0.5*sqrt( 1.0 + (T + delta0/T)/3.0 ) -x12 = sqrt(-4.0*S*S + 3.0 + 1.0/S) +T = ((delta1 + delta)/2.0_DP)**(1.0_DP/3.0_DP) +S = 0.5*sqrt( 1.0_DP + (T + delta0/T)/3.0_DP ) +x12 = sqrt(-4.0_DP*S*S + 3.0_DP + 1.0_DP/S) x(1) = 0.5 - S + 0.5*x12 x(2) = 0.5 - S - 0.5*x12 -x34 = sqrt(-4.0*S*S + 3.0 - 1.0/S) -x(3) = 0.5 + S + 0.5*x34 -x(4) = 0.5 + S - 0.5*x34 +x34 = max(-4*S**2 + 3.0_DP - 1.0_DP/S,0.0_DP) +x34 = sqrt(x34) +x(3) = 0.5_DP + S + 0.5_DP*x34 +x(4) = 0.5_DP + S - 0.5_DP*x34 do i = 1,4 - if (x(i) <= 1.0 .and. x(i) >= 0.0) then + if (x(i) <= 1.0_DP .and. x(i) >= 0.0_DP) then costheta = x(i) end if end do !if (present(theta)) then - z = costheta +z = costheta !else ! z = r * (1.0 - costheta) * costheta !end if