Energy ramp

atintegrators/BndMPoleSymplectic4QuantPass.c

@@ -253,7 +253,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
             Elem->fringeIntM0, Elem->fringeIntP0,
             Elem->T1, Elem->T2, Elem->R1, Elem->R2,
             Elem->RApertures, Elem->EApertures,
-            Elem->KickAngle, Elem->Scaling, Elem->Energy,
+            Elem->KickAngle, Elem->Scaling, Param->energy,
             Param->thread_rng, num_particles);
     return Elem;
 }

atintegrators/BndMPoleSymplectic4RadPass.c

@@ -59,6 +59,8 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do
     double B0 = B[0];
     double A0 = A[0];
 
+    /*printf("In a dipole. E0 = %f eV",E0);*/
+
     if (KickAngle) {   /* Convert corrector component to polynomial coefficients */
         B[0] -= sin(KickAngle[0])/le;
         A[0] += sin(KickAngle[1])/le;
@@ -143,7 +145,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
         BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error();
         EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error();
         ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error();
-        Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error();
+        /*Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error();*/
+        Energy=Param->energy;
         /*optional fields*/
         FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error();
         FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error();
@@ -201,7 +204,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
             Elem->fringeIntM0, Elem->fringeIntP0,
             Elem->T1, Elem->T2, Elem->R1, Elem->R2,
             Elem->RApertures, Elem->EApertures,
-            Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles);
+            Elem->KickAngle, Elem->Scaling, Param->energy, num_particles);
     return Elem;
 }
 

atintegrators/EnergyRampPass.c

