Simplified radiation integrals in wigglers

atmat/atphysics/LinearOptics/amat.m

@@ -49,19 +49,14 @@
 % n2x n2y n2z
 % n3x n3y n3z
 nn=0.5*abs(sqrt(-1i*Vn'*jm{dms}*Vxyz{dms}));
-[~,ind]=max(nn(select,select));
 
-indS=ind;
-%sort ind for the case of repeated indices
-%[~,indS]=sort(ind,2);
-%logic still not quite correct- should find best projector on x then y then
-%z.  Can the one vector project best on both x and y?  Is this projection
-%definition correct?
-%14 Oct, 2014- more testing needed before using this additional sort.
-%B. Nash
-
-%reorder pairs
-V_ordered=Vn(:,2*indS-1);
+rows=select;
+V_ordered=[];
+for ixz=select
+    [~,ind]=max(nn(rows,ixz));
+    V_ordered=[V_ordered Vn(:,rows(ind))]; %#ok<AGROW>
+    rows(ind)=[];
+end
 
 %build a matrix
 a=reshape([real(V_ordered);imag(V_ordered)],nv,nv);

atmat/atphysics/ParameterSummaryFunctions/atsummary.m

@@ -31,7 +31,7 @@
 params   = atgetcells(ring,'Class','RingParam');
 if any(params)
     LatticeName=ring{find(params,1)}.FamName;
-elseif isfield(E0,'LatticeFile')
+elseif isfield(GLOBVAL,'LatticeFile')
     LatticeName=GLOBVAL.LatticeFile;
 else
     LatticeName='';
@@ -53,14 +53,16 @@
 smm.revTime       = smm.circumference / cspeed;
 smm.revFreq       = cspeed / smm.circumference;
 smm.gamma         = 1.e9 * smm.e0 / e_mass;
-smm.beta          = sqrt(1 - 1/smm.gamma);
+smm.beta          = sqrt(1 - 1/smm.gamma/smm.gamma);
+beta_c            = smm.beta * cspeed;
 
 [TD, smm.tunes, smm.chromaticity] = twissring(ring, 0, 1:length(ring)+1, 'chrom', 1e-8);
 smm.compactionFactor = mcf(ring);
 
 % For calculating the synchrotron integrals
-[I1,I2,I3,I4,I5,I6,~] = DipoleRadiation(ring,TD);
-smm.integrals=[I1,I2,I3,I4,I5,I6];
+[I1d,I2d,I3d,I4d,I5d,I6,~] = DipoleRadiation(ring,TD);
+[I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,TD);
+smm.integrals=[I1d+I1w,I2d+I2w,I3d+I3w,I4d+I4w,I5d+I5w,I6];
 
 % Damping numbers
 % Use Robinson's Theorem
@@ -87,11 +89,15 @@
 % references:  H. Winick, "Synchrotron Radiation Sources: A Primer",
 % World Scientific Publishing, Singapore, pp92-95. (1995)
 % Wiedemann, pp290,350. Chao, pp189.
-smm.syncphase = pi - asin(1/smm.overvoltage);
+if smm.overvoltage > 1
+    smm.syncphase = pi - asin(1/smm.overvoltage);
+else
+    smm.syncphase = NaN;
+end
 smm.energyacceptance = sqrt(v_cav*sin(smm.syncphase)*2*(sqrt(smm.overvoltage^2-1) - acos(1/smm.overvoltage))/(pi*smm.harmon*abs(smm.etac)*smm.e0*1e9));
-smm.synctune = sqrt((smm.etac*smm.harmon*v_cav*cos(smm.syncphase))/(2*pi*smm.e0*1e9));
-smm.bunchlength = smm.beta*cspeed*abs(smm.etac)*smm.naturalEnergySpread/(smm.synctune*smm.revFreq*2*pi);
-
+smm.synctune = sqrt((smm.etac*smm.harmon*v_cav*cos(smm.syncphase))/(2*pi*smm.e0*1e9))/smm.beta;
+bunchtime = abs(smm.etac)*smm.naturalEnergySpread/(smm.synctune*smm.revFreq*2*pi);
+smm.bunchlength = beta_c * bunchtime;
 % optics
 % [bx by] = modelbeta;
 % [ax ay] = modelalpha;
@@ -146,7 +152,7 @@
     fprintf('   Synchronous Phase:  \t\t% 4.5f [rad] (%4.5f [deg])\n', smm.syncphase, smm.syncphase*180/pi);
     fprintf('   Linear Energy Acceptance:  \t% 4.3f %%\n', smm.energyacceptance*100);
     fprintf('   Synchrotron Tune:   \t\t% 4.5f (%4.5f kHz or %4.2f turns) \n', smm.synctune, smm.synctune/smm.revTime*1e-3, 1/smm.synctune);
-    fprintf('   Bunch Length:       \t\t% 4.5f [mm] (%4.5f ps)\n', smm.bunchlength*1e3, smm.bunchlength/3e8*1e12);
+    fprintf('   Bunch Length:       \t\t% 4.5f [mm] (%4.5f ps)\n', smm.bunchlength*1e3, bunchtime*1e12);
     fprintf(SeparatorString);
     fprintf('   Injection:\n');
     fprintf('   H: beta = %06.3f [m] alpha = %+04.1e eta = %+04.3f [m] eta'' = %+04.1e \n', bx(1), ax(1), etax(1), etaprimex(1));

atmat/atphysics/ParameterSummaryFunctions/atx.m

@@ -81,7 +81,7 @@
     function [linusr,pm]=atx2(ring,energy,periods,voltage,eloss,varargin)
 
         c = PhysConstant.speed_of_light_in_vacuum.value;
-        [dpp,refusr]=parseargs({0,1:length(ring)},varargin);
+        [dpp,refusr]=getargs(varargin,0,1:length(ring));
         if islogical(refusr)
             refusr(end+1,length(ring)+1)=false;
         else
@@ -252,8 +252,12 @@
 
             function [chi,nu]=decode(range)
                 matr=Rmat(range,range);
-                nu=atan2(matr(1,2)-matr(2,1),matr(1,1)+matr(2,2))/2/pi;
-                chi=-log(sqrt(det(matr)));
+%                 nu=mod(atan2(matr(1,2)-matr(2,1),matr(1,1)+matr(2,2))/2/pi,1);
+%                 chi=-log(sqrt(det(matr)));
+                ev=eig(matr);
+                ev1log=log(ev(1));
+                chi=-real(ev1log);
+                nu=mod(imag(ev1log)/2/pi,1);
             end
         end
 

atmat/atphysics/ParameterSummaryFunctions/ringpara.m

@@ -30,6 +30,9 @@
 %radiation effects removed until repaired...
 %Analytical computation of radiation integrals
 
+%Modified by Laurent Farvacque, 2020-09-24
+%wiggler radiation effect reintroduced
+
 global THERING
 
 e_mass=PhysConstant.electron_mass_energy_equivalent_in_MeV.value*1e6;	% eV
@@ -58,15 +61,16 @@
 end
 
 gamma = energy/e_mass;
+betar = sqrt(1-1/gamma/gamma);
+beta_c = betar*cspeed;
 [lindata,tune,chrom]=atlinopt(ring,0.0,1:length(ring)+1);
-[I1,I2,I3,I4,I5,I6,Iv] = DipoleRadiation(ring,lindata);
-
-% [I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,lindata);
-% I1=I1+I1w;
-% I2=I2+I2w;
-% I3=I3+I3w;
-% I4=I4+I4w;
-% I5=I5+I5w;
+[I1d,I2d,I3d,I4d,I5d,~,Iv] = DipoleRadiation(ring,lindata);
+[I1w,I2w,I3w,I4w,I5w] = WigglerRadiation(ring,lindata);
+I1=I1d+I1w;
+I2=I2d+I2w;
+I3=I3d+I3w;
+I4=I4d+I4w;
+I5=I5d+I5w;
 
 spos=findspos(ring,1:length(ring)+1);
 circ=spos(end);
@@ -100,16 +104,8 @@
 rp.dampingtime = 1./[alphax, alphay, alphaE];
 
 
-%compute coupled damping times
-%[nu,chi]=atTunesAndDampingRatesFromMat(findm66(atradon(ring)));
-try
-    [~,chi]=atTunesAndDampingRatesFromMat(findm66((ring)));
-catch
-    warning('failed coupled damping times computation');
-    chi=[NaN,NaN,NaN];
-end
-
-rp.coupleddampingtime=T0./chi;
+%computation of coupled damping time removed since radiation must be off
+%coupled damping times are available in atx
 
 rp.dampingJ = [Jx,Jy,Je];
 
@@ -126,12 +122,17 @@
     freq_rf = 476.314e6;
 end
 
-phi_s = pi-asin(rp.U0/Vrf);
-nus = sqrt(harm*Vrf*abs(rp.etac*cos(phi_s))/2/pi/rp.E0);
+if rp.U0 <= Vrf
+    phi_s = pi-asin(rp.U0/Vrf);
+else
+    phi_s = NaN;
+end
+nus = sqrt(harm*Vrf*abs(rp.etac*cos(phi_s))/2/pi/rp.E0)/betar;
 rp.nus = nus;
 rp.phi_s = phi_s;
 rp.harm = harm;
-rp.bunchlength = rp.sigma_E*harm*abs(rp.etac)/nus/2/pi/freq_rf*cspeed; % rms bunchlength in meter
+bunchtime=rp.sigma_E*harm*abs(rp.etac)/nus/2/pi/freq_rf;    % [s]
+rp.bunchlength = beta_c*bunchtime;                          % [m]
 delta_max = sqrt(2*U0/pi/alphac/harm/rp.E0)*sqrt( sqrt((Vrf/U0).^2-1) - acos(U0./Vrf));
 rp.delta_max = delta_max;
 
@@ -160,7 +161,7 @@
 
 if nargout == 0
     fprintf('\n');
-    fprintf('   ******** AT Ring Parmater Summary ********\n');
+    fprintf('   ******** AT Ring Parameter Summary ********\n');
     fprintf('   Energy: \t\t\t%4.5f [GeV]\n', rp.E0/1E9);
     fprintf('   Circumference: \t\t%4.5f [m]\n', rp.R*2*pi);
     fprintf('   Revolution time: \t\t%4.5f [ns] (%4.5f [MHz]) \n', rp.T0*1e9,1./rp.T0*1e-6);
@@ -191,7 +192,7 @@
     fprintf('   Synchronous Phase:  %4.5f [rad] (%4.5f [deg])\n', rp.phi_s, rp.phi_s*180/pi);
     fprintf('   Linear Energy Acceptance:  %4.5f %%\n', rp.delta_max*100);
     fprintf('   Synchrotron Tune:   %4.5f (%4.5f kHz or %4.2f turns) \n', rp.nus, rp.nus/rp.T0*1e-3, 1/rp.nus);
-    fprintf('   Bunch Length:       %4.5f [mm], %4.5f [ps]\n', rp.bunchlength*1e3, rp.bunchlength/cspeed*1e12);
+    fprintf('   Bunch Length:       %4.5f [mm], %4.5f [ps]\n', rp.bunchlength*1e3, bunchtime*1e12);
     fprintf('\n');
 %     fprintf('   Vertical Emittance:  %4.5f [nm]\n', rp.emitty*1e9);
 %     fprintf('   Emitty from Dy:  %4.5f [nm], from linear coupling: %4.5f\n', rp.emitty_d*1e9,rp.emitty_c*1e9);

atmat/atphysics/Radiation/WigglerRadiation.m

@@ -0,0 +1,114 @@
+function [I1,I2,I3,I4,I5] = WigglerRadiation(ring,lindata)
+%WIGGLERRADIATION	Compute the radiation integrals in wigglers
+%
+%[I1,I2,I3,I4,I5] = WIGGLERRADIATION(RING,LINDATA)
+%
+%RING       Lattice structure
+%LINDATA    Output of atlinopt for all lattice elements
+%
+%WigglerRadiation computes the radiation integrals for all wigglers with
+%the following approximations:
+%
+%- The self-induced dispersion is neglected in I4 and I5, but is is used as
+%  a lower limit for the I5 contribution
+%
+% I1, I2 are integrated analytically
+% I3 is integrated analytically for a single harmonic, numerically otherwise
+
+nstep=60;
+e_mass=PhysConstant.electron_mass_energy_equivalent_in_MeV.value*1e6;	% eV
+cspeed = PhysConstant.speed_of_light_in_vacuum.value;                   % m/s
+
+iswiggler=@(elem) strcmp(elem.Class,'Wiggler') && ~strcmp(elem.PassMethod,'DriftPass');
+wigglers=cellfun(iswiggler, ring);
+if any(wigglers)
+    energy=unique(atgetfieldvalues(ring(wigglers),'Energy'));
+    if length(energy) > 1
+        error('AT:NoEnergy','Energy field not equal for all elements');
+    end
+    Brho=sqrt(energy*energy - e_mass*e_mass)/cspeed;
+    vini=lindata([wigglers;false])';
+    [di1,di2,di3,di4,di5]=cellfun(@wigrad,ring(wigglers),num2cell(vini));
+    I1=sum(di1);
+    I2=sum(di2);
+    I3=sum(di3);
+    I4=sum(di4);
+    I5=sum(di5);
+else
+    [I1,I2,I3,I4,I5]=deal(0);
+end
+
+    function [di1,di2,di3,di4,di5]=wigrad(elem,dini)
+        le=elem.Length;
+        alphax0=dini.alpha(1);
+        betax0=dini.beta(1);
+        gammax0=(alphax0.*alphax0+1)./betax0;
+        eta0 = dini.Dispersion(1);
+        etap0 = dini.Dispersion(2);
+        H0=gammax0*eta0*eta0 + 2*alphax0*eta0*etap0 + betax0*etap0*etap0;
+        avebetax=betax0+alphax0*le+gammax0*le*le/3;
+
+        kw=2*pi/elem.Lw;
+        rhoinv=elem.Bmax/Brho;
+        coefh=elem.By(2,:)*rhoinv;
+        coefv=elem.Bx(2,:)*rhoinv;
+        coef2=[coefh coefv];
+        if length(coef2)==1     % Analytical I3
+            di3=coef2^3*4*le/3/pi;
+        else                    % Numerical I3
+            [bx,bz]=Baxis(elem,linspace(0,elem.Lw,nstep+1));
+            B2=bx.*bx+bz.*bz;
+            rinv=sqrt(B2)/Brho;
+            di3=trapz(rinv.^3)*le/nstep;
+        end
+        di2=elem.Length*(coefh*coefh'+coefv*coefv')/2;
+        di1=-di2/kw/kw;
+        di4=0;
+        if ~isempty(coefh)
+            d5lim=4*avebetax*le*coefh(1)^5/15/pi/kw/kw;
+        else
+            d5lim = 0;
+        end
+        di5=max(H0*di3,d5lim);
+%         fprintf('%s\t%e\t%e\t%e\t%e\t%e\t(%e,%e,%e)\n',elem.FamName,di1,di2,di3,di4,di5,d5lim,H0*di3,dini.Dispersion(1));
+%         avebetaz=dini.beta(2)+dini.alpha(2)*le+(1+dini.alpha(2)^2)/dini.beta(2)*le*le/3;
+%         dnuz=avebetaz*di2/8/pi;
+%         fprintf('%s\t%e\t%e\n',elem.FamName,dnuz,avebetaz)
+    end
+
+    function [bx,bz] = Baxis(wig,s)
+        %Bwig   Compute the field on the axis of a generic wiggler
+        %
+        %The field of an horizontal wiggler is represented by a sum of harmonics,
+        %each being described by a 6x1 column in a matrix By such as:
+        %
+        %   Bz/Bmax = -By2 * cos(By5*kw*s + By6)
+        %
+        %The field of an vertical wiggler is represented by a sum of harmonics,
+        %each being described by a 6x1 column in a matrix Bx such as:
+        %
+        %   Bx/Bmax =  Bx2 * cos(Bx5*kw*s + Bx6)
+
+        kw=2*pi/wig.Lw;
+        Bmax=wig.Bmax;
+        kws=kw*s;
+        [bxh,bzh]=cellfun(@harmh,num2cell(wig.By,1),'UniformOutput',false);
+        [bxv,bzv]=cellfun(@harmv,num2cell(wig.Bx,1),'UniformOutput',false);
+        bx=sum(cat(3,bxh{:},bxv{:}),3);
+        bz=sum(cat(3,bzh{:},bzv{:}),3);
+
+        function [bx,bz]=harm(pb)
+            bz=-Bmax*pb(2)*cos(pb(5)*kws+pb(6));
+            bx=zeros(size(kws));
+        end
+        function [bx,bz]=harmh(pb)
+            [bx,bz]=harm(pb);
+        end
+        function [bx,bz]=harmv(pb)
+            [bz,fz]=harm(pb);
+            bx=-fz;
+        end
+    end
+
+end
+

atmat/atphysics/Radiation/atenergy.m

@@ -88,13 +88,32 @@
 end
 
 % get energy loss per turn
-if nargout >= 5 
+if nargout >= 5
     % Losses = Cgamma/2/pi*EGeV^4*I2
-    cgamma=4e9*pi*PhysConstant.classical_electron_radius.value/3/...
-        PhysConstant.electron_mass_energy_equivalent_in_MeV.value^3; % [m/GeV^3]
+    e_mass=PhysConstant.electron_mass_energy_equivalent_in_MeV.value*1e6;	% eV
+    cspeed = PhysConstant.speed_of_light_in_vacuum.value;                   % m/s
+    e_radius=PhysConstant.classical_electron_radius.value;                  % m
+    Cgamma=4.0E27*pi*e_radius/3/e_mass^3;                       % [m/GeV^3]
+    Brho=sqrt(energy*energy - e_mass*e_mass)/cspeed;
+
+    % Dipole radiation
     lendp=atgetfieldvalues(ring(dipoles),'Length');
-    losses=atgetfieldvalues(ring(atgetcells(ring,'I2')),'I2');
-    I2=nbper*(sum(abs(theta.*theta./lendp))+sum(losses));            % [m-1]
-    U0=cgamma/2/pi*(energy*1.e-9)^4*I2*1e9;                          % [eV]
+    I2d=sum(abs(theta.*theta./lendp));
+    % Wiggler radiation
+    iswiggler=@(elem) strcmp(elem.Class,'Wiggler') && ~strcmp(elem.PassMethod,'DriftPass');
+    wigglers=cellfun(iswiggler, ring);
+    I2w=sum(cellfun(@wiggler_i2,ring(wigglers)));
+    % Additional radiation
+    I2x=sum(atgetfieldvalues(ring(atgetcells(ring,'I2')),'I2'));
+    I2=nbper*(I2d+I2w+I2x);                                     % [m-1]
+    U0=Cgamma/2/pi*(energy*1.e-9)^4*I2*1e9;                     % [eV]
 end
+
+    function i2=wiggler_i2(elem)
+        rhoinv=elem.Bmax/Brho;
+        coefh=elem.By(2,:);
+        coefv=elem.Bx(2,:);
+        i2=elem.Length*(coefh*coefh'+coefv*coefv')*rhoinv*rhoinv/2;
+    end
+
 end

atmat/atphysics/Radiation/atradoff.m

@@ -4,6 +4,8 @@
 %   [RING2,RADINDEX,CAVINDEX] = ATRADOFF(RING,CAVIPASS,BENDPASS,QUADPASS)
 %    Changes passmethods to turn off radiation damping. 
 %
+%   The default is to turn cavities OFF and remove radiation in all magnets.
+%
 %  INPUTS
 %  1. RING      initial AT structure
 %  2. CAVIPASS  pass method for cavities
@@ -26,6 +28,7 @@
 %   'quadpass'      pass method for quadrupoles. Default 'auto'
 %   'sextupass'     pass method for sextupoles. Default 'auto'
 %   'octupass'      pass method for bending magnets. Default 'auto'
+%   'wigglerpass'	pass method for wigglers. Default 'auto'
 %
 %   OUPUTS
 %   1. RING2     Output ring
@@ -37,7 +40,7 @@
 [octupass,varargs]=getoption(varargin,'octupass','auto');
 [sextupass,varargs]=getoption(varargs,'sextupass','auto');
 [quadpass,varargs]=getoption(varargs,'quadpass','auto');
-%[wigglerpass,varargs]=getoption(varargs,'wigglerpass','auto');
+[wigglerpass,varargs]=getoption(varargs,'wigglerpass','auto');
 [bendpass,varargs]=getoption(varargs,'bendpass','auto');
 [cavipass,varargs]=getoption(varargs,'cavipass','auto');
 [cavipass,bendpass,quadpass]=getargs(varargs,cavipass,bendpass,quadpass);
@@ -52,7 +55,7 @@
 
 [ring,octupoles]=changepass(ring,octupass,@(rg) selmpole(rg,4),@autoMultiPolePass,'Octu');
 
-%[ring,wigglers]=changepass(ring,wigglerpass,@selwiggler,@autoMultiPolePass,'Wiggler');
+[ring,wigglers]=changepass(ring,wigglerpass,@(rg) atgetcells(rg,'Class','Wiggler'),@autoMultiPolePass,'Wiggler');
 
 cavitiesIndex=atgetcells(ring,'PassMethod',@(elem,pass) endsWith(pass,'CavityPass'));
 radelemIndex=atgetcells(ring,'PassMethod',@(elem,pass) endsWith(pass,'RadPass'));
@@ -61,7 +64,7 @@
     atdisplay(1,['Cavities modified at position ' num2str(find(cavities)')]);
 end
 
-radnum=sum(dipoles|quadrupoles|sextupoles|octupoles);
+radnum=sum(dipoles|quadrupoles|sextupoles|octupoles|wigglers);
 if radnum > 0
     atdisplay(1,[num2str(radnum) ' elements switched to include radiation']);
 end