ringparam always returns an emittance

atmat/atphysics/ParameterSummaryFunctions/atsummary.m

@@ -73,19 +73,18 @@
         dp=TD(1).ClosedOrbit(5);
 
         % For calculating the synchrotron integrals
-        [I1d,I2d,I3d,I4d,I5d,I6,~] = DipoleRadiation(ring,TD);
+        [I1d,I2d,I3d,I4d,I5d,I6,~] = ElementRadiation(ring,TD);
         [I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,TD);
         smm.integrals=[I1d+I1w,I2d+I2w,I3d+I3w,I4d+I4w,I5d+I5w,I6];
-
+        
         if is6d
             alphac=mcf(atdisable_6d(ring),dp);
             eloss=atgetU0(ring,'method','tracking');            % eV
-            offmom = (abs(rev_frequency*harm_number - freq) > 2);
         else
             alphac=mcf(ring,dp);
             eloss=1.0e9*Cgamma/2/pi*smm.e0.^4*smm.integrals(2); % eV
-            offmom = (abs(dp) > 1.e-4);
         end
+
         etac=1/gamma/gamma- alphac;
         smm.compactionFactor = alphac;
         smm.etac = etac;
@@ -100,19 +99,9 @@
         smm.damping(3) = 2 + smm.integrals(4)/smm.integrals(2);
 
         smm.radiation = 1.0e-9*eloss;   % GeV
-        if offmom
-            warning('AT:NoEmit','\n%s\n%s\n%s',...
-                'For off-momentum lattices, the equilibrium emittance',...
-                'cannot be computed because the contribution of quadrupoles is missing.',...
-                'To avoid this warning, use ">> warning(''off'',''AT:NoEmit'')"');
-            smm.naturalEnergySpread = NaN;
-            smm.naturalEmittance    = NaN;
-            alphac2=NaN;
-        else
-            smm.naturalEnergySpread = gamma*sqrt(Cq*smm.integrals(3)/(2*smm.integrals(2) + smm.integrals(4)));
-            smm.naturalEmittance    = Cq*gamma.^2*smm.integrals(5)/(smm.integrals(2)-smm.integrals(4));
-            alphac2 = smm.integrals(1)/smm.circumference;
-        end
+        smm.naturalEnergySpread = gamma*sqrt(Cq*smm.integrals(3)/(2*smm.integrals(2) + smm.integrals(4)));
+        smm.naturalEmittance    = Cq*gamma.^2*smm.integrals(5)/(smm.integrals(2)-smm.integrals(4));
+        alphac2 = smm.integrals(1)/smm.circumference;
 
         % Damping times
         cf=smm.radiation/2/smm.revTime/smm.e0;

atmat/atphysics/ParameterSummaryFunctions/ringpara.m

@@ -60,12 +60,7 @@
 
         beta_c = betar*cspeed;
         [ringdata,lindata]=atlinopt6(ring,1:length(ring)+1,varargin{:},'get_chrom');
-        if is6d
-            offmom = (abs(frev*harm - freq_rf) > 2);
-        else
-            offmom = (abs(lindata(1).ClosedOrbit(5)) > 1.e-4);
-        end
-        [I1d,I2d,I3d,I4d,I5d,~,Iv] = DipoleRadiation(ring,lindata);
+        [I1d,I2d,I3d,I4d,I5d,~,Iv] = ElementRadiation(ring,lindata);
         [I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,lindata);
         I1=I1d+I1w;
         I2=I2d+I2w;
@@ -74,7 +69,7 @@
         I5=I5d+I5w;
         if ~isempty(Ux)
             U0 = Ux*1e6; %convert MeV to eV
-            fprintf('dipole radiation loss:  %4.5f keV\n', U0/1000.);
+            fprintf('Radiation loss:  %4.5f keV\n', U0/1000.);
         elseif is6d
             U0 = atgetU0(ring,'method','tracking');
         else
@@ -83,20 +78,9 @@
         Jx = 1-I4/I2;
         Jy = 1.00;
         Je = 2+I4/I2;
-        if offmom
-            warning('AT:NoEmit','\n%s\n%s\n%s',...
-                'For off-momentum lattices, the equilibrium emittance',...
-                'cannot be computed because the contribution of quadrupoles is missing.',...
-                'To avoid this warning, use ">> warning(''off'',''AT:NoEmit'')"');
-            emittx = NaN;
-            sigma_E = NaN;
-            alphac = NaN;
-        else
-            emittx = Cq*gamma^2*I5/(I2-I4);
-            sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4));
-            alphac = I1/circ;
-        end
-
+        emittx = Cq*gamma^2*I5/(I2-I4);
+        sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4));
+        alphac = I1/circ;
         % minimum emittance due to radiation 1/gamma cone (Handbook, Chao, eq23, pag 211)
         emitty_lim = 13/55*Cq/Jy/I2*Iv;
 