@@ -0,0 +1,163 @@
+/* 
+ *  EnergyRampPass.c
+ *  Accelerator Toolbox 
+ *  25/07/2024
+ *  Nicola Carmignani
+ */
+
+#include "atconstants.h"
+#include "atelem.c"
+
+
+struct elem
+{
+    double E0;
+    int NPointsRamp;
+    double* TurnsRamp;
+    double* EnergyRamp;
+    /* internal variable */
+    int CurrentRampIndex;
+};
+
+void EnergyRampPass(double *r_in, struct elem *Elem, int nturn, 
+        double *Energy, int num_particles)
+{
+    double E_entrance=*Energy;
+    double E_exit;
+    int t1, t2, endramp=0;
+    double E1, E2;
+    if (Elem->CurrentRampIndex>=Elem->NPointsRamp)
+        endramp=1;
+    if (nturn >= Elem->TurnsRamp[Elem->CurrentRampIndex] && !endramp)
+        Elem->CurrentRampIndex++;
+    if (endramp)
+        E_exit = Elem->EnergyRamp[Elem->NPointsRamp-1];
+    else
+    {
+        /* linear interpolation */
+        if (Elem->CurrentRampIndex==0)
+        {
+            t1=0;
+            E1 = Elem->E0;
+        }
+        else
+        {
+            t1 = Elem->TurnsRamp[Elem->CurrentRampIndex-1];
+            E1 = Elem->EnergyRamp[Elem->CurrentRampIndex-1];
+        }
+        t2 = Elem->TurnsRamp[Elem->CurrentRampIndex];
+        E2 = Elem->EnergyRamp[Elem->CurrentRampIndex];
+        E_exit = E1 + ((E2-E1)/(t2-t1))*(nturn-t1);
+    }
+    *Energy=E_exit;
+
+    /* loop in the particles to update delta coordinate */
+    for (int c = 0; c<num_particles; c++) { /*Loop over particles  */
+        /* Get the pointer to the current particle's state vector. */
+        double *r6 = r_in + 6*c;
+        if (!atIsNaN(r6[0])) {
+            double p_norm;
+            double xpr, ypr;
+
+            /* calculate angles from tranverse momenta 	*/
+            p_norm = 1.0 + r6[4];
+            xpr = r6[1]/p_norm;
+            ypr = r6[3]/p_norm;
+
+            /* change of delta coordinate */
+            r6[4] = (E_entrance * (1.0 + r6[4]) - E_exit) / E_exit;
+
+            /* recalculate momenta from angles after change of delta */
+            p_norm = 1.0 + r6[4];
+            r6[1] = xpr*p_norm;
+            r6[3] = ypr*p_norm;
+        }
+    }
+}
+
+#if defined(MATLAB_MEX_FILE) || defined(PYAT)
+ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
+			      double *r_in, int num_particles, struct parameters *Param)
+{
+    double Energy_entrance=Param->energy;
+    double Energy_exit, E0;
+    int nturn=Param->nturn;
+
+    if (!Elem) {
+        int NPointsRamp;
+        double *TurnsRamp, *EnergyRamp;
+        int CurrentRampIndex=0;
+
+        E0 = atGetDouble(ElemData,"E0"); check_error();
+        EnergyRamp = atGetDoubleArray(ElemData,"EnergyRamp"); check_error();
+        NPointsRamp = atGetLong(ElemData,"NPointsRamp"); check_error();
+        TurnsRamp = atGetDoubleArray(ElemData,"TurnsRamp"); check_error();
+
+        Elem = (struct elem*)atMalloc(sizeof(struct elem));
+
+        Elem->NPointsRamp=NPointsRamp;
+        Elem->TurnsRamp=TurnsRamp;
+        Elem->EnergyRamp=EnergyRamp;
+        Elem->CurrentRampIndex=CurrentRampIndex;
+    }
+    EnergyRampPass(r_in, Elem, nturn, &Param->energy,
+            num_particles);
+    return Elem;
+}
+
+MODULE_DEF(EnergyRampPass)        /* Dummy module initialisation */
+
+#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/
+
+#if defined(MATLAB_MEX_FILE)
+
+void mexFunction(	int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{ 	
+  if(nrhs == 2)
+    {
+      double *r_in;
+      const mxArray *ElemData = prhs[0];
+      int num_particles = mxGetN(prhs[1]);
+      double *EnergyRamp, *TurnsRamp, E0;
+      int NPointsRamp;
+      struct elem *Elem;
+      EnergyRamp = atGetDoubleArray(ElemData,"EnergyRamp"); check_error();
+      NPointsRamp = atGetLong(ElemData,"NPointsRamp"); check_error();
+      E0 = atGetDouble(ElemData,"E0"); check_error();
+      TurnsRamp = atGetDoubleArray(ElemData,"TurnsRamp"); check_error();  
+      Elem = (struct elem*)atMalloc(sizeof(struct elem));
+
+      Elem->NPointsRamp=NPointsRamp;
+      Elem->TurnsRamp=TurnsRamp;
+      Elem->EnergyRamp=EnergyRamp;
+      Elem->CurrentRampIndex=0;
+      Elem->E0=E0;
+
+      int nturn=0;
+      double Energy=200e6;
+
+      if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix");
+      /* ALLOCATE memory for the output array of the same size as the input  */
+      plhs[0] = mxDuplicateArray(prhs[1]);
+      r_in = mxGetDoubles(plhs[0]);
+      EnergyRampPass(r_in, Elem, nturn, &Energy, num_particles);
+    }
+  else if (nrhs == 0)
+  {   /* return list of required fields */
+      plhs[0] = mxCreateCellMatrix(4,1);
+      mxSetCell(plhs[0],0,mxCreateString("E0"));
+      mxSetCell(plhs[0],1,mxCreateString("EnergyRamp"));
+      mxSetCell(plhs[0],2,mxCreateString("NPointsRamp"));
+      mxSetCell(plhs[0],3,mxCreateString("TurnsRamp"));
+      if(nlhs>1) /* optional fields */
+      {
+          plhs[1] = mxCreateCellMatrix(0,1);
+      }
+  }
+  else
+  {
+      mexErrMsgIdAndTxt("AT:WrongArg","Needs 0 or 2 arguments");
+  }
+
+}
+#endif /* MATLAB_MEX_FILE */

atintegrators/RFCavityPass.c

@@ -18,16 +18,55 @@ struct elem
   double HarmNumber;
   double TimeLag;
   double PhaseLag;
+  /* variables for voltage ramp */
+  int NPointsRamp;
+  double* TurnsRamp;
+  double* VoltageRamp;
+  /* internal variable */
+  int CurrentRampIndex;
 };
 
