From 3726ec07c709d732cc43df08e0043baa6fdc155f Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Sat, 11 Feb 2017 00:52:48 +0100 Subject: [PATCH] protect against sin(theta_b) being zero --- src/mssm_twoloophiggs_impl.f | 71 +++++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/src/mssm_twoloophiggs_impl.f b/src/mssm_twoloophiggs_impl.f index a3e79a877..3becfd395 100644 --- a/src/mssm_twoloophiggs_impl.f +++ b/src/mssm_twoloophiggs_impl.f @@ -41,6 +41,8 @@ subroutine DDSHiggs(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, $ ,FAp double precision, parameter :: eps_st = 1d-5 double precision, parameter :: eps_t1 = 1d-5 + double precision, parameter :: eps_sb = 1d-5 + double precision, parameter :: eps_b1 = 1d-5 pi = 3.14159265897d0 @@ -67,6 +69,27 @@ subroutine DDSHiggs(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, endif endif c end of addition by ALEX + +c ADDED by ALEX: guards against NANs when sin theta is zero! + if (dabs(sb).lt.eps_sb) then + if (sb.ge.0.0d0) then + sb = eps_sb + else + sb = -eps_sb + endif + endif + cb = sqrt(1.d0-sb*sb) +c end of addition by ALEX + +c ADDED by ALEX: guards against NANs when B1 == B2 + if (dabs((B1-B2)/B1).lt.eps_b1) then + if (B1.gt.B2) then + B1 = B2 / (1d0 - eps_b1) + else + B1 = B2 / (1d0 + eps_b1) + endif + endif +c end of addition by ALEX s2t = 2d0*ct*st s2b = 2d0*cb*sb @@ -129,6 +152,8 @@ subroutine DDSodd(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, $ ,FAp double precision, parameter :: eps_st = 1d-5 double precision, parameter :: eps_t1 = 1d-5 + double precision, parameter :: eps_sb = 1d-5 + double precision, parameter :: eps_b1 = 1d-5 pi = 3.14159265897d0 @@ -156,6 +181,27 @@ subroutine DDSodd(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, endif c end of addition by ALEX +c ADDED by ALEX: guards against NANs when sin theta is zero! + if (dabs(sb).lt.eps_sb) then + if (sb.ge.0.0d0) then + sb = eps_sb + else + sb = -eps_sb + endif + endif + cb = sqrt(1.d0-sb*sb) +c end of addition by ALEX + +c ADDED by ALEX: guards against NANs when B1 == B2 + if (dabs((B1-B2)/B1).lt.eps_b1) then + if (B1.gt.B2) then + B1 = B2 / (1d0 - eps_b1) + else + B1 = B2 / (1d0 + eps_b1) + endif + endif +c end of addition by ALEX + s2t = 2d0*ct*st s2b = 2d0*cb*sb c2t = ct**2 - st**2 @@ -202,6 +248,8 @@ subroutine DDStad(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, double precision v1,v2,S1,S2 double precision, parameter :: eps_st = 1d-10 double precision, parameter :: eps_t1 = 1d-8 + double precision, parameter :: eps_sb = 1d-10 + double precision, parameter :: eps_b1 = 1d-8 pi = 3.14159265897d0 @@ -228,6 +276,27 @@ subroutine DDStad(t,b,A0,T1,T2,B1,B2,st,ct,sb,cb,q,mu,tanb,vv, endif endif c end of addition by ALEX + +c ADDED by ALEX: guards against NANs when sin theta is zero! + if (dabs(sb).lt.eps_sb) then + if (sb.ge.0.0d0) then + sb = eps_sb + else + sb = -eps_sb + endif + endif + cb = sqrt(1.d0-sb*sb) +c end of addition by ALEX + +c ADDED by ALEX: guards against NANs when B1 == B2 + if (dabs((B1-B2)/B1).lt.eps_b1) then + if (B1.gt.B2) then + B1 = B2 / (1d0 - eps_b1) + else + B1 = B2 / (1d0 + eps_b1) + endif + endif +c end of addition by ALEX s2t = 2d0*ct*st s2b = 2d0*cb*sb @@ -308,7 +377,7 @@ subroutine makefuncs(t,b,A0,T1,T2,B1,B2,s2t,c2t,s2b,c2b, DB1c2b = DT1c2t DB2c2b = DT2c2t Dtc2b = Dbc2t - + call makederiv(t,b,A0,T1,T2,B1,B2,s2t,c2t,s2b,c2b,q,mu,vv,tanb) F1t = Dtt + DT1T1 + DT2T2 + 2d0*(DT1t + DT2t + DT1T2)