diff --git a/atintegrators/VariableThinMPolePass.c b/atintegrators/VariableThinMPolePass.c
index 8ba2f88fa..a9d0e3da1 100644
--- a/atintegrators/VariableThinMPolePass.c
+++ b/atintegrators/VariableThinMPolePass.c
@@ -3,32 +3,78 @@
    S.White
 */
 
+/* 2024nov19 oblanco at ALBA CELLS. modified to be compatible with
+                                    pyat and AT matlab */
+
+/*
+ This pass method implements three modes.
+ Mode 0 : a sin function with frequency and phase
+ Mode 1 : a random kick
+ Mode 2 : a custom function defined in every turn
+
+ All modes could be ramped using the flag `ramp:.
+
+ The value of the polynoms A and B are calculated per turn during tracking,
+ and also, in modes 0 and 2, per particle to take in to account
+ the individual delay time.
+
+ In mode 0, the sin function could be limited to be between any value above
+ `Sinabove`. For example a positive half-sin function would be obtained by
+ setting `Sinabove` to zero.
+
+ In mode 1 the stream of pseudo-random values is taken from the
+ parameters structure in attypes.h. For more details about the random
+ generator in AT:
+ https://github.com/atcollab/at/discussions/879
+
+ In mode 2, the values could be periodically applied or not.
+*/
+
+
 #include "atconstants.h"
 #include "atelem.c"
 #include "atlalib.c"
 #include "atrandom.c"
 #include "driftkick.c"
 
+/* This struct contains the values to set one of the two
+   polynoms: A or B , as a function of the element parameters. */
 struct elemab {
     double* Amplitude;
     double Frequency;
+    double Sinabove;
     double Phase;
-    int NSamples;
     double* Func;
+    double FuncTimeDelay;
+    int NSamples;
+    int Ktaylor;
+    double* Buffer;
+    int BufferSize;
 };
 
+/* This struct contains the parameters of the element.
+ * It uses elemab struct */
 struct elem {
-    double* PolynomA;
-    double* PolynomB;
     struct elemab* ElemA;
     struct elemab* ElemB;
-    int Seed;
-    int Mode;
+    double* PolynomA;        // calculated on every turn
+    double* PolynomB;        // calculated on every turn
     int MaxOrder;
+    int Mode;
     double* Ramps;
     int Periodic;
+    double *R1;
+    double *R2;
+    double *T1;
+    double *T2;
+    double *EApertures;
+    double *RApertures;
 };
 
+
+/* 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.*/
 double get_amp(double amp, double* ramps, double t)
 {
     double ampt = amp;
@@ -48,32 +94,88 @@ double get_amp(double amp, double* ramps, double t)
     return ampt;
 }
 
