VariableMultipole compatibility between matlab AT and pyat

atintegrators/interpolate.c

@@ -2,12 +2,12 @@
    Accelerator Toolbox 
    Created: 14/11/08
    Z.Martí
- 
-  
-  
- Numerical recipies functions for 2D interpolation adapted to Matlab 
- * leanguagge
- 
+
+
+
+ Numerical recipies functions for 2D interpolation adapted to Matlab
+ * language
+
 */
 
 #include "at.h"

atmat/Contents.m

@@ -378,7 +378,7 @@
 %   atskewquad                          - Creates a skew quadrupole element with Class 'Multipole'
 %   atsolenoid                          - Creates a new solenoid element with Class 'Solenoid'
 %   atthinmultipole                     - Creates a thin multipole element
-%   atvariablemultipole                 - Creates a variable thin multipole element
+%   atvariablethinmultipole             - Creates a variable thin multipole element
 %   atwiggler                           - Creates a wiggler
 %   
 %   ATMAT/LATTICE/ELEMENT_CREATION/PRIVATE

atmat/atphysics/Radiation/atdiffmat.m

@@ -32,6 +32,16 @@
     orbit=num2cell(linepass(ring, orb0, 1:NumElements),1)';
 end
 
+% check for non-supported pass methods
+nonsupported_methods = {'VariableThinMultipole'};
+for i = 1:length(nonsupported_methods)
+    thepassm = nonsupported_methods{i};
+    lidx = find(atgetcells(ring, 'Class', thepassm));
+    if ~isempty(lidx)
+        ring = atsetfieldvalues(ring, lidx, 'PassMethod', 'IdentityPass');
+    end
+end
+
 zr=zeros(6,6);
 % Calculate cumulative Radiation-Diffusion matrix for the ring
 BCUM = zeros(6,6);

atmat/atutils/atinspectvariablethinmultipole.m