atmat/atphysics/Radiation/DipoleRadiation.m

@@ -7,105 +7,5 @@
 % R.H. Helm, M.J. Lee, P.L. Morton and M. Sands
 % SLAC-PUB-1193, March 1973
 
-
-angle=atgetfieldvalues(ring,'BendingAngle');
-isdipole=isfinite(angle) & (angle~=0);
-vini=lindata([isdipole;false])';
-vend=lindata([false;isdipole])';
-[di1,di2,di3,di4,di5,di6,div]=cellfun(@diprad,ring(isdipole),num2cell(vini),num2cell(vend));
-I1=sum(di1);
-I2=sum(di2);
-I3=sum(di3);
-I4=sum(di4);
-I5=sum(di5);
-I6=sum(di6);
-Iv=sum(div);
-
-    function [di1,di2,di3,di4,di5,di6,div]=diprad(elem,dini,dend)
-        betax0 = dini.beta(1);
-        betaz0 = dini.beta(2);
-        alphax0 = dini.alpha(1);
-        eta0 = dini.Dispersion(1);
-        etap0 = dini.Dispersion(2);
-
-        ll = elem.Length;
-        theta = elem.BendingAngle;
-        rho = ll / theta;
-        rho2 = rho * rho;
-        K = getfoc(elem);
-        kx2 = K + 1.0/rho2;
-        kz2 = -K;
-        eps1 = tan(elem.EntranceAngle) / rho;
-        eps2 = tan(elem.ExitAngle) / rho;
-
-        eta3 = dend.Dispersion(1);
-        alphax1 = alphax0 - betax0*eps1;
-        gammax1 = (1.0 + alphax1 * alphax1) / betax0;
-        alphaz1 = dini.alpha(2) + betaz0*eps1;
-        alphaz2 = dend.alpha(2) - dend.beta(2)*eps2;
-        gammaz1 = (1.0 + alphaz1*alphaz1) / betaz0;
-        etap1 = etap0 + eta0*eps1;
-        etap2 = dend.Dispersion(2) - eta3*eps2;
-
-        h0 = gammax1*eta0*eta0 + 2.0*alphax1*eta0*etap1 + betax0*etap1*etap1;
-
-        if kx2 ~= 0.0
-            if kx2 > 0.0         % Focusing
-                kl = ll * sqrt(kx2);
-                ss = sin(kl) / kl;
-                cc = cos(kl);
-            else                % Defocusing
-                kl = ll * sqrt(-kx2);
-                ss = sinh(kl) / kl;
-                cc = cosh(kl);
-            end
-            eta_ave = (theta - (etap2 - etap1)) / kx2 / ll;
-            bb = 2.0 * (alphax1*eta0 + betax0*etap1) * rho;
-            aa = -2.0 * (alphax1*etap1 + gammax1*eta0) * rho;
-            h_ave = h0 + (aa * (1.0-ss) + bb * (1.0-cc) / ll ...
-                          + gammax1 * (3.0-4.0*ss+ss*cc) / 2.0 / kx2 ...
-                          - alphax1 * (1.0-cc)^2 / kx2 / ll ...
-                          + betax0 * (1.0-ss*cc) / 2.0 ...
-                          ) / kx2 / rho2;
-        else
-            eta_ave = 0.5 * (eta0 + eta3) - ll*ll / 12.0 / rho;
-            hp0 = 2.0 * (alphax1 * eta0 + betax0 * etap1) / rho;
-            h2p0 = 2.0 * (-alphax1*etap1 + betax0/rho - gammax1*eta0) / rho;
-            h_ave = h0 + hp0*ll/2.0 + h2p0*ll*ll/6.0 ...
-                - alphax1*ll^3/4.0/rho2 ...
-                + gammax1*ll^4/20.0/rho2;
-        end
-        if kz2 ~= 0.0
-            bz_ave=(gammaz1 + kz2*betaz0 + (alphaz2-alphaz1)/ll)/2/kz2;
-        else
-            bz_ave=betaz0-alphaz1*ll + gammaz1*ll*ll/3;
-        end
-
-        di1 = eta_ave * ll / rho;
-        di2 = ll / rho2;
-        di3 = ll / abs(rho) / rho2;
-        di4 = eta_ave * ll * (2.0*K+1.0/rho2) / rho ...
-            - (eta0*eps1 + eta3*eps2)/rho;
-        di5 = h_ave * ll / abs(rho) / rho2;
-        di6 = kz2^2*eta_ave^2*ll;
-        div = abs(theta)/rho2*bz_ave;
-    end
-
-    function K=getfoc(elem)
-        if isfield(elem,'PolynomB') && length(elem.PolynomB) >= 2
-            K = elem.PolynomB(2);
-            if isfield(elem,'K') && elem.K ~= K
-                warning('AT:InconsistentK',...
-                'Values in K and PolynomB(2) are different. Using PolynomB(2)');
-            end
-        elseif isfield(elem,'K')
-            K = elem.K;
-            if K ~= 0
-                warning('AT:InconsistentK','PolynomB(2) is missing. Using K');
-            end
-        else
-            K = 0;
-        end
-    end
-end
+[I1,I2,I3,I4,I5,I6,Iv] = ElementRadiation(ring,lindata,'UseQuadrupole', false);
 