-double get_pol(struct elemab* elem, double* ramps, int mode,
-    double t, int turn, int seed, int order, int periodic)
+
+/* get_pol calculates the PolynomA/B per turn and mode.
+   If the mode is 0 or 2, the polynom depends on the particle
+   time delay*/
+double get_pol(
+    struct elemab* elem,
+    double* ramps,
+    int mode,
+    double t,           // time
+    int turn,
+    int order,          //max(length of any of the polynoms) - 1
+    int periodic,
+    pcg32_random_t* rng
+    )
 {
-    int idx;
-    double ampt, freq, ph, val;
+    int turnidx,i; // index
+    double ampt;
+
+    // sin mode parameters
+    double whole_sin_above = elem->Sinabove;
+    double freq, ph, sinval;
+
+    // random mode parameters
+    double random_value;
+
+    // custom mode parameters
     double* func;
     double* amp = elem->Amplitude;
+    double functdelay;
+    double functot;
+    double tpow,thefactorial; // variables in for cycles
+    int Ktaylor;
+
     if (!amp) {
         return 0.0;
     }
     ampt = get_amp(amp[order], ramps, turn);
+
     switch (mode) {
     case 0:
         freq = elem->Frequency;
         ph = elem->Phase;
-        ampt *= sin(TWOPI * freq * t + ph);
+        sinval = sin(TWOPI * freq * t + ph);
+        // branchless if
+        // sinval >= wholeabove  -> sinval
+        // else                  -> 0
+        ampt = ampt*sinval*(sinval >= whole_sin_above);
         return ampt;
     case 1:
-        val = atrandn(0.0, 1.0);
-        ampt *= val;
+        random_value = atrandn_r(rng, 0, 1);
+        ampt *= random_value;
+        /* save random value into buffer */
+        if (turn < elem->BufferSize){
+            elem->Buffer[turn] = random_value;
+        }
         return ampt;
     case 2:
         if (periodic || turn < elem->NSamples) {
+            /* get the function, the delay, and turns */
             func = elem->Func;
-            idx = turn % elem->NSamples;
-            ampt *= func[idx];
+            functdelay = elem->FuncTimeDelay;
+            turnidx = turn % elem->NSamples;
+            Ktaylor = elem->Ktaylor;
+
+            /* calculate the amplitude as a Taylor expansion */
+            /* first order is taken directly from the function table */
+            t = t - functdelay;
+            functot = func[Ktaylor*turnidx];
+            tpow = 1;
+            thefactorial = 1;
+            /* do a taylor expansion if Ktaylor is more than first order */
+            for (i=1;i<Ktaylor;i++){
+              tpow = tpow * t;
+              thefactorial = thefactorial * i;
+              /* indexing is fortran-like. We start with columns.
+                 cols are taylor components.
+                 rows are turn samples.
+              */
+              functot =  functot + tpow / thefactorial * func[i + Ktaylor*turnidx];
+            };
+            ampt = ampt * functot;
+
             return ampt;
         } else {
             return 0.0;
@@ -83,41 +185,93 @@ double get_pol(struct elemab* elem, double* ramps, int mode,
     }
 }
 
-void VariableThinMPolePass(double* r, struct elem* Elem, double t0, int turn, int num_particles)
+/* This function tracks a particle through a thin element with
+   variable PolynomB and PolynomA per turn.*/
+void VariableThinMPolePass(
+    double* r, // the particle in 6D
+    struct elem* Elem,  // the variable element
+    double t0, // time of one revolution
+    int turn,  // number of the current turn
+    int num_particles, //number of particles
+    pcg32_random_t* rng // pcg32 random stream
+    )
 {
 
-    int i, c;
-    double* r6;
-    double t = t0 * turn;
+    int i; // order of the polynom, counter
+    int c; // particle, counter
+    double* r6; // particle 6D coordinates
+    double tpart; // particle time delay
 
-    int maxorder = Elem->MaxOrder;
-    int periodic = Elem->Periodic;
-    double* pola = Elem->PolynomA;
-    double* polb = Elem->PolynomB;
-    int seed = Elem->Seed;
+    // particle time offset for the mode
+    double time_in_this_mode = 0;
+
+    // setting the element properties
     int mode = Elem->Mode;
     struct elemab* ElemA = Elem->ElemA;
     struct elemab* ElemB = Elem->ElemB;
+    double* pola = Elem->PolynomA;
+    double* polb = Elem->PolynomB;
+    int maxorder = Elem->MaxOrder;
     double* ramps = Elem->Ramps;
 
-    if (mode != 0) {
+    // custom function is periodic
+    int periodic = Elem->Periodic;
+
+    // offsets at input and output
+    double *T1 = Elem->T1;
+    double *T2 = Elem->T2;
+    // rotations at input and output
+    double *R1 = Elem->R1;
+    double *R2 = Elem->R2;
+    // apertures
+    double *RApertures = Elem->RApertures;
+    double *EApertures = Elem->EApertures;
+
+    /* mode 0 : sin function */
+    /* mode 1 : random value applied to all particles */
+    /* mode 2 : custom function */
+
+    if (mode == 1) {
+         /* calculate the polynom to apply on all particles */
+         /* the particle delay time is not used */
         for (i = 0; i < maxorder + 1; i++) {
-            pola[i] = get_pol(ElemA, ramps, mode, t, turn, seed, i, periodic);
-            polb[i] = get_pol(ElemB, ramps, mode, t, turn, seed, i, periodic);
+            pola[i] = get_pol(ElemA, ramps, mode, 0, turn, i, periodic, rng);
+            polb[i] = get_pol(ElemB, ramps, mode, 0, turn, i, periodic, rng);
         };
     };
 
+    // offset the time when applying the sin function
+    // branchless if
+    // mode == 0 -> time_in_this_mode = t0 *turn
+    // else      -> 0
+    time_in_this_mode = t0 * turn * (mode == 0);
+
+    /* cycle over all particles */
     for (c = 0; c < num_particles; c++) {
         r6 = r + c * 6;
+        /* check if the particle is alive */
         if (!atIsNaN(r6[0])) {
-            if (mode == 0) {
-                double tpart = t + r6[5] / C0;
+            /* mode 0 and mode 2 take into account the particle delay time */
+            if (mode == 0 || mode == 2) {
+                /* modify the time of delay of the particle  */
+                tpart = time_in_this_mode + r6[5] / C0;
+                /* calculate the polynom A and B components seen by the particle */
                 for (i = 0; i < maxorder + 1; i++) {
-                    pola[i] = get_pol(ElemA, ramps, mode, tpart, turn, seed, i, periodic);
-                    polb[i] = get_pol(ElemB, ramps, mode, tpart, turn, seed, i, periodic);
+                    pola[i] = get_pol(ElemA, ramps, mode, tpart, turn, i, periodic, rng);
+                    polb[i] = get_pol(ElemB, ramps, mode, tpart, turn, i, periodic, rng);
                 };
             };
+            /*  misalignment at entrance  */
+            if (T1) ATaddvv(r6,T1);
+            if (R1) ATmultmv(r6,R1);
+            /* Check physical apertures at the entrance of the magnet */
+            if (RApertures) checkiflostRectangularAp(r6,RApertures);
+            if (EApertures) checkiflostEllipticalAp(r6,EApertures);
+            /* track */
             strthinkick(r6, pola, polb, 1.0, maxorder);
+            /* Misalignment at exit */
+            if (R2) ATmultmv(r6,R2);
+            if (T2) ATaddvv(r6,T2);
         }
     }
 }
@@ -126,13 +280,24 @@ void VariableThinMPolePass(double* r, struct elem* Elem, double t0, int turn, in
 ExportMode struct elem* trackFunction(const atElem* ElemData, struct elem* Elem,
     double* r_in, int num_particles, struct parameters* Param)
 {
+    /* Initialize element */
     if (!Elem) {
-        int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic;
+        int MaxOrder, Mode, NSamplesA, NSamplesB, KtaylorA, KtaylorB, Periodic;
+        int BufferSizeA, BufferSizeB;
+        double *R1, *R2, *T1, *T2, *EApertures, *RApertures;
         double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB;
-        double *Ramps, *FuncA, *FuncB;
+        double *Ramps, *FuncA, *FuncB, *BufferA, *BufferB;
+        double FuncATimeDelay, FuncBTimeDelay;
         double FrequencyA, FrequencyB;
         double PhaseA, PhaseB;
+        double SinAabove, SinBabove;
         struct elemab *ElemA, *ElemB;
+        R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error();
+        R2=atGetOptionalDoubleArray(ElemData,"R2"); check_error();
+        T1=atGetOptionalDoubleArray(ElemData,"T1"); check_error();
+        T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error();
+        EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error();
+        RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error();
         MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error();
         Mode=atGetLong(ElemData,"Mode"); check_error();
         PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error();
@@ -143,20 +308,34 @@ ExportMode struct elem* trackFunction(const atElem* ElemData, struct elem* Elem,
         FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error();
         PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error();
         PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error();
+        SinAabove=atGetOptionalDouble(ElemData,"SinAabove", 0); check_error();
+        SinBabove=atGetOptionalDouble(ElemData,"SinBabove", 0); check_error();
         Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error();
-        Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error();
         NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 1); check_error();
         NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 1); check_error();
+        KtaylorA=atGetOptionalLong(ElemData, "KtaylorA", 1); check_error();
+        KtaylorB=atGetOptionalLong(ElemData, "KtaylorB", 1); check_error();
         FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error();
         FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error();
-        Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error();
+        BufferA=atGetOptionalDoubleArray(ElemData,"BufferA"); check_error();
+        BufferB=atGetOptionalDoubleArray(ElemData,"BufferB"); check_error();
+        BufferSizeA=atGetOptionalLong(ElemData, "BufferSizeA", 0); check_error();
+        BufferSizeB=atGetOptionalLong(ElemData, "BufferSizeB", 0); check_error();
+        FuncATimeDelay=atGetOptionalDouble(ElemData,"FuncATimeDelay", 0); check_error();
+        FuncBTimeDelay=atGetOptionalDouble(ElemData,"FuncBTimeDelay", 0); check_error();
+        Periodic=atGetOptionalLong(ElemData,"Periodic", 0); check_error();
         Elem = (struct elem*)atMalloc(sizeof(struct elem));
         ElemA = (struct elemab*)atMalloc(sizeof(struct elemab));
         ElemB = (struct elemab*)atMalloc(sizeof(struct elemab));
+        Elem->R1=R1;
+        Elem->R2=R2;
+        Elem->T1=T1;
+        Elem->T2=T2;
+        Elem->EApertures=EApertures;
+        Elem->RApertures=RApertures;
         Elem->PolynomA = PolynomA;
         Elem->PolynomB = PolynomB;
         Elem->Ramps = Ramps;
-        Elem->Seed = Seed;
         Elem->Mode = Mode;
         Elem->MaxOrder = MaxOrder;
         Elem->Periodic = Periodic;
@@ -166,16 +345,35 @@ ExportMode struct elem* trackFunction(const atElem* ElemData, struct elem* Elem,
         ElemB->Frequency = FrequencyB;
         ElemA->Phase = PhaseA;
         ElemB->Phase = PhaseB;
+        ElemA->Sinabove = SinAabove;
+        ElemB->Sinabove = SinBabove;
         ElemA->NSamples = NSamplesA;
         ElemB->NSamples = NSamplesB;
+        ElemA->Ktaylor = KtaylorA;
+        ElemB->Ktaylor = KtaylorB;
         ElemA->Func = FuncA;
         ElemB->Func = FuncB;
+        ElemA->Buffer = BufferA;
+        ElemB->Buffer = BufferB;
+        ElemA->BufferSize = BufferSizeA;
+        ElemB->BufferSize = BufferSizeB;
+        ElemA->FuncTimeDelay = FuncATimeDelay;
+        ElemB->FuncTimeDelay = FuncBTimeDelay;
         Elem->ElemA = ElemA;
         Elem->ElemB = ElemB;
     }
-    double t0 = Param->T0;
-    int turn = Param->nturn;
-    VariableThinMPolePass(r_in, Elem, t0, turn, num_particles);
+    double t0 = Param->T0;  // revolution time of the nominal ring
+    int turn = Param->nturn;  // current turn
+    int i;
+
+    /* track */
+    VariableThinMPolePass(r_in, Elem, t0, turn, num_particles, Param->thread_rng);
+    /* reset the polynom values.
+       This is done to save the ring in an unperturbed configuration. */
+    for (i = 0; i < Elem->MaxOrder + 1; i++){
+        Elem->PolynomA[i] = 0;
+        Elem->PolynomB[i] = 0;
+    };
     return Elem;
 }
 
@@ -190,14 +388,24 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
         double* r_in;
         const mxArray* ElemData = prhs[0];
         int num_particles = mxGetN(prhs[1]);
-        int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic;
+        int MaxOrder, Mode, NSamplesA, NSamplesB, KtaylorA, KtaylorB, Periodic;
+        int BufferSizeA, BufferSizeB;
+        double *R1, *R2, *T1, *T2, *EApertures, *RApertures;
         double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB;
-        double *Ramps, *FuncA, *FuncB;
+        double *Ramps, *FuncA, *FuncB, *BufferA, *BufferB;
+        double FuncATimeDelay, FuncBTimeDelay;
         double FrequencyA, FrequencyB;
         double PhaseA, PhaseB;
+        double SinAabove, SinBabove;
         struct elemab ElA, *ElemA = &ElA;
         struct elemab ElB, *ElemB = &ElB;
         struct elem El, *Elem = &El;
+        R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error();
+        R2=atGetOptionalDoubleArray(ElemData,"R2"); check_error();
+        T1=atGetOptionalDoubleArray(ElemData,"T1"); check_error();
+        T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error();
+        EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error();
+        RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error();
         MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error();
         Mode=atGetLong(ElemData,"Mode"); check_error();
         PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error();
@@ -208,17 +416,25 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
         FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error();
         PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error();
         PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error();
+        SinAabove=atGetOptionalDouble(ElemData,"SinAabove", 0); check_error();
+        SinBabove=atGetOptionalDouble(ElemData,"SinBabove", 0); check_error();
         Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error();
-        Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error();
         NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 0); check_error();
         NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 0); check_error();
+        KtaylorA=atGetOptionalLong(ElemData, "KtaylorA", 0); check_error();
+        KtaylorB=atGetOptionalLong(ElemData, "KtaylorB", 0); check_error();
         FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error();
         FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error();
-        Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error();
+        BufferA=atGetOptionalDoubleArray(ElemData,"BufferA"); check_error();
+        BufferB=atGetOptionalDoubleArray(ElemData,"BufferB"); check_error();
+        BufferSizeA=atGetOptionalLong(ElemData, "BufferSizeA", 0); check_error();
+        BufferSizeB=atGetOptionalLong(ElemData, "BufferSizeB", 0); check_error();
+        FuncATimeDelay=atGetOptionalDouble(ElemData,"FuncATimeDelay", 0); check_error();
+        FuncBTimeDelay=atGetOptionalDouble(ElemData,"FuncBTimeDelay", 0); check_error();
+        Periodic=atGetOptionalLong(ElemData,"Periodic", 0); check_error();
         Elem->PolynomA = PolynomA;
         Elem->PolynomB = PolynomB;
         Elem->Ramps = Ramps;
-        Elem->Seed = Seed;
         Elem->Mode = Mode;
         Elem->MaxOrder = MaxOrder;
         Elem->Periodic = Periodic;
@@ -228,16 +444,27 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
         ElemB->Frequency = FrequencyB;
         ElemA->Phase = PhaseA;
         ElemB->Phase = PhaseB;
+        ElemA->Sinabove = SinAabove;
+        ElemB->Sinabove = SinBabove;
         ElemA->NSamples = NSamplesA;
         ElemB->NSamples = NSamplesB;
+        ElemA->Ktaylor = KtaylorA;
+        ElemB->Ktaylor = KtaylorB;
         ElemA->Func = FuncA;
         ElemB->Func = FuncB;
+        ElemA->Buffer = BufferA;
+        ElemB->Buffer = BufferB;
+        ElemA->BufferSize = BufferSizeA;
+        ElemB->BufferSize = BufferSizeB;
         Elem->ElemA = ElemA;
         Elem->ElemB = ElemB;
+        ElemA->FuncTimeDelay = FuncATimeDelay;
+        ElemB->FuncTimeDelay = FuncBTimeDelay;
         /* ALLOCATE memory for the output array of the same size as the input  */
         plhs[0] = mxDuplicateArray(prhs[1]);
         r_in = mxGetDoubles(plhs[0]);
-        VariableThinMPolePass(r_in, Elem, 0, 0, num_particles);
+
+        VariableThinMPolePass(r_in, Elem, 0, 0, num_particles, &pcg32_global);
     } else if (nrhs == 0) {
         /* list of required fields */
         plhs[0] = mxCreateCellMatrix(4, 1);
@@ -245,25 +472,40 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
         mxSetCell(plhs[0], 1, mxCreateString("Mode"));
         mxSetCell(plhs[0], 2, mxCreateString("PolynomA"));
         mxSetCell(plhs[0], 3, mxCreateString("PolynomB"));
-        if (nlhs > 1) {
+        if (nlhs > 4) {
             /* list of optional fields */
-            plhs[1] = mxCreateCellMatrix(13, 1);
+            plhs[1] = mxCreateCellMatrix(28, 1);
             mxSetCell(plhs[1], 0, mxCreateString("AmplitudeA"));
             mxSetCell(plhs[1], 1, mxCreateString("AmplitudeB"));
             mxSetCell(plhs[1], 2, mxCreateString("FrequencyA"));
             mxSetCell(plhs[1], 3, mxCreateString("FrequencyB"));
             mxSetCell(plhs[1], 4, mxCreateString("PhaseA"));
             mxSetCell(plhs[1], 5, mxCreateString("PhaseB"));
-            mxSetCell(plhs[1], 6, mxCreateString("Ramps"));
-            mxSetCell(plhs[1], 7, mxCreateString("Seed"));
-            mxSetCell(plhs[1], 8, mxCreateString("FuncA"));
-            mxSetCell(plhs[1], 9, mxCreateString("FuncB"));
-            mxSetCell(plhs[1], 10, mxCreateString("NSamplesA"));
-            mxSetCell(plhs[1], 11, mxCreateString("NSamplesB"));
-            mxSetCell(plhs[1], 12, mxCreateString("Periodic"));
+            mxSetCell(plhs[1], 6, mxCreateString("SinAabove"));
+            mxSetCell(plhs[1], 7, mxCreateString("SinBabove"));
+            mxSetCell(plhs[1], 8, mxCreateString("Ramps"));
+            mxSetCell(plhs[1], 9, mxCreateString("FuncA"));
+            mxSetCell(plhs[1], 10,mxCreateString("FuncB"));
+            mxSetCell(plhs[1], 11,mxCreateString("BufferA"));
+            mxSetCell(plhs[1], 12,mxCreateString("BufferB"));
+            mxSetCell(plhs[1], 13,mxCreateString("BufferSizeA"));
+            mxSetCell(plhs[1], 14,mxCreateString("BufferSizeB"));
+            mxSetCell(plhs[1], 15,mxCreateString("FuncATimeDelay"));
+            mxSetCell(plhs[1], 16,mxCreateString("FuncBTimeDelay"));
+            mxSetCell(plhs[1], 17,mxCreateString("NSamplesA"));
+            mxSetCell(plhs[1], 18,mxCreateString("NSamplesB"));
+            mxSetCell(plhs[1], 19,mxCreateString("KtaylorA"));
+            mxSetCell(plhs[1], 20,mxCreateString("KtaylorB"));
+            mxSetCell(plhs[1], 21,mxCreateString("Periodic"));
+            mxSetCell(plhs[1], 22,mxCreateString("T1"));
+            mxSetCell(plhs[1], 23,mxCreateString("T2"));
+            mxSetCell(plhs[1], 24,mxCreateString("R1"));
+            mxSetCell(plhs[1], 25,mxCreateString("R2"));
+            mxSetCell(plhs[1], 26,mxCreateString("RApertures"));
+            mxSetCell(plhs[1], 27,mxCreateString("EApertures"));
         }
     } else {
-        mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments");
+        mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 4 arguments");
     }
 }
 #endif /* MATLAB_MEX_FILE */
diff --git a/atintegrators/interpolate.c b/atintegrators/interpolate.c
index 09dc5f9da..9284757f7 100644
--- a/atintegrators/interpolate.c
+++ b/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"
diff --git a/atmat/Contents.m b/atmat/Contents.m
index 534c2687f..2ad28ba6c 100644
--- a/atmat/Contents.m
+++ b/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
diff --git a/atmat/atphysics/Radiation/atdiffmat.m b/atmat/atphysics/Radiation/atdiffmat.m
index 8546cb39d..bea31808a 100644
--- a/atmat/atphysics/Radiation/atdiffmat.m
+++ b/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);
diff --git a/atmat/atutils/atinspectvariablethinmultipole.m b/atmat/atutils/atinspectvariablethinmultipole.m
new file mode 100644
index 000000000..4b272e3ae
--- /dev/null
+++ b/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
diff --git a/atmat/lattice/at2str.m b/atmat/lattice/at2str.m
index d1dbb7406..516cfa97d 100644
--- a/atmat/lattice/at2str.m
+++ b/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);
diff --git a/atmat/lattice/atguessclass.m b/atmat/lattice/atguessclass.m
index 36d22c0b8..40aa5df0c 100644
--- a/atmat/lattice/atguessclass.m
+++ b/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;
diff --git a/atmat/lattice/element_creation/atvariablemultipole.m b/atmat/lattice/element_creation/atvariablemultipole.m
deleted file mode 100644
index bdf83b4d7..000000000
--- a/atmat/lattice/element_creation/atvariablemultipole.m
+++ /dev/null
@@ -1,129 +0,0 @@
-function elem=atvariablemultipole(fname,varargin)
-%ATVARIABLEMULTIPOLE Creates a variable thin multipole element
-%
-%  ATVARIABLEMULTIPOLE(FAMNAME,MODE,PASSMETHOD,[KEY,VALUE]...)
-%	
-%  INPUTS
-%    FNAME          Family name 
-%    MODE           Excitation mode: 'SINE', 'WHITENOISE' or 'ARBITRARY'.
-%                   Default: 'SINE'
-%    PASSMETHOD     Tracking function. Default: 'VariableThinMPolePass'
-%
-%  OPTIONS (order does not matter)
-%    AMPLITUDEA     Vector or scalar to define the excitation amplitude for
-%                   PolynomA
-%    AMPLITUDEB     Vector or scalar to define the excitation amplitude for
-%                   PolynomA
-%    FREQUENCYA     Frequency of SINE excitation for PolynomA
-%    FREQUENCYB     Frequency of SINE excitation for PolynomB
-%    PHASEA         Phase of SINE excitation for PolynomA
-%    PHASEB         Phase of SINE excitation for PolynomB
-%	 MAXORDER       Order of the multipole for a scalar amplitude
-%    SEED           Input seed for the random number generator
-%    FUNCA          ARBITRARY excitation turn-by-turn kick list for PolynomA
-%    FUNCB          ARBITRARY excitation turn-by-turn kick list for PolynomB
-%    PERIODIC       If true (default) the user input kick list is repeated
-%    RAMPS          Vector (t0, t1, t2, t3) in turn number to define the ramping of the excitation
-%                   * t<t0: excitation amlpitude is zero
-%                   * t0<t<t1: exciation amplitude is linearly ramped up
-%                   * t1<t<t2: exciation amplitude is constant             
-%                   * t2<t<t3: exciation amplitude is linearly ramped down
-%                   * t3<t: exciation amplitude is zero 
-%
-%  OUTPUTS
-%  1. ELEM - Structure with the AT element
-%
-%  NOTES
-%    1. For all excitation modes at least one amplitude (A or B) is
-%    required. The default excitation is SINE
-%    2. For SINE excitation modes the FREQUENCY corresponding to the input
-%    AMPLITUDE is required
-%    3. For ARBITRARY excitation modes the FUNC corresponding to the input
-%    AMPLITUDE is required
-%
-%  EXAMPLES
-%
-% % Create a sinusoidal dipole with amplitude 0.1 mrad and frequency 1 kHz
-% >> atvariablemultipole('ACM','SINE','AmplitudeB',1.e-4,'FrequencyB',1.e3);
-%
-% % Create a white noise dipole excitation of amplitude 0.1 mrad
-% >> atvariablemultipole('ACM','WHITENOISE','AmplitudeB',1.e-4);
-
-% Input parser for option
-[mode,rsrc]=getargs(varargin,'SINE','check',@(arg) any(strcmpi(arg,{'SINE','WHITENOISE','ARBITRARY'})));
-[method,rsrc]=getargs(rsrc,'VariableThinMPolePass','check',@(arg) (ischar(arg) || isstring(arg)) && endsWith(arg,'Pass'));
-[mode,rsrc]                       = getoption(rsrc,'Mode',mode);
-[method,rsrc]                     = getoption(rsrc,'PassMethod',method);
-[cl,rsrc]                         = getoption(rsrc,'Class','VariableMultipole');
-[maxorder,rsrc]                   = getoption(rsrc,'MaxOrder',0);
-rsrc                              = struct(rsrc{:});
-rsrc.MaxOrder                     = maxorder;
-
-m=struct('SINE',0,'WHITENOISE',1,'ARBITRARY',2);
-
-if ~any(isfield(rsrc,{'AmplitudeA','AmplitudeB'}))
-    error("Please provide at least one amplitude for A or B")
-end
-rsrc = setparams(rsrc,mode,'A');
-rsrc = setparams(rsrc,mode,'B');
-rsrc = setmaxorder(rsrc);
-
-% Build the element
-% rsrc =namedargs2cell(rsrc);   % introduced in R2019b
-rsrc=reshape([fieldnames(rsrc) struct2cell(rsrc)]',1,[]);
-elem=atbaselem(fname,method,'Class',cl,'Length',0,'Mode',m.(upper(mode)),...
-               'PolynomA',[],'PolynomB',[],rsrc{:});
-
-
-    function setsine(rsrc, ab)
-        if ~isfield(rsrc,strcat('Frequency',ab))
-            error(strcat('Please provide a value for Frequency',ab))
-        end
-    end
-
-    function rsrc = setarb(rsrc, ab)
-        funcarg=strcat('Func',ab);
-        if ~isfield(rsrc,funcarg)
-            error(strcat('Please provide a value for Func',ab))
-        end
-        rsrc.(strcat('NSamples',ab))=length(rsrc.(funcarg));
-    end
-
-    function rsrc = setparams(rsrc,mode,ab)
-        amplarg=strcat('Amplitude',ab);
-        if isfield(rsrc,amplarg)
-            amp=rsrc.(amplarg);
-            if isscalar(amp)
-                rsrc.(amplarg)=[zeros(1,rsrc.MaxOrder) amp];
-            end
-            if strcmpi(mode,'SINE')
-                setsine(rsrc,ab);
-            end
-            if strcmpi(mode,'ARBITRARY')
-                rsrc = setarb(rsrc,ab);
-            end        
-        end
-    end
-
-    function rsrc = setmaxorder(rsrc)
-        if isfield(rsrc,'AmplitudeA')
-            mxa=find(abs(rsrc.AmplitudeA)>0,1,'last');
-        else
-            mxa=0;
-        end
-        if isfield(rsrc,'AmplitudeB')
-            mxb=find(abs(rsrc.AmplitudeB)>0,1,'last');
-        else
-            mxb=0;
-        end
-        mxab=max([mxa,mxb,rsrc.MaxOrder-1]);
-        rsrc.MaxOrder=mxab-1;
-        if isfield(rsrc,'AmplitudeA')
-            rsrc.AmplitudeA(mxa+1:mxab)=0;
-        end
-        if isfield(rsrc,'AmplitudeB')
-            rsrc.AmplitudeB(mxb+1:mxab)=0;
-        end               
-    end
-    
-end
diff --git a/atmat/lattice/element_creation/atvariablethinmultipole.m b/atmat/lattice/element_creation/atvariablethinmultipole.m
new file mode 100644
index 000000000..b3aa7faf4
--- /dev/null
+++ b/atmat/lattice/element_creation/atvariablethinmultipole.m
@@ -0,0 +1,309 @@
+function elem=atvariablethinmultipole(fname, inmode, varargin)
+% ATVARIABLETHINMULTIPOLE Creates a variable thin multipole element
+%
+%  ATVARIABLETHINMULTIPOLE(FAMNAME,MODE)
+%  ATVARIABLETHINMULTIPOLE(FAMNAME,MODE,PASSMETHOD,[KEY,VALUE]...)
+%
+%  INPUTS
+%    FNAME          Family name 
+%    MODE           Excitation mode: 'SINE', 'WHITENOISE' or 'ARBITRARY'.
+%
+%  OPTIONS
+%    PASSMETHOD     Tracking function. Default: 'VariableThinMPolePass'
+%
+%  MORE OPTIONS (order does not matter)
+%    AmplitudeA     Vector or scalar to define the excitation amplitude for
+%                   PolynomA
+%    AmplitudeB     Vector or scalar to define the excitation amplitude for
+%                   PolynomB
+%    FrequencyA     Frequency of SINE excitation for PolynomA. Unit Hz
+%    FrequencyB     Frequency of SINE excitation for PolynomB. Unit Hz.
+%    PhaseA         Phase of SINE excitation for PolynomA
+%    PhaseB         Phase of SINE excitation for PolynomB
+%    SinAabove      Limit the sine function to be above. Default -1.1
+%    SinBabove      Limit the sine function to be above. Default -1.1
+%    BufferSizeA    Set the buffer length in WHITENOISE mode.
+%    BufferSizeB    Set the buffer length in WHITENOISE mode.
+%    FuncA          ARBITRARY excitation turn-by-turn (tbt) list for
+%                   PolynomA
+%    FuncB          ARBITRARY excitation turn-by-turn (tbt) list for
+%                   PolynomB
+%    FuncATimeDelay TimeDelay to generate a small time offset on the
+%                   function FUNC. It only has an effect if any of the
+%                   derivatives is not zero. Default: 0.
+%    FuncBTimeDelay TimeDelay to generate a small time offset on the
+%                   function FUNC. It only has an effect if any of the
+%                   derivatives is not zero. Default 0.
+%    Periodic       If 1 the user input list is repeated. Default: 0.
+%    Ramps          Vector (t0, t1, t2, t3) in turn number to define the
+%                   ramping of the excitation:
+%                   * t<=t0: excitation amplitude is zero.
+%                   * t0<t<=t1: excitation amplitude is linearly ramped up.
+%                   * t1<t<=t2: excitation amplitude is constant.
+%                   * t2<t<=t3: excitation amplitude is linearly ramped down.
+%                   * t3<t:    excitation amplitude is zero.
+%
+%  OUTPUTS
+%  1. ELEM - Structure with the AT element
+%
+%  NOTES
+%    1. For all excitation modes at least one amplitude (A or B) is
+%    required.
+%    2. For SINE excitation modes the FREQUENCY is required.
+%    3. For ARBITRARY excitation modes the FUNC is required.
+%
+%  EXAMPLES
+%
+% % Create a sinusoidal skew quadrupole with amplitude 0.0001  and frequency 1 kHz
+% >> atvariablethinmultipole('ACSKEW','SINE','AmplitudeA',[0,1.e-4],'FrequencyA',1.e3);
+%
+% % Create a white noise horizontal dipole excitation of amplitude 0.1 mrad
+% >> atvariablethinmultipole('ACKICK','WHITENOISE','AmplitudeB',1.e-4);
+%
+% % Create a vertical kick in the first turn and the opposite kick in the second
+% % turn.
+% >> funca = [1 -1 0];
+% >> atvariablethinmultipole('CUSTOMFUNC','ARBITRARY','AmplitudeA',1e-4, ...
+%                     'FuncA',funca,'Periodic', 1);
+%
+%
+% MORE DETAILS
+%
+% This Class creates a thin multipole of any order (dipole kick, quadrupole,
+% sextupole, etc.) and type (Normal or Skew) defined by AmplitudeA and/or
+% AmplitudeB components; the polynoms PolynomA and PolynomB are calculated
+% on every turn depending on the chosen mode, and for some modes on the
+% particle time delay. All modes could be ramped.
+% Default pass method: ``VariableThinMPolePass``.
+%
+% Keep in mind that as this element varies on every turn, and at the end of
+% the tracking PolynomA and PolynomB are set to zero, this could interfere
+% with optics calculations. Therefore, it is preferred to set the pass method
+% to IdentityPass when the element should not be active.
+%
+% There are three different modes that could be set,
+% SINE: sine function
+% WHITENOISE: gaussian white noise
+% ARBITRARY: user defined turn-by-turn kick list
+% SINE = 0, WHITENOISE = 1 and ARBITRARY = 2, when the structure is
+% displayed, or in a file.
+%
+% For example, use SINE or mode = 0 to create an sine function element.
+%
+% The **SINE** mode requires amplitude and frequency for A and/or B.
+% The value of the jth component of the polynom (A or B) at the nth turn
+% is given by
+%
+% Amplitude[j]*sin[TWOPI*frequency*(n*T0 + \tau_p) + phase],
+%
+% where T0 is the revolution period of the ideal ring, and \tau_p is the delay
+% of the pth particle i.e. the sixth coordinate over the speed of light. Also,
+% note that the position of the element on the ring has no effect, the phase
+% could be used to add any delay due to the position along s. The following is
+% an example of the SINE mode of an skew quad
+%
+% eleskew = atvariablethinmultipole('VAR_SKEW',"SINE", ...
+% "AmplitudeA",[0,skewa2],"FrequencyA",freqA,"PhaseA",phaseA)
+%
+% The values of the sine function could be limited to be above a defined
+% threshold using ``Sin[AB]above``. For example, you could create a positive
+% half-sine by setting ``Sin[AB]above`` to zero. You could also create a
+% negative half-sine function by setting the amplitude to -1.
+%
+% The **WHITENOISE** mode requires the amplitude of A and/or B. The gaussian
+% distribution is generated with zero-mean and one standard deviation from
+% a pseudo-random stream pcg32. The pcg32 seed is fixed by the tracking
+% function, therefore using the same stream on all trackings (sequencial or
+% parallel). See https://github.com/atcollab/at/discussions/879 for more
+% details on the pseudo random stream. For example,
+%
+% elenoise = atvariablethinmultipole('MYNOISE',"WHITENOISE", ...
+% "AmplitudeA",[noiseA1])
+%
+% creates a vertical kick as gaussian noise of amplitude noiseA1.
+%
+% Additionally, one could use the parameter BufferSizeA or BufferSizeB
+% to create buffer to store the random numbers used in the element settings
+% during tracking. The buffer could be smaller or larger than the tracking.
+% If the buffer is smaller, it will store from start to the moment of filling
+% completely the buffer. If the buffer is larger, all non used places are set
+% to zero.
+%
+% Care should be taken when changing the buffer and buffer size. The only check
+% done is at the moment of initialization, for example, when reading a lattice
+% from a file. Buffer and buffer size could be replaced at any moment, just
+% ensure that the actual buffer length and buffersize agree.
+%
+% The **ARBITRARY** mode requires the definition of a custom discrete function
+% to be sampled at every turn. The function and its Taylor expansion with
+% respect to \tau up to any given order is
+%
+% value = f(turn) + f'(turn)*tau + 0.5*f''(turn)*tau**2
+% + 1/6*f'''(turn)*tau**3 + 1/24*f''''(turn)*tau**4 ...
+%
+% f is a row array of values, f',f'',f''',f'''', are row arrays containing
+% the derivatives wrt \tau, and \tau is the time delay of the particle, i.e.
+% the the sixth coordinate divided by the speed of light. Therefore, the
+% function func is a 2D-array with a given column contains f and its
+% derivatives, and the ith column has the values to be used in the ith turn.
+% For example, a positive value on the first turn f(1)=1 with positive
+% derivative f'(1)=0.1 followed by a negative value in the second turn
+% f(2)=-1 with negative derivative f'(2)=-0.2, and f(3) = 2, would be
+%   Func=[[1,  -1   2];
+%         [0.1 -0.2 0]].
+%
+% The time tau could be offset using ``FuncATimeDelay`` or ``FuncBTimeDelay``.
+%
+% tau = tau - Func[AB]TimeDelay
+%
+% The function value is then **multiplied by Amplitude A and/or B**.
+% Use FuncA or FuncB, and AmplitudeA or AmplitudeB accordingly.
+% For example, the following is a positive vertical kick in the first turn,
+% negative on the second turn, and zero on the third turn.
+%
+% elesinglekick = atvariablethinmultipole('CUSTOMFUNC',"ARBITRARY", ...
+% "AmplitudeA",1e-4,"FuncA",[1,-1,0],"Periodic",1)
+%
+% As already mentioned, one could include the derivatives of the function
+% from a Taylor expansion in a 2D array. First row for the values of f on
+% every turn, second row for the values of f' on every turn, and so on.
+%
+% By default the array is assumed non periodic, the function has no effect
+% on the particle in turns exceeding the function definition. If
+% ``Periodic`` is set to 1, the sequence is repeated.
+%
+% One could use the function atinspectvariablethinmultiple to check
+% the polynom values used in every turn before actually doing the tracking.
+%
+% Any mode could be ramped. The ramp is defined by four values (t0,t1,t2,t3):
+%     * ``t<=t0``: excitation amplitude is zero.
+%     * ``t0<t<=t1``: excitation amplitude is linearly ramped up.
+%     * ``t1<t<=t2``: excitation amplitude is constant.
+%     * ``t2<t<=t3``: excitation amplitude is linearly ramped down.
+%     * ``t3<t``: excitation amplitude is zero.
+% See also:
+%
+%  *ATINSPECTVARIABLETHINMULTIPOLE*
+
+% Input parser for option
+[method,rsrc]=getargs(varargin,'VariableThinMPolePass', ...
+    'check',@(arg) (ischar(arg) || isstring(arg)) && endsWith(arg,'Pass'));
+[method,rsrc] = getoption(rsrc,'PassMethod',method);
+[cl,rsrc] = getoption(rsrc,'Class','VariableThinMultipole');
+rsrc = struct(rsrc{:});
+
+m=struct('SINE',0,'WHITENOISE',1,'ARBITRARY',2);
+
+if ~any(isfield(rsrc,{'AmplitudeA','AmplitudeB'}))
+    rsrc.AmplitudeB = 0;
+    rsrc.AmplitudeA = 0;
+end
+
+% accept name or number mode
+if isnumeric(inmode)
+    mode = inmode;
+else
+    mode = m.(upper(inmode));
+end
+
+thepols = {'A','B'};
+for i = 1:length(thepols)
+    rsrc = setparams(rsrc,mode,thepols{i});
+    thefieldname = strcat('Polynom',thepols{i});
+    if ~isfield(rsrc, thefieldname)
+        rsrc.(thefieldname)=[];
+    end
+end
+rsrc = setmaxorder(rsrc);
+
+% Build the element
+% rsrc =namedargs2cell(rsrc);   % introduced in R2019b
+rsrc=reshape([fieldnames(rsrc) struct2cell(rsrc)]',1,[]);
+elem=atbaselem(fname,method,'Class',cl,'Length',0,'Mode',mode,rsrc{:});
+
+
+    function rsrc = setsine(rsrc, ab)
+        funcarg = strcat('Frequency',ab);
+        if ~isfield(rsrc,funcarg)
+            rsrc.(funcarg) = 0;
+        end
+        funcarg=strcat('Phase',ab);
+        if ~isfield(rsrc,funcarg)
+            rsrc.(funcarg) = 0;
+        end
+        funcarg=strcat('Sin',ab,'above');
+        if ~isfield(rsrc,funcarg)
+            rsrc.(funcarg) = -1.1;
+        end
+    end
+
+    function rsrc = setarb(rsrc, ab)
+        funcarg=strcat('Func',ab);
+        if ~isfield(rsrc,funcarg)
+            rsrc.(funcarg) = zeros(1,1);
+        end
+        [ktaylor, nsamples] = size(rsrc.(funcarg));
+        rsrc.(strcat('NSamples',ab)) = nsamples;
+        rsrc.(strcat('Ktaylor', ab)) = ktaylor;
+        funcarg = strcat('Func',ab,'TimeDelay');
+        if ~isfield(rsrc,funcarg)
+            rsrc.(funcarg) = 0;
+        end
+    end
+
+    function rsrc = setwhitenoise(rsrc, ab)
+      funcarg = strcat('BufferSize',ab);
+      if ~isfield(rsrc,funcarg)
+          rsrc.(funcarg) = 0;
+      end
+      funcarg = strcat('Buffer',ab);
+      if ~isfield(rsrc,funcarg)
+          bfsize = rsrc.(strcat('BufferSize',ab));
+          rsrc.(funcarg) = zeros(1,bfsize);
+      end
+    end
+
+    function rsrc = setparams(rsrc,mode,ab)
+        amplarg=strcat('Amplitude',ab);
+        if isfield(rsrc,amplarg)
+            switch mode
+            case m.SINE
+                rsrc = setsine(rsrc,ab);
+            case m.ARBITRARY
+                rsrc = setarb(rsrc,ab);
+                if ~isfield(rsrc,'Periodic')
+                    rsrc.Periodic = 0;
+                end
+            case m.WHITENOISE
+                rsrc = setwhitenoise(rsrc,ab);
+            end
+        end
+    end
+
+    function rsrc = setmaxorder(rsrc)
+        if isfield(rsrc,'AmplitudeA')
+            mxa=find(abs(rsrc.AmplitudeA)>0,1,'last');
+            if isempty(mxa)
+                mxa=1;
+            end
+        else
+            mxa=0;
+        end
+        if isfield(rsrc,'AmplitudeB')
+            mxb=find(abs(rsrc.AmplitudeB)>0,1,'last');
+            if isempty(mxb)
+                mxb=1;
+            end
+        else
+            mxb=0;
+        end
+        mxab=max([mxa,mxb,1]);
+        rsrc.MaxOrder=mxab-1;
+        if isfield(rsrc,'AmplitudeA')
+            rsrc.AmplitudeA(mxa+1:mxab)=0;
+        end
+        if isfield(rsrc,'AmplitudeB')
+            rsrc.AmplitudeB(mxb+1:mxab)=0;
+        end
+    end
+end
diff --git a/pyat/at/lattice/elements/variable_elements.py b/pyat/at/lattice/elements/variable_elements.py
index d0bc28425..381b27e20 100644
--- a/pyat/at/lattice/elements/variable_elements.py
+++ b/pyat/at/lattice/elements/variable_elements.py
@@ -1,163 +1,470 @@
-"""Time-dependent Multipole"""
+"""Time-dependent Multipole."""
 
 from __future__ import annotations
 
-from datetime import datetime
 from enum import IntEnum
 
-import numpy
+import numpy as np
 
+from ..exceptions import AtError, AtWarning
 from .conversions import _array
 from .element_object import Element
-from ..exceptions import AtError
 
 
 class ACMode(IntEnum):
-    """Class to define the excitation types"""
+    """Class to define the excitation types."""
 
     SINE = 0
     WHITENOISE = 1
     ARBITRARY = 2
 
 
-class VariableMultipole(Element):
-    """Class to generate an AT variable thin multipole element"""
+def _anyarray(value: np.ndarray) -> np.ndarray:
+    # Ensure proper ordering(F) and alignment(A) for "C" access in integrators
+    return np.require(value, dtype=np.float64, requirements=["F", "A"])
 
-    _BUILD_ATTRIBUTES = Element._BUILD_ATTRIBUTES
+
+class VariableThinMultipole(Element):
+    """Turn by turn variable thin multipole."""
+
+    _BUILD_ATTRIBUTES = Element._BUILD_ATTRIBUTES + ["Mode"]
     _conversions = dict(
         Element._conversions,
-        Mode=int,
         AmplitudeA=_array,
         AmplitudeB=_array,
         FrequencyA=float,
         FrequencyB=float,
         PhaseA=float,
         PhaseB=float,
-        Seed=int,
-        NSamplesA=int,
-        NSamplesB=int,
-        FuncA=_array,
-        FuncB=_array,
+        SinAabove=float,
+        SinBabove=float,
+        NsamplesA=int,
+        NsamplesB=int,
+        KtaylorA=int,
+        KtaylorB=int,
+        FuncA=_anyarray,
+        FuncB=_anyarray,
+        BufferA=_anyarray,
+        BufferB=_anyarray,
+        BufferSizeA=int,
+        BufferSizeB=int,
         Ramps=_array,
-        Periodic=bool,
+        Periodic=int,
     )
 
     def __init__(
-        self, family_name, AmplitudeA=None, AmplitudeB=None, mode=ACMode.SINE, **kwargs
-    ):
-        # noinspection PyUnresolvedReferences,SpellCheckingInspection
-        r"""
-        Parameters:
-            family_name(str):    Element name
+        self: VariableThinMultipole,
+        family_name: str,
+        mode: int or ACMode,
+        **kwargs: dict[any, any],
+    ) -> VariableThinMultipole:
+        r"""Init VariableThinMultipole.
+
+        This Class creates a thin multipole of any order (dipole kick, quadrupole,
+        sextupole, etc.) and type (Normal or Skew) defined by AmplitudeA and/or
+        AmplitudeB components; the polynoms PolynomA and PolynomB are calculated
+        on every turn depending on the chosen mode, and for some modes on the
+        particle time delay. All modes could be ramped.
+        Default pass method: ``VariableThinMPolePass``.
+
+        Keep in mind that as this element varies on every turn, and at the end of
+        the tracking PolynomA and PolynomB are set to zero, this could interfere
+        with optics calculations. Therefore, it is preferred to set the pass method
+        to IdentityPass when the element should not be active.
+
+        There are three different modes that could be set,
+        :py:attr:`.ACMode.SINE`: sine function
+        :py:attr:`.ACMode.WHITENOISE`: gaussian white noise
+        :py:attr:`.ACMode.ARBITRARY`: user defined turn-by-turn kick list
+        SINE = 0, WHITENOISE = 1 and ARBITRARY = 2. See ACMode.
+
+        For example, use at.ACMode.SINE or mode = 0 to create an sin function element.
+
+        The **SINE** mode requires amplitude and frequency for A and/or B.
+        The value of the jth component of the polynom (A or B) at the nth turn
+        is given by
+
+        Amplitude[j]*sin[TWOPI*frequency*(n*T0 + \tau_p) + phase],
+
+        where T0 is the revolution period of the ideal ring, and \tau_p is the delay
+        of the pth particle i.e. the sixth coordinate over the speed of light. Also,
+        note that the position of the element on the ring has no effect, the phase
+        could be used to add any delay due to the position along s. The following is
+        an example of the SINE mode of an skew quad
+
+        eleskew = at.VariableThinMultipole('VAR_SKEW',at.ACMode.SINE,
+        AmplitudeA=[0,skewa2],FrequencyA=freqA,PhaseA=phaseA)
+
+        The values of the sine function could be limited to be above a defined
+        threshold using ``Sin[AB]above``. For example, you could create a positive
+        half-sine by setting ``Sin[AB]above`` to zero. You could also create a
+        negative half-sine function by setting the amplitude to -1.
+
+        The **WHITENOISE** mode requires the amplitude of A and/or B. The gaussian
+        distribution is generated with zero-mean and one standard deviation from
+        a pseudo-random stream pcg32. The pcg32 seed is fixed by the tracking
+        function, therefore using the same stream on all trackings (sequencial or
+        parallel). See https://github.com/atcollab/at/discussions/879 for more
+        details on the pseudo random stream. For example,
+
+        elenoise = at.VariableThinMultipole('MYNOISE',at.ACMode.WHITENOISE,
+        AmplitudeA=[noiseA1])
+
+        creates a vertical kick as gaussian noise of amplitude noiseA1.
+
+        Additionally, one could use the parameter BufferSizeA or BufferSizeB
+        to create buffer to store the random numbers used in the element settings
+        during tracking. The buffer could be smaller or larger than the tracking.
+        If the buffer is smaller, it will store from start to the moment of filling
+        completely the buffer. If the buffer is larger, all non used places are set
+        to zero.
+
+        Care should be taken when changing the buffer and buffer size. The only check
+        done is at the moment of initialization, for example, when reading a lattice
+        from a file. Buffer and buffer size could be replaced at any moment, just
+        ensure that the actual buffer length and buffersize agree.
+
+        The **ARBITRARY** mode requires the definition of a custom discrete function
+        to be sampled at every turn. The function and its Taylor expansion with
+        respect to \tau up to any given order is
+
+        value = f(turn) + f'(turn)*tau + 0.5*f''(turn)*tau**2
+        + 1/6*f'''(turn)*tau**3 + 1/24*f''''(turn)*tau**4 ...
+
+        f is a row array of values, f',f'',f''',f'''', are row arrays containing
+        the derivatives wrt \tau, and \tau is the time delay of the particle, i.e.
+        the the sixth coordinate divided by the speed of light. Therefore, the
+        function func is a 2D-array with a given column contains f and its
+        derivatives, and the ith column has the values to be used in the ith turn.
+        For example, a positive value on the first turn f(1)=1 with positive
+        derivative f'(1)=0.1 followed by a negative value in the second turn
+        f(2)=-1 with negative derivative f'(2)=-0.2, and f(3) = 2, would be
+          Func=[[1,  -1   2],
+                [0.1 -0.2 0]].
+
+        The time tau could be offset using ``FuncATimeDelay`` or ``FuncBTimeDelay``.
+
+        tau = tau - Func[AB]TimeDelay
+
+        The function value is then **multiplied by Amplitude A and/or B**.
+        Use FuncA or FuncB, and AmplitudeA or AmplitudeB accordingly.
+        For example, the following is a positive vertical kick in the first turn,
+        negative on the second turn, and zero on the third turn.
+
+        elesinglekick = at.VariableThinMultipole('CUSTOMFUNC',at.ACMODE.ARBITRARY,
+        AmplitudeA=1e-4,FuncA=[1,-1,0],Periodic=True)
+
+        As already mentioned, one could include the derivatives of the function
+        from a Taylor expansion in a 2D array. First row for the values of f on
+        every turn, second row for the values of f' on every turn, and so on.
+
+        By default the array is assumed non periodic, the function has no effect
+        on the particle in turns exceeding the function definition. If
+        ``Periodic`` is set to 1, the sequence is repeated.
+
+        One could use the method inspect_polynom_values to check the polynom values
+        used in every turn.
+
+        Any mode could be ramped. The ramp is defined by four values (t0,t1,t2,t3):
+            * ``t<=t0``: excitation amplitude is zero.
+            * ``t0<t<=t1``: excitation amplitude is linearly ramped up.
+            * ``t1<t<=t2``: excitation amplitude is constant.
+            * ``t2<t<=t3``: excitation amplitude is linearly ramped down.
+            * ``t3<t``: excitation amplitude is zero.
+
+        Arguments:
+            family_name(str):  Element name
+            mode(at.ACMode): defines the mode. Default ACMode.SINE
+            kwargs: described as follow.
 
         Keyword Arguments:
             AmplitudeA(list,float): Amplitude of the excitation for PolynomA.
-              Default None
+            Default None
             AmplitudeB(list,float): Amplitude of the excitation for PolynomB.
-              Default None
-            mode(ACMode): defines the evaluation grid. Default ACMode.SINE
-
-              * :py:attr:`.ACMode.SINE`: sine function
-              * :py:attr:`.ACMode.WHITENOISE`: gaussian white noise
-              * :py:attr:`.GridMode.ARBITRARY`: user defined turn-by-turn kick list
-            FrequencyA(float): Frequency of the sine excitation for PolynomA
-            FrequencyB(float): Frequency of the sine excitation for PolynomB
-            PhaseA(float): Phase of the sine excitation for PolynomA. Default 0
-            PhaseB(float): Phase of the sine excitation for PolynomB. Default 0
-            MaxOrder(int): Order of the multipole for scalar amplitude. Default 0
-            Seed(int): Seed of the random number generator for white
-                       noise excitation. Default datetime.now()
-            FuncA(list): User defined tbt kick list for PolynomA
-            FuncB(list): User defined tbt kick list for PolynomB
-            Periodic(bool): If True (default) the user defined kick is repeated
-            Ramps(list): Vector (t0, t1, t2, t3) in turn number to define the ramping
-                         of the excitation
-
-              * ``t<t0``: excitation amplitude is zero
-              * ``t0<t<t1``: exciation amplitude is linearly ramped up
-              * ``t1<t<t2``: exciation amplitude is constant
-              * ``t2<t<t3``: exciation amplitude is linearly ramped down
-              * ``t3<t``: exciation amplitude is zero
+            Default None
+            FrequencyA(float): Frequency of the PolynomA sine excitation. Unit Hz
+            FrequencyB(float): Frequency of the PolynomB sine excitation. Unit Hz
+            PhaseA(float): Phase of the sine excitation for PolynomA. Default 0 rad
+            PhaseB(float): Phase of the sine excitation for PolynomB. Default 0 rad
+            SinAabove(float): Limit the sin function to be above. Default -1.1
+            SinBabove(float): Limit the sin function to be above. Default -1.1
+            BufferSizeA(int): Sets the buffersize in WHITENOISE mode.
+            BufferSizeB(int): Sets the buffersize in WHITENOISE mode.
+            FuncA(Sequence[float]): User defined turn by turn list for PolynomA
+            FuncB(Sequence[float]): User defined turn by turn list for PolynomB
+            FuncATimeDelay(float): generate a time offset on the function FUNCA.
+            It only has an effect if any of the derivatives is not zero.
+            FuncBTimeDelay(float): generate a time offset on the function FUNCB.
+            It only has an effect if any of the derivatives is not zero.
+            Periodic(bool): If True the user definition is repeated. Default False.
+            Ramps(list): Four values (t0,t1,t2,t3) defining the ramping steps:
 
-        Examples:
+        Raises:
+        :raise AtError: if none of ``AmplitudeA`` or ``AmplitudeB`` is passed.
+        :raise AtError: if ramp is not vector of length 4 when using ``Ramps``.
+        :raise AtError: when frequency is not defined if using Mode ``ACMode.SINE``.
+        :raise AtError: when the custom function is not defined if using
+        mode ``ACMode.ARBITRARY``.
 
-            >>> acmpole = at.VariableMultipole(
-            ...     "ACMPOLE", AmplitudeB=amp, FrequencyB=frequency
-            ... )
-            >>> acmpole = at.VariableMultipole(
-            ...     "ACMPOLE", AmplitudeB=amp, mode=at.ACMode.WHITENOISE
-            ... )
-            >>> acmpole = at.VariableMultipole(
-            ...     "ACMPOLE", AmplitudeB=amp, FuncB=fun, mode=at.ACMode.ARBITRARY
-            ... )
-
-        .. note::
-
-            * For all excitation modes at least one amplitude has to be provided.
-              The default excitation is ``ACMode.SINE``
-            * For ``mode=ACMode.SINE`` the ``Frequency(A,B)`` corresponding to the
-              ``Amplitude(A,B)`` has to be provided
-            * For ``mode=ACMode.ARBITRARY`` the ``Func(A,B)`` corresponding to the
-              ``Amplitude(A,B)`` has to be provided
+        Examples:
+            >>> acmpole = at.VariableThinMultipole('ACSKEW',
+            at.ACMode.SINE, AmplitudeA=[0,amp], FrequencyA=freq)
+            >>> acmpole = at.VariableThinMultipole('ACHKICK',
+            at.ACMode.WHITENOISE, AmplitudeB=amp)
+            >>> acmpole = at.VariableThinMultipole('ACMPOLE',
+            at.ACMode.ARBITRARY, AmplitudeB=[0,0,amp], FuncB=fun)
         """
-        kwargs.setdefault("PassMethod", "VariableThinMPolePass")
-        self.MaxOrder = kwargs.pop("MaxOrder", 0)
-        self.Periodic = kwargs.pop("Periodic", True)
-        self.Mode = int(mode)
-        if AmplitudeA is None and AmplitudeB is None:
-            raise AtError("Please provide at least one amplitude for A or B")
-        AmplitudeB = self._set_params(AmplitudeB, mode, "B", **kwargs)
-        AmplitudeA = self._set_params(AmplitudeA, mode, "A", **kwargs)
-        self._setmaxorder(AmplitudeA, AmplitudeB)
-        if mode == ACMode.WHITENOISE:
-            self.Seed = kwargs.pop("Seed", datetime.now().timestamp())
-        self.PolynomA = numpy.zeros(self.MaxOrder + 1)
-        self.PolynomB = numpy.zeros(self.MaxOrder + 1)
-        ramps = kwargs.pop("Ramps", None)
-        if ramps is not None:
-            assert len(ramps) == 4, "Ramps has to be a vector with 4 elements"
-            self.Ramps = ramps
-        super().__init__(family_name, **kwargs)
 
-    def _setmaxorder(self, ampa, ampb):
-        mxa, mxb = 0, 0
-        if ampa is not None:
-            mxa = numpy.max(numpy.nonzero(ampa))
-        if ampb is not None:
-            mxb = numpy.max(numpy.nonzero(ampb))
-        self.MaxOrder = max(mxa, mxb)
-        if ampa is not None:
-            delta = self.MaxOrder - len(ampa)
-            if delta > 0:
-                ampa = numpy.pad(ampa, (0, delta))
-            self.AmplitudeA = ampa
-        if ampb is not None:
-            delta = self.MaxOrder + 1 - len(ampb)
-            if delta > 0:
-                ampb = numpy.pad(ampb, (0, delta))
-            self.AmplitudeB = ampb
-
-    def _set_params(self, amplitude, mode, ab, **kwargs):
-        if amplitude is not None:
-            if numpy.isscalar(amplitude):
-                amp = numpy.zeros(self.MaxOrder)
-                amplitude = numpy.append(amp, amplitude)
+        def _check_amplitudes(**kwargs: dict[any, any]) -> dict[str, any]:
+            amp_aux = {"A": None, "B": None}
+            all_amplitudes_are_none = True
+            for key in amp_aux:
+                if "amplitude" + key in kwargs:
+                    raise AtWarning(
+                        "amplitude" + key + " should be Amplitude" + key + "."
+                    )
+                lower_case_kwargs = {k.lower(): v for k, v in kwargs.items()}
+                amp_aux[key] = lower_case_kwargs.pop("amplitude" + key.lower(), None)
+                if amp_aux[key] is not None:
+                    all_amplitudes_are_none = False
+            if all_amplitudes_are_none:
+                raise AtError("Please provide at least one amplitude for A or B")
+            return amp_aux
+
+        def _getmaxorder(ampa: np.ndarray, ampb: np.ndarray) -> int:
+            mxa, mxb = 0, 0
+            if ampa is not None:
+                mxa = np.max(np.append(np.nonzero(ampa), 0))
+            if ampb is not None:
+                mxb = np.max(np.append(np.nonzero(ampb), 0))
+            return max(mxa, mxb)
+
+        def _set_thismode(
+            mode: int, a_b: str, **kwargs: dict[any, any]
+        ) -> dict[str, any]:
             if mode == ACMode.SINE:
-                self._set_sine(ab, **kwargs)
+                kwargs = _set_sine(a_b, **kwargs)
             if mode == ACMode.ARBITRARY:
-                self._set_arb(ab, **kwargs)
-        return amplitude
-
-    def _set_sine(self, ab, **kwargs):
-        frequency = kwargs.pop("Frequency" + ab, None)
-        phase = kwargs.pop("Phase" + ab, 0)
-        assert frequency is not None, "Please provide a value for Frequency" + ab
-        setattr(self, "Frequency" + ab, frequency)
-        setattr(self, "Phase" + ab, phase)
-
-    def _set_arb(self, ab, **kwargs):
-        func = kwargs.pop("Func" + ab, None)
-        nsamp = len(func)
-        assert func is not None, "Please provide a value for Func" + ab
-        setattr(self, "Func" + ab, func)
-        setattr(self, "NSamples" + ab, nsamp)
+                kwargs = _set_arb(a_b, **kwargs)
+            if mode == ACMode.WHITENOISE:
+                kwargs = _set_white_noise(a_b, **kwargs)
+            return kwargs
+
+        def _set_sine(a_b: str, **kwargs: dict[any, any]) -> dict[str, any]:
+            frequency = kwargs.get("Frequency" + a_b, None)
+            if frequency is None:
+                raise AtError("Please provide a value for Frequency" + a_b)
+            kwargs.setdefault("Phase" + a_b, 0)
+            kwargs.setdefault("Sin" + a_b + "above", -1.1)
+            return kwargs
+
+        def _set_arb(a_b: str, **kwargs: dict[any, any]) -> dict[str, any]:
+            func = kwargs.get("Func" + a_b, None)
+            if func is None:
+                raise AtError("Please provide a value for Func" + a_b)
+            if len(np.squeeze(func.shape)) == 1:
+                nsamples = len(func)
+                ktaylor = 1
+            else:
+                ktaylor, nsamples = func.shape
+            kwargs.setdefault("Func" + a_b + "TimeDelay", 0)
+            kwargs["NSamples" + a_b] = nsamples
+            kwargs["Ktaylor" + a_b] = ktaylor
+            kwargs.setdefault("Periodic", False)
+            return kwargs
+
+        def _set_white_noise(a_b: str, **kwargs: dict[any, any]) -> dict[any, any]:
+            bfsize = int(kwargs.setdefault("BufferSize" + a_b, 0))
+            kwargs.setdefault("Buffer" + a_b, np.zeros((bfsize)))
+            kwargs["BufferSize" + a_b] = bfsize
+            # we check the case when Buffer and BufferSize don't match
+            the_userbuffershape = kwargs["Buffer" + a_b].shape
+            if the_userbuffershape[0] != kwargs["BufferSize" + a_b]:
+                raise AtError(("Buffer and BufferSizeMismatch"))
+            return kwargs
+
+        def _check_ramp(**kwargs: dict[any, any]) -> _array:
+            ramps = kwargs.get("Ramps", None)
+            if ramps is not None:
+                if len(ramps) != 4:
+                    raise AtError("Ramps has to be a vector with 4 elements")
+                ramps = np.asarray(ramps)
+            return ramps
+
+        # init begins
+        kwargs.setdefault("Mode", mode)
+        kwargs.setdefault("PassMethod", "VariableThinMPolePass")
+        if len(kwargs) > 2:
+            amp_aux = _check_amplitudes(**kwargs)
+            for key, value in amp_aux.items():
+                if value is not None:
+                    kwargs["Amplitude" + key] = value
+                    kwargs = _set_thismode(mode, key, **kwargs)
+            maxorder = _getmaxorder(amp_aux["A"], amp_aux["B"])
+            kwargs["MaxOrder"] = kwargs.get("MaxOrder", maxorder)
+            for key in amp_aux:
+                kwargs["Polynom" + key] = kwargs.get(
+                    "Polynom" + key, np.zeros(maxorder + 1)
+                )
+            ramps = _check_ramp(**kwargs)
+            if ramps is not None:
+                kwargs["Ramps"] = ramps
+        super().__init__(family_name, **kwargs)
+
+    def inspect_polynom_values(
+        self: VariableThinMultipole, **kwargs: dict[any, any]
+    ) -> dict[str, list]:
+        """
+        Get the polynom values per turn.
+
+        Translations (T1,T2) and Rotations (R1,R2) in the element are ignored.
+
+        Arguments:
+            kwargs: as follow.
+
+        Keyword Arguments:
+            turns(int): Default 1. Number of turns to calculate.
+            T0(float): revolution time in seconds. Use only in SINE mode.
+            tparticle(float): Default 0. Time of the particle in seconds.
+            kwargs: as needed by the mode.
+
+        Returns:
+            Dictionary with a list of PolynomA and PolynomB per turn.
+        """
+        # initialize in case of early return
+        listpola = []
+        listpolb = []
+
+        turns = kwargs.setdefault("turns", 1)
+        mode = self.Mode
+        timeoffset = 0
+        if mode == 0:
+            # revolution time
+            trevol = float(kwargs["T0"])
+            tparticle = float(kwargs.setdefault("tparticle", 0))
+            timeoffset = trevol + tparticle
+        elif mode == 3:
+            print("Random values are not available here.")
+        elif mode == 2:
+            # particle time
+            timeoffset = float(kwargs.setdefault("tparticle", 0))
+        ramps = getattr(self, "Ramps", 0)
+        periodic = getattr(self, "Periodic", False)
+        maxorder = self.MaxOrder
+
+        pola = np.full(maxorder + 1, np.nan)
+        polb = np.full(maxorder + 1, np.nan)
+
+        for turn in range(turns):
+            for order in range(maxorder + 1):
+                if hasattr(self, "AmplitudeA"):
+                    pola[order] = self._get_pol(
+                        "A", ramps, mode, timeoffset * turn, turn, order, periodic
+                    )
+                if hasattr(self, "AmplitudeB"):
+                    polb[order] = self._get_pol(
+                        "B", ramps, mode, timeoffset * turn, turn, order, periodic
+                    )
+            listpola.append(np.copy(pola))
+            listpolb.append(np.copy(polb))
+        return {"PolynomA": listpola, "PolynomB": listpolb}
+
+    def _get_amp(
+        self: VariableThinMultipole, amp: float, ramps: _array, _time: float
+    ) -> float:
+        """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.
+
+        Parameters:
+            amp: amplitude component.
+            ramps: array containing the turns that define the ramp
+            _time: turn
+
+        Returns:
+            amp if no ramp.
+            amp multiplied by the ramp state.
+        """
+        ampt = amp
+        if ramps != 0:
+            if _time <= ramps[0]:
+                ampt = 0.0
+            elif _time <= ramps[1]:
+                ampt = amp * (_time - ramps[0]) / (ramps[1] - ramps[0])
+            elif _time <= ramps[2]:
+                ampt = amp
+            elif _time <= ramps[3]:
+                ampt = amp - amp * (_time - ramps[2]) / (ramps[3] - ramps[2])
+            else:
+                ampt = 0.0
+        return ampt
+
+    def _get_pol(
+        self: VariableThinMultipole,
+        a_b: str,
+        ramps: _array,
+        mode: int,
+        _time: float,
+        turn: int,
+        order: int,
+        periodic: bool,
+    ) -> float:
+        """
+        Return the polynom component of a given order.
+
+        Parameters:
+            a_b: either 'A' or 'B' indicating the polynom.
+            ramps: array containing the ramp definition.
+            mode: value to specify the type of variable element.
+            _time: time for this mode
+            turn: turn to check
+            order: order of the polynom
+            periodic: whether the sequence is periodic or not.
+
+        Returns:
+            the amplitude for the polynom component
+        """
+        allamp = getattr(self, "Amplitude" + a_b)
+        amp = allamp[order]
+        ampout = 0
+        # check if amp is zero
+        if amp == 0:
+            return ampout
+
+        # get the ramp value
+        ampout = self._get_amp(amp, ramps, turn)
+
+        if mode == 0:
+            # sin mode parameters
+            whole_sin_above = getattr(self, "Sin" + a_b + "above")
+            freq = getattr(self, "Frequency" + a_b)
+            phase = getattr(self, "Phase" + a_b)
+            sinval = np.sin(2 * np.pi * freq * _time + phase)
+            if sinval >= whole_sin_above:
+                ampout = ampout * sinval
+            else:
+                ampout = 0
+        elif mode == 1:
+            ampout = np.nan
+        elif mode == 2:
+            nsamples = getattr(self, "NSamples" + a_b)
+            if periodic or turn < nsamples:
+                func = getattr(self, "Func" + a_b)
+                functdelay = float(getattr(self, "Func" + a_b + "TimeDelay"))
+                turnidx = np.mod(turn, nsamples)
+                ktaylor = int(getattr(self, "Ktaylor" + a_b))
+                _time = _time - functdelay
+                functot = func[0, turnidx]
+                thefactorial = 1
+                tpow = 1
+                for i in range(1, ktaylor):
+                    tpow = _time * tpow
+                    thefactorial = thefactorial * i
+                    functot = functot + tpow / thefactorial * func[i, turnidx]
+                ampout = ampout * functot
+            else:
+                ampout = 0.0
+        else:
+            ampout = 0.0
+        return ampout
diff --git a/pyat/at/load/utils.py b/pyat/at/load/utils.py
index e6ccb9347..45803aa6a 100644
--- a/pyat/at/load/utils.py
+++ b/pyat/at/load/utils.py
@@ -29,6 +29,8 @@
 from at.lattice import AtWarning
 from at.lattice import elements as elt
 from at.lattice import Lattice, Particle, Element, Marker
+from at.lattice import idtable_element
+from at.lattice import variable_elements
 
 _ext_suffix = sysconfig.get_config_var("EXT_SUFFIX")
 _plh = "placeholder"
@@ -124,6 +126,7 @@ def __init__(
     "AperturePass": elt.Aperture,
     "IdTablePass": elt.InsertionDeviceKickMap,
     "GWigSymplecticPass": elt.Wiggler,
+    "VariableThinMPolePass": variable_elements.VariableThinMultipole,
 }
 
 # Lattice attributes which must be dropped when writing a file
diff --git a/pyat/at/tracking/utils.py b/pyat/at/tracking/utils.py
index 9f854ee5d..dde80f0c7 100644
--- a/pyat/at/tracking/utils.py
+++ b/pyat/at/tracking/utils.py
@@ -8,7 +8,7 @@
 from typing import Optional
 from ..lattice import Lattice, Element
 from ..lattice import BeamMoments, Collective, QuantumDiffusion
-from ..lattice import SimpleQuantDiff, VariableMultipole
+from ..lattice import SimpleQuantDiff, VariableThinMultipole
 from ..lattice import elements, refpts_iterator, set_value_refpts
 from ..lattice import DConstant, checktype, checkattr, get_bool_index
 
@@ -20,7 +20,7 @@
 
 DIMENSION_ERROR = 'Input to lattice_pass() must be a 6xN array.'
 _COLLECTIVE_ELEMS = (BeamMoments, Collective)
-_VAR_ELEMS = (QuantumDiffusion, SimpleQuantDiff, VariableMultipole)
+_VAR_ELEMS = (QuantumDiffusion, SimpleQuantDiff, VariableThinMultipole)
 _DISABLE_ELEMS = _COLLECTIVE_ELEMS + _VAR_ELEMS