diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index b542a70f9..bce3be029 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/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; } diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index acb47b660..eb65ffca7 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/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; } diff --git a/atintegrators/EnergyRampPass.c b/atintegrators/EnergyRampPass.c new file mode 100644 index 000000000..02500d429 --- /dev/null +++ b/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; cenergy; + 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 */ diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index c2ac702f6..322259c82 100755 --- a/atintegrators/RFCavityPass.c +++ b/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 diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index 83fa35448..a885ca320 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/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; } diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 88fb00470..b9ab3775b 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/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(); diff --git a/atmat/lattice/element_creation/atEnergyRamp.m b/atmat/lattice/element_creation/atEnergyRamp.m new file mode 100644 index 000000000..6dd6a2212 --- /dev/null +++ b/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