atmat/atphysics/Radiation/ElementRadiation.m

@@ -0,0 +1,147 @@
+function [I1,I2,I3,I4,I5,I6,Iv] = ElementRadiation(ring,lindata,varargin)
+%ELEMENTRADIATION	Compute the radiation integrals in dipoles
+%and quadrupoles
+
+% Analytical integration from:
+%
+% EVALUATION OF SYNCHROTRON RADIATION INTEGRALS
+% R.H. Helm, M.J. Lee, P.L. Morton and M. Sands
+% SLAC-PUB-1193, March 1973
+
+[quadon, ~]=getoption(varargin,'UseQuadrupoles', true);
+angle=atgetfieldvalues(ring,'BendingAngle');
+isdipole=isfinite(angle) & (angle~=0);
+
+if quadon
+    class=atgetcells(ring,'Class','Quadrupole');
+    pb = atgetfieldvalues(ring,'PolynomB',{2});
+    isquadrupole=class & (pb~=0);
+    iselement = isdipole | isquadrupole;
+else
+    iselement = isdipole;
+end
+vini=lindata([iselement;false])';
+vend=lindata([false;iselement])';
+
+[di1,di2,di3,di4,di5,di6,div]=cellfun(@elrad,ring(iselement),num2cell(vini),num2cell(vend));
+I1=sum(di1);
+I2=sum(di2);
+I3=sum(di3);
+I4=sum(di4);
+I5=sum(di5);
+I6=sum(di6);
+Iv=sum(div);
+
+    function [di1,di2,di3,di4,di5,di6,div]=elrad(elem,dini,dend)
+        betax0 = dini.beta(1);
+        betaz0 = dini.beta(2);
+        alphax0 = dini.alpha(1);
+        eta0 = dini.Dispersion(1);
+        etap0 = dini.Dispersion(2);
+        try
+            theta = elem.BendingAngle;
+        catch
+            xpi = dini.ClosedOrbit(2)/(1+dini.ClosedOrbit(5));
+            xpo = dend.ClosedOrbit(2)/(1+dini.ClosedOrbit(5));
+            theta = xpi-xpo;
+        end
+        if theta==0.0
+            di1 = 0.0;
+            di2 = 0.0;
+            di3 = 0.0;
+            di4 = 0.0;
+            di5 = 0.0;
+            di6 = 0.0;
+            div = 0.0;
+            return
+        end
+        try
+            ange = elem.EntranceAngle;
+        catch
+            ange = 0.0;
+        end
+        try
+            ango = elem.ExitAngle;
+        catch
+            ango = 0.0;
+        end
+        ll = elem.Length;       
+        rho = ll / theta;
+        rho2 = rho * rho;
+        K = getfoc(elem);
+        kx2 = K + 1.0/rho2;
+        kz2 = -K;
+        eps1 = tan(ange) / rho;
+        eps2 = tan(ango) / rho;
+
+        eta3 = dend.Dispersion(1);
+        alphax1 = alphax0 - betax0*eps1;
+        gammax1 = (1.0 + alphax1 * alphax1) / betax0;
+        alphaz1 = dini.alpha(2) + betaz0*eps1;
+        alphaz2 = dend.alpha(2) - dend.beta(2)*eps2;
+        gammaz1 = (1.0 + alphaz1*alphaz1) / betaz0;
+        etap1 = etap0 + eta0*eps1;
+        etap2 = dend.Dispersion(2) - eta3*eps2;
+
+        h0 = gammax1*eta0*eta0 + 2.0*alphax1*eta0*etap1 + betax0*etap1*etap1;
+
+        if kx2 ~= 0.0
+            if kx2 > 0.0         % Focusing
+                kl = ll * sqrt(kx2);
+                ss = sin(kl) / kl;
+                cc = cos(kl);
+            else                % Defocusing
+                kl = ll * sqrt(-kx2);
+                ss = sinh(kl) / kl;
+                cc = cosh(kl);
+            end
+            eta_ave = (theta - (etap2 - etap1)) / kx2 / ll;
+            bb = 2.0 * (alphax1*eta0 + betax0*etap1) * rho;
+            aa = -2.0 * (alphax1*etap1 + gammax1*eta0) * rho;
+            h_ave = h0 + (aa * (1.0-ss) + bb * (1.0-cc) / ll ...
+                          + gammax1 * (3.0-4.0*ss+ss*cc) / 2.0 / kx2 ...
+                          - alphax1 * (1.0-cc)^2 / kx2 / ll ...
+                          + betax0 * (1.0-ss*cc) / 2.0 ...
+                          ) / kx2 / rho2;
+        else
+            eta_ave = 0.5 * (eta0 + eta3) - ll*ll / 12.0 / rho;
+            hp0 = 2.0 * (alphax1 * eta0 + betax0 * etap1) / rho;
+            h2p0 = 2.0 * (-alphax1*etap1 + betax0/rho - gammax1*eta0) / rho;
+            h_ave = h0 + hp0*ll/2.0 + h2p0*ll*ll/6.0 ...
+                - alphax1*ll^3/4.0/rho2 ...
+                + gammax1*ll^4/20.0/rho2;
+        end
+        if kz2 ~= 0.0
+            bz_ave=(gammaz1 + kz2*betaz0 + (alphaz2-alphaz1)/ll)/2/kz2;
+        else
+            bz_ave=betaz0-alphaz1*ll + gammaz1*ll*ll/3;
+        end
+
+        di1 = eta_ave * ll / rho;
+        di2 = ll / rho2;
+        di3 = ll / abs(rho) / rho2;
+        di4 = eta_ave * ll * (2.0*K+1.0/rho2) / rho ...
+            - (eta0*eps1 + eta3*eps2)/rho;
+        di5 = h_ave * ll / abs(rho) / rho2;
+        di6 = kz2^2*eta_ave^2*ll;
+        div = abs(theta)/rho2*bz_ave;
+    end
+
+    function K=getfoc(elem)
+        if isfield(elem,'PolynomB') && length(elem.PolynomB) >= 2
+            K = elem.PolynomB(2);
+            if isfield(elem,'K') && elem.K ~= K
+                warning('AT:InconsistentK',...
+                'Values in K and PolynomB(2) are different. Using PolynomB(2)');
+            end
+        elseif isfield(elem,'K')
+            K = elem.K;
+            if K ~= 0
+                warning('AT:InconsistentK','PolynomB(2) is missing. Using K');
+            end
+        else
+            K = 0;
+        end
+    end
+end
+

