From 05be0276f05670d36e1b187062bc0904cec91463 Mon Sep 17 00:00:00 2001
From: ydamour <yann.damour@hotmail.fr>
Date: Thu, 3 Nov 2022 14:30:30 +0100
Subject: [PATCH] 2 by 2 diag

---
 src/utils/util.irp.f | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f
index 84593031..f295aea7 100644
--- a/src/utils/util.irp.f
+++ b/src/utils/util.irp.f
@@ -411,3 +411,28 @@ subroutine lowercase(txt,n)
   enddo
 end
 
+subroutine v2_over_x(v,x,res)
+
+  !BEGIN_DOC
+  ! Two by two diagonalization to avoid the divergence in v^2/x when x goes to 0
+  !END_DOC
+
+  implicit none
+
+  double precision, intent(in)  :: v, x
+  double precision, intent(out) :: res
+
+  double precision :: delta_E, tmp, val
+
+  res = 0d0
+  delta_E = x
+  if (v == 0.d0) return
+
+  val = 2d0 * v
+  tmp = dsqrt(delta_E * delta_E + val * val)
+  if (delta_E < 0.d0) then
+      tmp = -tmp
+  endif
+  res = 0.5d0 * (tmp - delta_E)
+
+end