@@ -0,0 +1,166 @@
+function [listpola,listpolb] = atinspectvariablethinmultipole(element,varargin)
+% [pola_array,polb_array] = ATINSPECTVARIABLETHINMULTIPOLE(element)
+%
+% Get the ATVARIABLETHINMULTIPOLE polynom values per turn.
+%
+% Translations (T1,T2) and rotations (R1,R2) in the element are ignored.
+%
+% ARGUMENTS:
+%    element: an atvariablethinmultipole element.
+% OPTIONAL ARGUMENTS:
+%    turns(int): Default 1. Number of turns to calculate.
+%    T0(float): revolution time in seconds. Only needed in SINE mode.
+%    tparticle(float): Default 0. Time of the particle in seconds.
+%
+% RETURNS
+%    Two cell arrays with the values of PolynomA and PolynomB per turn.
+%
+% See Also:
+%   ATVARIABLETHINMULTIPOLE
+%
+    p = inputParser;
+    addOptional(p,'turns',1,@(x) isnumeric(x) && isscalar(x) && (x > 0));
+    if element.Mode == 0
+        if ~any(strcmp(varargin,'T0'))
+            error('T0 is required.');
+        end
+        addOptional(p,'T0',0,@(x) isnumeric(x) && (x > 0));
+    end
+    addOptional(p,'tparticle',0,@(x) isnumeric(x));
+    parse(p,varargin{:});
+    par = p.Results;
+
+    turns = par.turns;
+    if element.Mode == 0
+        T0 = par.T0;
+    else
+        T0 = 0;
+    end
+    tparticle = par.tparticle;
+    mode = element.Mode;
+
+    switch mode
+      case 0
+        % revolution time
+        timeoffset = T0 + tparticle;
+      case 2
+        % particle time
+        timeoffset = tparticle;
+      otherwise
+        timeoffset = 0;
+    end
+
+    if isfield(element,'Ramps')
+      ramps = element.Ramps;
+    else
+      ramps = 0;
+    end
+
+    if isfield(element,'Periodic')
+      periodic = element.Periodic;
+    else
+      periodic = 0;
+    end
+
+    maxorder = element.MaxOrder;
+
+    pola = nan(1,maxorder+1);
+    polb = nan(1,maxorder+1);
+
+    listpola = cell(1, turns);
+    listpolb = cell(1, turns);
+
+    for turn = 0:turns-1
+      for order = 0:maxorder
+        if isfield(element,'AmplitudeA')
+          pola(order+1) = get_pol(element, 'A', ramps, mode,  ...
+              timeoffset * turn, turn, order, periodic);
+        end
+        if isfield(element,'AmplitudeB')
+          pola(order+1) = get_pol(element, 'B', ramps, mode, ...
+              timeoffset * turn, turn, order, periodic);
+        end
+      listpola{turn+1} = pola;
+      listpolb{turn+1} = polb;
+      end
+    end
+
+end
+
+function ampt = get_amp(amp, ramps, t)
+%    get_amp(ele, amp, ramps, t)
+%
+%    get_amp returns the input value `amp` when ramps is False.
+%    If ramps is True, it returns a value linearly interpolated
+%    accoding to the ramping turn.
+    ampt = amp;
+    if length(ramps) == 4
+        if t <= ramps(0)
+            ampt = 0.0;
+        elseif t <= ramps(1)
+            ampt = amp * (t - ramps(0)) / (ramps(1) - ramps(0));
+        elseif t <= ramps(2)
+            ampt = amp;
+        elseif t <= ramps(3)
+            ampt = amp - amp * (t - ramps(2))/ (ramps(3) - ramps(2));
+        else
+            ampt = 0.0;
+        end
+    end
+end
+
+
+function ampout = get_pol(element , ab, ramps, mode, t, turn, order, periodic)
+%   get_pol(element , ab, ramps, mode, t, turn, order, periodic)
+%
+%   get_pol returns the polynom a or b for a given mode, turn, order,
+%   time and periodicity.
+    allamp = element.(strcat("Amplitude",ab) );
+    amp = allamp(order+1);
+    ampout = 0;
+
+    if amp == 0
+        return
+    end
+
+    % get the ramp value
+    ampoutaux = get_amp(amp, ramps, turn);
+
+    switch mode 
+      case 0
+        % sin mode parameters
+        whole_sin_above = element.(strcat("Sin",ab,"above"));
+        freq = element.(strcat("Frequency",ab));
+        ph = element.(strcat("Phase",ab));
+        sinval = sin(2 * pi * freq * t + ph);
+        if sinval >= whole_sin_above
+            ampout = ampoutaux * sinval;
+        else
+            ampout = 0;
+        end
+      case 1
+        ampout = nan;
+      case 2
+        nsamples = element.(strcat("NSamples",ab));
+        if periodic || turn < nsamples
+            func = element.(strcat("Func",ab));
+            funcdelay = element.(strcat("Func",ab,"TimeDelay"));
+            [ktaylor, nsamples] = size(func);
+            turnidx = mod(turn, nsamples)+1;
+            t = t - funcdelay;
+            functot = func(1,turnidx);
+            thefactorial = 1;
+            tpow = 1;
+            for i = 1:(ktaylor-1)
+                tpow = t *tpow;
+                thefactorial = thefactorial *i;
+                functot = functot + tpow/thefactorial * func(i+1,turnidx);
+            end
+            ampout = ampoutaux * functot;
+        else
+            ampout = 0.0;
+        end
+      otherwise
+        ampout = 0.0;
+    end
+end

atmat/lattice/at2str.m

@@ -119,6 +119,12 @@
                             'xtable', ...
                             'ytable' ...
                             });
+    case 'VariableThinMultipole'
+        create=@atvariablethinmultipole;
+        [options,args]=doptions(elem,create,{'Mode'});
+        if isfield(options,'Energy')
+            options=rmfield(options,'Energy');
+        end
     otherwise %'Marker'
         create=@atmarker;
         [options,args]=doptions(elem,create);

atmat/lattice/atguessclass.m

@@ -8,7 +8,7 @@
 %
 %  NOTES
 %  1. atclass=atguessclass(atelem,'useclass')
-%  By default, ATGUESSCLASS will default "useless" elements (PolynopmB==0)
+%  By default, ATGUESSCLASS will default "useless" elements (PolynomB==0)
 %  to 'Drift' or 'Marker', depending on 'Length'. When specifying
 %  'UseClass', it it will preserve the 'Class' field for those elements.
 %
@@ -30,6 +30,8 @@
     atclass='RingParam';
 elseif isfield(elem,'Limits')
     atclass='Aperture';
+elseif isfield(elem,'Mode')
+    atclass='VariableThinMultipole';
 elseif isfield(elem,'PolynomB')
     if useclass && isfield(elem,'Class')
         atclass=elem.Class;