pyat/at/lattice/lattice_object.py

@@ -15,7 +15,7 @@
 import copy
 import numpy
 import math
-from typing import Optional, Union, Tuple
+from typing import Optional, Union
 if sys.version_info.minor < 9:
     from typing import Callable, Iterable, Generator
     SupportsIndex = int
@@ -529,6 +529,31 @@ def copy(self) -> "Lattice":
     def deepcopy(self) -> "Lattice":
         """Returns a deep copy of the lattice"""
         return copy.deepcopy(self)
+
+    def slice_elements(self, refpts, slices: Optional[int] = 1) \
+            -> "Lattice":
+        """Create a new lattice by slicing the elements at refpts
+
+        Keyword arguments:
+            size=None:      Length of a slice. Default: computed from the
+              range and number of points: ``size = (s_max-s_min)/slices``.
+            slices=1:       Number of slices in the specified range. Ignored if
+              size is specified. Default: no slicing
+
+        Returns:
+            newring:    New Lattice object
+        """
+        def slice_generator(_):
+            check = get_bool_index(self, refpts)
+            for el, ok in zip(self, check):
+                if ok and (slices>1):
+                    frac = numpy.ones(slices) / slices
+                    for elem in el.divide(frac):
+                        yield elem
+                else:
+                    yield el
+
+        return Lattice(slice_generator, iterator=self.attrs_filter)        
 
     def slice(self, size: Optional[float] = None, slices: Optional[int] = 1) \
             -> "Lattice":