-void RFCavityPass(double *r_in, double le, double nv, double freq, double h, double lag, double philag,
-                  int nturn, double T0, int num_particles)
+void RFCavityPass(double *r_in, double le, double Voltage, double energy, 
+        double freq, double h, double lag, double philag, int nturn, 
+        double T0, int NPointsRamp, double *TurnsRamp, double *VoltageRamp, 
+        int* CurrentRampIndex, int num_particles)
 /* le - physical length
    nv - peak voltage (V) normalized to the design enegy (eV)
    r is a 6-by-N matrix of initial conditions reshaped into
    1-d array of 6*N elements
 */
 {
+    double nv;
+    if (NPointsRamp)
+    {
+        int  endramp=0;
+        double t1, t2, V1, V2;
+        if (*CurrentRampIndex>=NPointsRamp)
+            endramp=1;
+        if (nturn >= TurnsRamp[*CurrentRampIndex] && !endramp)
+            *CurrentRampIndex=*CurrentRampIndex+1;
+        if (endramp)
+            Voltage = VoltageRamp[NPointsRamp-1];
+        else
+        {
+            /* linear interpolation */
+            if (*CurrentRampIndex==0)
+            {
+                t1=0;
+                V1 = Voltage;
+            }
+            else
+            {
+                t1 = TurnsRamp[*CurrentRampIndex-1];
+                V1 = VoltageRamp[*CurrentRampIndex-1];
+            }
+            t2 = TurnsRamp[*CurrentRampIndex];
+            V2 = VoltageRamp[*CurrentRampIndex];
+            Voltage = V1 + ((V2-V1)/(t2-t1))*(nturn-t1);
+        }
+    }
+
+    nv=Voltage/energy;
     trackRFCavity(r_in, le, nv, freq, h, lag, philag, nturn, T0, num_particles);
 }
 
@@ -38,24 +77,35 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
     int nturn=Param->nturn;
     double T0=Param->T0;
     if (!Elem) {
-        double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag;
+        double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag, *VoltageRamp, *TurnsRamp;
+        int NPointsRamp;
         Length=atGetDouble(ElemData,"Length"); check_error();
         Voltage=atGetDouble(ElemData,"Voltage"); check_error();
         Energy=atGetDouble(ElemData,"Energy"); check_error();
         Frequency=atGetDouble(ElemData,"Frequency"); check_error();
         TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error();
         PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error();
+        VoltageRamp = atGetOptionalDoubleArray(ElemData,"VoltageRamp"); check_error();
+        NPointsRamp = atGetOptionalLong(ElemData,"NPointsRamp",0); check_error();
+        TurnsRamp = atGetOptionalDoubleArray(ElemData,"TurnsRamp"); check_error();
         Elem = (struct elem*)atMalloc(sizeof(struct elem));
         Elem->Length=Length;
         Elem->Voltage=Voltage;
-        Elem->Energy=Energy;
+        Elem->Energy=Param->energy;
         Elem->Frequency=Frequency;
         Elem->HarmNumber=round(Frequency*T0);
         Elem->TimeLag=TimeLag;
         Elem->PhaseLag=PhaseLag;
+        Elem->NPointsRamp=NPointsRamp;
+        Elem->TurnsRamp=TurnsRamp;
+        Elem->VoltageRamp=VoltageRamp;
+        Elem->CurrentRampIndex=0;
     }
-    RFCavityPass(r_in, Elem->Length, Elem->Voltage/Elem->Energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag,
-                 Elem->PhaseLag, nturn, T0, num_particles);
+
+    RFCavityPass(r_in, Elem->Length, Elem->Voltage, Param->energy, Elem->Frequency, 
+            Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, 
+            Elem->NPointsRamp, Elem->TurnsRamp, Elem->VoltageRamp, 
+            &Elem->CurrentRampIndex,num_particles);
     return Elem;
 }
 
