From dcb292706b174071b0efa7fd95d2d7db2dc9e3a2 Mon Sep 17 00:00:00 2001
From: Alexis Gamelin <alexis.gamelin@synchrotron-soleil.fr>
Date: Tue, 23 Jan 2024 10:16:37 +0100
Subject: [PATCH] Fix Resonator

Fix resonator wake function for negative time and update docstrings.
---
 mbtrack2/impedance/resonator.py | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/mbtrack2/impedance/resonator.py b/mbtrack2/impedance/resonator.py
index 640a52d..d8dc346 100644
--- a/mbtrack2/impedance/resonator.py
+++ b/mbtrack2/impedance/resonator.py
@@ -14,6 +14,9 @@ class Resonator(WakeField):
         """
         Resonator model WakeField element which computes the impedance and the 
         wake function in both longitudinal and transverse case.
+        
+        ! For transverse case, if Q < 2, there is a kick factor mismatch 
+        between time and frequency domain.
 
         Parameters
         ----------
@@ -22,7 +25,7 @@ class Resonator(WakeField):
         frequency : array of float
             Frequency points where the impedance will be evaluated in [Hz].
         Rs : float
-            Shunt impedance in [ohm].
+            Shunt impedance in [ohm] in longitudinal or [ohm/m] in transverse.
         fr : float
             Resonance frequency in [Hz].
         Q : float
@@ -94,6 +97,8 @@ class Resonator(WakeField):
                    (2 * self.Q_p)))
         if np.any(np.abs(t) < atol):
             wl[np.abs(t) < atol] = wl[np.abs(t) < atol] / 2
+        elif np.any(t < -atol):
+            wl[t < -atol] = 0
         return wl
 
     def long_impedance(self, f):
@@ -105,13 +110,16 @@ class Resonator(WakeField):
 
     def transverse_wake_function(self, t):
         if self.Q >= 0.5:
-            return (self.wr * self.Rs / self.Q_p *
+            wt = (self.wr * self.Rs / self.Q_p *
                     np.exp(-1 * t * self.wr / 2 / self.Q_p) *
                     np.sin(self.wr_p * t))
         else:
-            return (self.wr * self.Rs / self.Q_p *
+            wt = (self.wr * self.Rs / self.Q_p *
                     np.exp(-1 * t * self.wr / 2 / self.Q_p) *
                     np.sinh(self.wr_p * t))
+        if np.any(t < 0):
+            wt[t<0] = 0
+        return wt
 
 
 class PureInductive(WakeField):
-- 
GitLab