Merge pull request #50 from aoeftiger/feature/sc-linedensity-interpolation

Feature/sc linedensity interpolation

Author: rdemaria

examples/spacecharge/001_prepare_spacecharge_lattice_and_particles.py

@@ -17,12 +17,14 @@
 
 if particle_type == "protons":
 
-    # For bunched
+    # Space charge:
     number_of_particles = 2e11
+
+    # For bunched
     bunchlength_rms = 0.22
 
     # For coasting
-    line_density = 2e11 / 0.5
+    circumference = 0.5
 
     mass = pysixtrack.Particles.pmass
     p0c = 25.92e9
@@ -37,20 +39,22 @@
 
 if particle_type == "ions":
 
-    # For bunched
+    # Space charge:
     number_of_particles = 3.5e8
+
+    # For bunched
     bunchlength_rms = 0.22
 
     # For coasting
-    line_density = 3.5e8 / 0.5
+    circumference = 0.5
 
     mass = 193.7e9
     p0c = 1402.406299e9
     charge_state = 82.0
-    neps_x = 1.63e-6 
-    neps_y = 0.86e-6 
+    neps_x = 1.63e-6
+    neps_y = 0.86e-6
     delta_rms = 1.0e-3
-    V_RF_MV = 3 
+    V_RF_MV = 3
     lag_RF_deg = 0.0
     n_SCkicks = 250
     length_fuzzy = 1.5
@@ -90,11 +94,11 @@
 # Check consistency
 if sc_mode == "Bunched":
     sc_elements, sc_names = line.get_elements_of_type(
-        pysixtrack.elements.SpaceChargeBunched
+        pysixtrack.elements.SCQGaussProfile
     )
 elif sc_mode == "Coasting":
     sc_elements, sc_names = line.get_elements_of_type(
-        pysixtrack.elements.SpaceChargeCoasting
+        pysixtrack.elements.SCCoasting
     )
 else:
     raise ValueError("mode not understood")
@@ -148,7 +152,8 @@
         sc_lengths,
         sc_twdata,
         betagamma,
-        line_density,
+        number_of_particles,
+        circumference,
         delta_rms,
         neps_x,
         neps_y,
@@ -177,7 +182,7 @@
 with open("particle_on_CO.pkl", "rb") as fid:
     part_on_CO = pysixtrack.Particles.from_dict(pickle.load(fid))
 
-part = part_on_CO.copy()  
+part = part_on_CO.copy()
 
 # get beta functions from twiss table
 with open("twiss_at_start.pkl", "rb") as fid:

pysixtrack/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.0.5"
+__version__ = "0.0.6"
 
 from . import base_classes
 from . import elements

pysixtrack/be_beamfields/qgauss.py

@@ -0,0 +1,75 @@
+from pysixtrack.mathlibs import MathlibDefault
+
+
+class QGauss(object):
+    @staticmethod
+    def calc_cq(q, mathlib=MathlibDefault, EPS=1e-6):
+        assert q < 3
+        cq = mathlib.sqrt(mathlib.pi)
+        if q >= (1 + EPS):
+            cq *= mathlib.gamma((3 - q) / (2 * q - 2))
+            cq /= mathlib.sqrt((q - 1)) * mathlib.gamma(1 / (q - 1))
+        elif q <= (1 - EPS):
+            cq *= 2 * mathlib.gamma(1 / (1 - q))
+            cq /= (
+                (3 - q)
+                * mathlib.sqrt(1 - q)
+                * mathlib.gamma((3 - q) / (2 - 2 * q))
+            )
+        return cq
+
+    @staticmethod
+    def sqrt_beta(sigma, mathlib=MathlibDefault):
+        assert sigma > 0
+        return 1 / (mathlib.sqrt(2) * sigma)
+
+    @staticmethod
+    def exp_q(x, q, mathlib=MathlibDefault, EPS=1e-6):
+        assert q < 3
+        if mathlib.abs(1 - q) > EPS:
+            u_plus = 1 + x * (1 - q)
+            if u_plus < 0:
+                u_plus = 0
+            return mathlib.pow(u_plus, 1 / (1 - q))
+        else:
+            return mathlib.exp(x)
+
+    def __init__(self, q=1.0, mathlib=MathlibDefault, cq_eps=1e-6):
+        assert q < 3
+        self._m = mathlib
+        self._q = q
+        self._cq = QGauss.calc_cq(q, mathlib=mathlib, EPS=cq_eps)
+
+    @property
+    def q(self):
+        return self._q
+
+    @q.setter
+    def q(self, q_value):
+        assert q_value < 3
+        self._q = q_value
+        self._cq = QGauss.calc_cq(self._q, self._m)
+
+    @property
+    def cq(self):
+        return self._cq
+
+    def min_support(self, sqrt_beta):
+        assert self._q < 3
+        assert sqrt_beta > 0
+        if self._q >= 1:
+            return -1e10
+        else:
+            return -1 / self._m.sqrt(sqrt_beta * sqrt_beta * (1 - self._q))
+
+    def max_support(self, sqrt_beta):
+        return -(self.min_support(sqrt_beta))
+
+    def eval(self, x, sqrt_beta, mu=0.0):
+        assert self._m.abs(self._cq) > 0
+        assert self._q < 3
+        assert sqrt_beta > 0
+        factor = sqrt_beta / self._cq
+        arg = sqrt_beta * sqrt_beta
+        arg *= (x - mu) * (x - mu)
+        return factor * QGauss.exp_q(-arg, self._q, self._m)