Skip to content

Commit

Permalink
Simplified radiation integrals in wigglers (#188)
Browse files Browse the repository at this point in the history
* Wiggler field, wiggler integrals

* Correct global variable name (#187)

* corrected Bwig

* Bug fixes

* Bug fixes

* Bug fixes

* Updated atsummary and removed debug info

* Analytical I3

* amat, atradon

* Fixed missing 'beta' in sync tune and bunch length

* Fixed missing 'beta' in sync tune and bunch length

* Removed temp files

* Fix import of wigglers from matlab

* Check overvoltage < 1

* Wiggler radiation in python

* Remove useless import, bug fix

* Improved docstring

* vector ordering in amat

* fix coupled tunes

* Handle wigglers in atradon and atradoff

* Ignore wigglers using DriftPass

Co-authored-by: Laurent Farvacque <laurent.farvacque@gmail.com>
  • Loading branch information
lfarv and lfarv committed Oct 12, 2020
1 parent 52cac25 commit e4c7a77
Show file tree
Hide file tree
Showing 13 changed files with 381 additions and 148 deletions.
19 changes: 7 additions & 12 deletions atmat/atphysics/LinearOptics/amat.m
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
22 changes: 14 additions & 8 deletions atmat/atphysics/ParameterSummaryFunctions/atsummary.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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));
Expand Down
10 changes: 7 additions & 3 deletions atmat/atphysics/ParameterSummaryFunctions/atx.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
47 changes: 24 additions & 23 deletions atmat/atphysics/ParameterSummaryFunctions/ringpara.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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];

Expand All @@ -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;

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
114 changes: 114 additions & 0 deletions atmat/atphysics/Radiation/WigglerRadiation.m
Original file line number Diff line number Diff line change
@@ -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

31 changes: 25 additions & 6 deletions atmat/atphysics/Radiation/atenergy.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
9 changes: 6 additions & 3 deletions atmat/atphysics/Radiation/atradoff.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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'));
Expand All @@ -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
Expand Down
Loading

0 comments on commit e4c7a77

Please sign in to comment.