feat: compute longitudinal RDT fluctuations. (#685)

* feat: compute longitudinal RDT fluctuations.
* fix: add input arguments validation.

atmat/atphysics/NonLinearDynamics/RDTbuildupFluct.m

@@ -0,0 +1,182 @@
+function [h21000s, h30000s, h10110s, h10020s, h10200s, h20001s, h00201s, h10002s,...
+          h22000s, h11110s, h00220s,...
+          h31000s, h40000s, h20110s, h11200s, h20020s, h20200s, h00310s, h00400s] = RDTbuildupFluct(betax,betay,...
+          etax,phix,phiy,b2L,b3L,b4L,nData)
+%BUILDUPRDTFLUCTUATION build-up fluctuation of RDTs
+%
+%   This function can show the build-up and cancellation of RDTs.
+%   DON'T use this function directly!!! 
+%   use computeRDTfluctuation()
+
+h21000s = complex(zeros(nData, 1));
+h30000s = complex(zeros(nData, 1));
+h10110s = complex(zeros(nData, 1));
+h10020s = complex(zeros(nData, 1));
+h10200s = complex(zeros(nData, 1));
+h20001s = complex(zeros(nData, 1));
+h00201s = complex(zeros(nData, 1));
+h10002s = complex(zeros(nData, 1));
+h22000s = complex(zeros(nData, 1));
+h11110s = complex(zeros(nData, 1));
+h00220s = complex(zeros(nData, 1));
+h31000s = complex(zeros(nData, 1));
+h40000s = complex(zeros(nData, 1));
+h20110s = complex(zeros(nData, 1));
+h11200s = complex(zeros(nData, 1));
+h20020s = complex(zeros(nData, 1));
+h20200s = complex(zeros(nData, 1));
+h00310s = complex(zeros(nData, 1));
+h00400s = complex(zeros(nData, 1));
+
+h21000 = 0;
+h30000 = 0;
+h10110 = 0;
+h10020 = 0;
+h10200 = 0;
+h20001 = 0;
+h00201 = 0;
+h10002 = 0;
+h22000 = 0;
+h11110 = 0;
+h00220 = 0;
+h31000 = 0;
+h40000 = 0;
+h20110 = 0;
+h11200 = 0;
+h20020 = 0;
+h20200 = 0;
+h00310 = 0;
+h00400 = 0;
+for ii=1:(nData-1)
+    h21000s(ii) = h21000;
+    h30000s(ii) = h30000;
+    h10110s(ii) = h10110;
+    h10020s(ii) = h10020;
+    h10200s(ii) = h10200;
+    h20001s(ii) = h20001;
+    h00201s(ii) = h00201;
+    h10002s(ii) = h10002;
+    h22000s(ii) = h22000;
+    h11110s(ii) = h11110;
+    h00220s(ii) = h00220;
+    h31000s(ii) = h31000;
+    h40000s(ii) = h40000;
+    h20110s(ii) = h20110;
+    h11200s(ii) = h11200;
+    h20020s(ii) = h20020;
+    h20200s(ii) = h20200;
+    h00310s(ii) = h00310;
+    h00400s(ii) = h00400;
+
+    betax_i = betax(ii);
+    betay_i = betay(ii);
+    etax_i = etax(ii);
+    phix_i = phix(ii);
+    phiy_i = phiy(ii);
+    if b2L(ii)~=0
+        b2l = b2L(ii);
+        h20001 = h20001 + b2l * betax_i * exp(0 + 2i * phix_i) / 8;
+        h00201 = h00201 - b2l * betay_i * exp(0 + 2i * phiy_i) / 8;
+        h10002 = h10002 + b2l * sqrt(betax_i) * etax_i * exp(0 + 1i * phix_i) / 2;
+    end
+    if b3L(ii)~=0
+        b3l = b3L(ii);
+        h20001 = h20001 - b3l * betax_i * etax_i * exp(0 + 1i * 2 * phix_i) / 4;
+        h00201 = h00201 + b3l * betay_i * etax_i * exp(0 + 1i * 2 * phiy_i) / 4;
+        h10002 = h10002 - b3l * sqrt(betax_i) * etax_i^2 * exp(0 + 1i * phix_i) / 2;
+
+        h21000j = -b3l * betax_i^1.5 * exp(1i * phix_i) / 8;
+        h30000j = -b3l * betax_i^1.5 * exp(1i * 3 * phix_i) / 24;
+        h10110j = b3l * sqrt(betax_i) * betay_i * exp(1i * phix_i) / 4;
+        h10020j = b3l * sqrt(betax_i) * betay_i * exp(1i * (phix_i - 2 * phiy_i)) / 8;
+        h10200j = b3l * sqrt(betax_i) * betay_i * exp(1i * (phix_i + 2 * phiy_i)) / 8;
+
+        h12000j = conj(h21000j);
+        h01110j = conj(h10110j);
+        h01200j = conj(h10020j);
+
+        h12000 = conj(h21000);
+        h01110 = conj(h10110);
+        h01200 = conj(h10020);
+
+        h22000 = h22000 + 1i * ((h21000 * h12000j - h12000 * h21000j) * 3 ...
+                  + (h30000 * conj(h30000j) - conj(h30000) * h30000j) * 9);
+
+        h11110 = h11110 + 1i * ((h21000 * h01110j - h01110 * h21000j) * 2 ...
+                  - (h12000 * h10110j - h10110 * h12000j) * 2 ...
+                  - (h10020 * h01200j - h01200 * h10020j) * 4 ...
+                  + (h10200 * conj(h10200j) - conj(h10200) * h10200j) * 4);
+
+        h00220 = h00220 + 1i * ((h10020 * h01200j - h01200 * h10020j) ...
+                  + (h10200 * conj(h10200j) - conj(h10200) * h10200j) ...
+                  + (h10110 * h01110j - h01110 * h10110j));
+
+        h31000 = h31000 + 1i * (h30000 * h12000j - h12000 * h30000j) * 6;
+
+        h40000 = h40000 + 1i * (h30000 * h21000j - h21000 * h30000j) * 3;
+
+        h20110 = h20110 + 1i * ((h30000 * h01110j - h01110 * h30000j) * 3 ...
+                  - (h21000 * h10110j - h10110 * h21000j) ...
+                  + (h10200 * h10020j - h10020 * h10200j) * 4);
+
+        h11200 = h11200 + 1i * ((h10200 * h12000j - h12000 * h10200j) * 2 ...
+                  + (h21000 * h01200j - h01200 * h21000j) * 2 ...
+                  + (h10200 * h01110j - h01110 * h10200j) * 2 ...
+                  - (h10110 * h01200j - h01200 * h10110j) * 2);
+
+        h20020 = h20020 + 1i * (-(h21000 * h10020j - h10020 * h21000j) ...
+                  + (h30000 * conj(h10200j) - conj(h10200) * h30000j) * 3 ...
+                  + (h10110 * h10020j - h10020 * h10110j) * 2);
+
+        h20200 = h20200 + 1i * ((h30000 * h01200j - h01200 * h30000j) * 3 ...
+                  + (h10200 * h21000j - h21000 * h10200j) ...
+                  - (h10110 * h10200j - h10200 * h10110j) * 2);
+
+        h00310 = h00310 + 1i * ((h10200 * h01110j - h01110 * h10200j) ...
+                  + (h10110 * h01200j - h01200 * h10110j));
+
+        h00400 = h00400 + 1i * (h10200 * h01200j - h01200 * h10200j);
+
+        h21000 = h21000 + h21000j;
+        h30000 = h30000 + h30000j;
+        h10110 = h10110 + h10110j;
+        h10020 = h10020 + h10020j;
+        h10200 = h10200 + h10200j;
+    end
+    if b4L(ii)~=0
+        b4l = b4L(ii);
+        h22000 = h22000 - 3 * b4l * betax_i^2 / 32;
+        h11110 = h11110 + 3 * b4l * betax_i * betay_i / 8;
+        h00220 = h00220 - 3 * b4l * betay_i^2 / 32;
+
+        h31000 = h31000 - b4l * betax_i^2 * exp(1i * 2 * phix_i) / 16;
+        h40000 = h40000 - b4l * betax_i^2 * exp(1i * 4 * phix_i) / 64;
+        h20110 = h20110 + 3 * b4l * betax_i * betay_i * exp(1i * 2 * phix_i) / 16;
+        h11200 = h11200 + 3 * b4l * betax_i * betay_i * exp(1i * 2 * phiy_i) / 16;
+        h20020 = h20020 + 3 * b4l * betax_i * betay_i * exp(1i * (2 * phix_i - 2 * phiy_i)) / 32;
+        h20200 = h20200 + 3 * b4l * betax_i * betay_i * exp(1i * (2 * phix_i + 2 * phiy_i)) / 32;
+        h00310 = h00310 - b4l * betay_i^2 * exp(1i * 2 * phiy_i) / 16;
+        h00400 = h00400 - b4l * betay_i^2 * exp(1i * 4 * phiy_i) / 64;
+    end
+end
+    h21000s(nData) = h21000;
+    h30000s(nData) = h30000;
+    h10110s(nData) = h10110;
+    h10020s(nData) = h10020;
+    h10200s(nData) = h10200;
+    h20001s(nData) = h20001;
+    h00201s(nData) = h00201;
+    h10002s(nData) = h10002;
+    h22000s(nData) = h22000;
+    h11110s(nData) = h11110;
+    h00220s(nData) = h00220;
+    h31000s(nData) = h31000;
+    h40000s(nData) = h40000;
+    h20110s(nData) = h20110;
+    h11200s(nData) = h11200;
+    h20020s(nData) = h20020;
+    h20200s(nData) = h20200;
+    h00310s(nData) = h00310;
+    h00400s(nData) = h00400;
+end
+

atmat/atphysics/NonLinearDynamics/RDTfluctuationIndicator.m

@@ -0,0 +1,56 @@
+function [h3ave,h4ave,h3chroave,adts]=RDTfluctuationIndicator(ring,varargin)
+%RDTFLUCTUATIONINDICATOR quantitative representation of RDT fluctuations
+%   This function calls computeRDTfluctuation(ring, varargin) to compute
+%   RDT fluctuations, and provides one example to quantitatively represents 
+%   the RDT fluctuations.
+%
+% [h3ave,h4ave,h3chroave,adts]=RDTfluctuationIndicator(ring,varargin)
+%
+%   ring is the AT lattice
+%   The additional argument:
+%     nslices: number of slices of each sextupole, which affects the
+%              crossing terms. default: 4.
+% 
+%   h3ave and h4ave quantitatively represents 3rd- and 4th-order geometric 
+%     RDT fluctuations, respectively. 
+%
+%   h3ave + w * h4ave can be used to represents geometric RDT fluctuations.
+%     The coefficient w can be estimated by action of particle at the
+%     anticipated DA, w ~ (2J_x)^0.5 = x / betax^0.5,  usually 0.01 is OK.
+% 
+%   h3chroave is the fluctuation of 3rd-order chromatic RDTs, 
+%     defined similarly to h3ave. 
+%
+%   adts is the sum of three ADTS terms, which are also used in onlinear optimization.
+%   The fluctuations of ADTS terms are not considered.
+%   It is calculated here to avoid duplicate computation.
+%
+% Noting:
+%    1.This function provides one example to quantitatively represents the 
+%      RDT fluctuations similar to that in Ref.[1]. But in Ref.[1], only
+%      the RDTs at the locations of nonlinear magnets are considered. 
+%      We think the differences are small and it's more important to 
+%      keep the function simple.
+%    2.People can call computeRDTfluctuation(ring, varargin) directly,
+%      and try other quantitative representations.
+%    3.The build-up RDT fluctuation can also be used, see Ref.[1].
+%      If the build-up RDT fluctuations are used, it is better to calculate
+%      the RDT build-up fluctuations for multiple periods to have
+%      better convergence of calculation.
+%
+% REFERENCE:
+%   [1] B. Wei, Z. Bai, J. Tan, L. Wang, and G. Feng, Phys. Rev. Accel. Beams 26, 084001 (2023)
+%
+
+[RDT, ~, natural] = computeRDTfluctuation(ring, 'nperiods', 1, varargin);
+h3ave = sqrt(mean(natural.f21000)^2 + mean(natural.f30000)^2 ...
+           + mean(natural.f10110)^2 + mean(natural.f10020)^2 ...
+           + mean(natural.f10200)^2);
+h3chroave = sqrt(mean(natural.f20001)^2 + mean(natural.f00201)^2 ...
+           + mean(natural.f10002)^2);
+h4ave = sqrt(mean(natural.f31000)^2 + mean(natural.f40000)^2 ...
+           + mean(natural.f20110)^2 + mean(natural.f11200)^2 ...
+           + mean(natural.f20020)^2 + mean(natural.f20200)^2 ...
+           + mean(natural.f00310)^2 + mean(natural.f00400)^2);
+adts = sqrt(RDT.dnux_dJx^2 + RDT.dnux_dJy^2 + RDT.dnuy_dJy^2);
+end