@@ -80,12 +130,16 @@ void mexFunction(	int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       double PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0);
       double T0=1.0/Frequency;      /* Does not matter since nturns == 0 */
       double HarmNumber=round(Frequency*T0);
+      double *VoltageRamp = atGetOptionalDoubleArray(ElemData,"VoltageRamp");
+      int NPointsRamp = atGetOptionalLong(ElemData,"NPointsRamp",0);
+      double *TurnsRamp = atGetOptionalDoubleArray(ElemData,"TurnsRamp");
       if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix");
       /* ALLOCATE memory for the output array of the same size as the input  */
       plhs[0] = mxDuplicateArray(prhs[1]);
       r_in = mxGetDoubles(plhs[0]);
-      RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles);
-
+      RFCavityPass(r_in, Length, Voltage, Energy, Frequency, HarmNumber,
+              TimeLag, PhaseLag, 0, T0, NPointsRamp, TurnsRamp, VoltageRamp,
+              0, num_particles);
     }
   else if (nrhs == 0)
   {   /* return list of required fields */
@@ -96,9 +150,12 @@ void mexFunction(	int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
       mxSetCell(plhs[0],3,mxCreateString("Frequency"));
       if(nlhs>1) /* optional fields */
       {
-          plhs[1] = mxCreateCellMatrix(2,1);
+          plhs[1] = mxCreateCellMatrix(5,1);
           mxSetCell(plhs[1],0,mxCreateString("TimeLag"));
           mxSetCell(plhs[1],1,mxCreateString("PhaseLag"));
+          mxSetCell(plhs[1],2,mxCreateString("NPointsRamp"));
+          mxSetCell(plhs[1],3,mxCreateString("TurnsRamp"));
+          mxSetCell(plhs[1],4,mxCreateString("VoltageRamp"));
       }
   }
   else

atintegrators/StrMPoleSymplectic4QuantPass.c

@@ -211,7 +211,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
             Elem->T1, Elem->T2, Elem->R1, Elem->R2,
             Elem->RApertures, Elem->EApertures,
             Elem->KickAngle, Elem->Scaling,
-            Elem->Energy, Param->thread_rng, num_particles);
+            Param->energy, Param->thread_rng, num_particles);
     return Elem;
 }
 

atintegrators/StrMPoleSymplectic4RadPass.c

@@ -118,7 +118,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
         PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error();
         MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error();
         NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error();
-        Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error();
+        /*Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error();*/
+        Energy=Param->energy;
         /*optional fields*/
         Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error();
         FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error();

atmat/lattice/element_creation/atEnergyRamp.m

@@ -0,0 +1,30 @@
+function elem=atEnergyRamp(fname,varargin)
+%ATENERGYRAMP Creates an energy ramp element with Class 'EnergyRamp'
+%ATENERGYRAMP(FAMNAME,TurnsRamp,EnergyRamp)
+%
+%  INPUTS
+%  1. FAMNAME	   - Family name
+%  2. TurnsRamp    - array (1,NPointsRamp) with turn numbers with change of
+%                    slope for theenergy ramp
+%  3. EnergyRamp   - Energy at the turn TurnsRamp
+%
+%  OUTPUTS
+%  1. ELEM - Structure with the AT element
+%
+%  EXAMPLES
+%  1. atEnergyRamp('EnRamp',[0,100,1000],[200e6,220e6,1000e6]);
+%
+
+
+[rsrc,turnsramp, energyramp] = decodeatargs({0,200e6},varargin);
+[method,rsrc]     = getoption(rsrc,'PassMethod','EnergyRampPass');
+[cl,rsrc]         = getoption(rsrc,'Class','EnergyRamp');
+
+NPointsRamp=min([length(turnsramp),length(energyramp)]);
+
+% Build the element
+elem=atbaselem(fname,method,'Class',cl,'Length',0,...
+    'TurnsRamp',turnsramp,'EnergyRamp',energyramp,...
+    'NPointsRamp',NPointsRamp,rsrc{:});
+
+end