Bmad to OPAL test

root=NotebookDirectory[];
scratchdir="~/Scratch/";

(* Set plot defaults *)
myPlotOptions={PlotRange->All,Frame->True,BaseStyle->{FontSize->15},ImageSize->{500}};
SetOptions[ListPlot,myPlotOptions];
SetOptions[Graphics,myPlotOptions];
SetOptions[ParametricPlot,myPlotOptions];
myPlotOptions={PlotRange->All,Frame->True,BaseStyle->{FontSize->15},ImageSize->{500}};
SetOptions[ListLogPlot,myPlotOptions];
SetOptions[Histogram,myPlotOptions];
SetOptions[ListDensityPlot,myPlotOptions];

(*Commonly used constants in SI units*)
cLightN=299792458.;
mec2N=0.510998928 (10^6);

Parse OPAL stat file utility

ParseColumnD[onecolumnD_]:=Module[{split},
(*parse a SDDS line, with some assumtions about the positions of things*)
split=StringSplit[onecolumnD, Characters[",=&"]][[1]];
{split[[2]], split[[6]], split[[8]]}
]
ImportOPALStat[statfile_String]:=Module[{begP, endP, nparam,opalstat, opalstatcolumns, opalstatparams},
opalstat=Check[Import[statfile],{ Print["Problem importing "<>statfile], Abort[]}];

(*Look for these header markers*)
begP=Position[opalstat,"SDDS1"][[1,1]];
endP=Position[opalstat,"&data mode=ascii &end"][[1,1]];

(*set global variables for data description, "D" stands for "description"*)
paramD=Select[opalstat[[begP;;endP]], StringMatchQ[#[[1]],"&parameter*"]&];
columnD=Select[opalstat[[begP;;endP]], StringMatchQ[#[[1]],"&column*"]&];

nparam=Length[paramD];


{statnames, statunits,statdescription} = Transpose[Table[ParseColumnD[columnD[[i]]], {i,1,Length[columnD]}]];
Do[OPALix[statnames[[i]]]=i, {i,1,Length[statnames]}];
Do[OPALunits[statnames[[i]]]= statunits[[i]], {i,1,Length[statnames]}];


(*set global variables for data*)
 opalstatparams=opalstat[[endP+1;;endP+1+nparam]];
opalstatcolumns=opalstat[[endP+1+nparam+1;;]];
Print["Parsed "<>statfile];
opalstatcolumns
];

OPALStatLabels[]:= Grid[{{"Column", "Name", "Unit", "Description"}}~Join~Table[Flatten[{i, ParseColumnD[columnD[[i]]]}], {i,1,Length[columnD]}], Frame->True, Alignment->"."];




Clear[PlotOPALStat];
PlotOPALStat[c1_, c2_, opalstat_, opts:OptionsPattern[]]:=Module[{dat, name1,name2, units1,units2,description1,description2},
dat={#[[c1]], #[[c2]]}&/@opalstat;
{name1, units1, description1}=ParseColumnD[columnD[[c1]]];
{name2, units2, description2}=ParseColumnD[columnD[[c2]]];


ListPlot[dat,Evaluate[FilterRules[{opts}, Options[ListPlot]]],FrameLabel->{name1<>" ("<>units1<>")",name2<>" ("<>units2<>")"}, Joined->True,PlotRange->All, PlotStyle->Black]
];

Clear[OPALStatArray];
(*Makes phase space array from slice*)
OPALStatArray[statVar_String, opalstat_]:=Module[{f},
OPALGetScaleAndLabel[statVar]; (*just used to check*)
f = Switch[statVar,
"energy",(mec2N + (10^6)#[[OPALix["energy"]]])&,
"kinetic_energy",( (10^6)#[[OPALix["energy"]]])&,

_, #[[OPALix[statVar]]]&
];
f/@ opalstat

];


Clear[OPALStatPlot];
OPALStatPlot[yVar_String, xVar_String,opalstat_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist, datx, daty},
{xscale,xlabel}=OPALGetScaleAndLabel[xVar];
{yscale,ylabel}=OPALGetScaleAndLabel[yVar];

datx=OPALStatArray[xVar, opalstat];
daty=OPALStatArray[yVar, opalstat];

ListPlot[Transpose[{xscale datx, yscale daty}],
FrameLabel->{xlabel, ylabel},Evaluate[FilterRules[{opts}, Options[ListPlot]]], Joined->True]


];
Options[OPALStatPlot] =  {PlotRange->All,  Joined->True };

Dark Current Tracker Analysis

Clear[BmadGetScaleAndLabel];
BmadGetScaleAndLabel[PSvariable_String]:=Module[{},
Switch[PSvariable,
	"x", {(10^3), "x (mm)"},
	"cp_x",{ 1, "c (Subscript[p, x]) (eV)"},
	"y",{(10^3), "y (mm)"},
	"cp_y",{1, "c (Subscript[p, y]) (eV)"},
	"s",{1, "s (m)"},
	"z",{1, "z (m)"},
	"cp_s", { (1/(10^6)), "c (Subscript[p, s]) (MeV)"},
	"cp_z", { (1/(10^6)), "c (Subscript[p, z]) (MeV)"},
	"t",   {(10^9), "t (ns)"},
	"E_x", { (10^(-6)), "(Subscript[E, x]) (MV/m)"},
	"E_y", { (10^(-6)), "(Subscript[E, y]) (MV/m)"},
	"E_s", { (10^(-6)), "(Subscript[E, s]) (MV/m)"},
	"E_z", { (10^(-6)), "(Subscript[E, z]) (MV/m)"},
	"B_x", { (10^3), "(Subscript[B, x]) (mT)"},
	"B_y", { (10^3), "(Subscript[B, y]) (mT)"},
	"B_s", { (10^3), "(Subscript[B, s]) (mT)"},
	"B_z", { (10^3), "(Subscript[B, z]) (mT)"},


	_, Print["Bad variable: "<>PSvariable];Abort[]
]
];

ImportBmadTimeParticles[particlefile_String]:=Module[{dat,ppos, timetracks, cnames, cunits},

Print["Importing Bmad time particles: "<>particlefile];
Check[dat=Import[particlefile, "Table"];, Print["Problem with: "<>particlefile];Abort[];];
cnames=dat[[1]]; (*names of columns*)
cunits=dat[[2]];(*Units of columns*)

Do[ix[cnames[[i]]]=i, {i,1,Length[cnames]}];
Do[units[cnames[[i]]]=cunits[[i]], {i,1,Length[cnames]}];

ppos=Position[dat, "particle"][[All,1]];
timetracks = Table[dat[[ppos[[i]]+1;;ppos[[i+1]]-1 ]], {i,1,Length[ppos]-1}];
timetracks  =timetracks~Join~{dat[[Last[ppos]+1;; ]]}; (*Append final particle*)
Print["Imported "<>ToString[Length[ppos]]<>" particles, "<>ToString[Length[dat]]<>" data points"];
timetracks
]

Clear[BmadPlot];
BmadPlot[xVar_String,yVar_String,timeTracks_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist,tracks, dat, fy},
{xscale,xlabel}=BmadGetScaleAndLabel[xVar];
{yscale,ylabel}=BmadGetScaleAndLabel[yVar];

tracks=timeTracks;

dat = Table[{xscale  #[[ix[xVar]]],yscale  #[[ix[yVar]]] }&/@tracks[[i]], {i,1,Length[tracks]}];

ListPlot[dat,
FrameLabel->{xlabel, ylabel},Evaluate[FilterRules[{opts}, Options[ListPlot]]], Joined->True]

];
Clear[BmadPlot3D];
BmadPlot3D[xVar_String,yVar_String,zVar_String, timeTracks_, opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist, zscale, zlabel, tracks, dat, fy},
{xscale,xlabel}=BmadGetScaleAndLabel[xVar];
{yscale,ylabel}=BmadGetScaleAndLabel[yVar];
{zscale,zlabel}=BmadGetScaleAndLabel[zVar];

tracks=timeTracks;

dat = Table[{xscale  #[[ix[xVar]]],yscale  #[[ix[yVar]]],zscale  #[[ix[zVar]]] }&/@tracks[[i]], {i,1,Length[tracks]}];

ListPointPlot3D[dat,
AxesLabel->{xlabel, ylabel, zlabel},Evaluate[FilterRules[{opts}, Options[ListPointPlot3D]]]]

];

Options[BmadPlot] =  {PlotRange->All,  Joined->True }

Options[BmadPlot]

OPAL Analysis

OPALH5Trajectories[h5file_, property_String,  step_Integer]:=Module[{idlist, xlist,stepstring, info, ix,ref},
(*Sort particles by their ID*)
stepstring= "/Step#"<>ToString[step];
Check[idlist=Import[h5file, {"Datasets",stepstring<>"/id"}];, Abort[]];
info= Import[h5file,{"Attributes",stepstring}];
(*ref =Riffle[info["RefPartR"],info["RefPartP"] ];*)
ref =Riffle[info["RefPartR"],{0,0,0} ];
ix["x"] = 1;ix["px"] = 2;ix["y"] = 3;ix["py"] = 4;ix["z"] = 5;ix["pz"] = 6;

xlist=Import[h5file, {"Datasets", "/Step#"<>ToString[step]<>"/"<>property<>""}] + ref[[ix[property]]];
Sort[Transpose[{idlist, xlist}]][[All,2]]
]


OPALPhaseSpaceNames = {"x", "px", "y", "py", "z", "pz"};

Clear[OPALGetScaleAndLabel];
OPALGetScaleAndLabel[PSvariable_String]:=Module[{},
Switch[PSvariable,
	"x", {(10^3), "x (mm)"},
	"px",{ mec2N, "c (Subscript[p, x]) (eV)"},
	"y",{(10^3), "y (mm)"},
	"py",{ mec2N, "c (Subscript[p, y]) (eV)"},
	"z",{1, "z (m)"},
	"pz", { (1/(10^6))mec2N, "c (Subscript[p, z]) (MeV)"},
	"Ex", { (10^(-6)), "(Subscript[E, x]) (MV/m)"},
	"Ey", { (10^(-6)), "(Subscript[E, y]) (MV/m)"},
	"Ez", { (10^(-6)), "(Subscript[E, z]) (MV/m)"},
	"t",   {(10^9), "t (ns)"},
	"Bx", { (10^3), "(Subscript[B, x]) (mT)"},
	"By", { (10^3), "(Subscript[B, y]) (mT)"},
	"Bz", { (10^3), "(Subscript[B, z]) (mT)"},

(*For stat *)
	"s",   {1, "s (m)"},
	"charge", {-(10^12), "-charge (pC)"},
	"numParticles",   {1, "(Subscript[N, particles])"},
	"mean_x",   {(10^6), "<x> (μm)"},
	"mean_y",   {(10^6), "<y> (μm)"},
	"mean_s",   {1, "<s> (m)"},
	"rms_x",   {(10^6), "(Subscript[σ, x]) (μm)"},
	"rms_y",   {(10^6), "(Subscript[σ, y]) (μm)"},
	"rms_s",   {(10^6), "(Subscript[σ, z]) (μm)"},
	"rms_px",   {mec2N, "c*(Subscript[σ, px]) (eV)"},
	"rms_py",   {mec2N, "c*(Subscript[σ, py]) (eV)"},
	"rms_ps",   {mec2N, "c*(Subscript[σ, pz]) (eV)"},
	"emit_x", {(10^6), "\!\(\*SubscriptBox[\(ϵ\), \(N, \\ x\)]\) (mm-mrad)"},
	"emit_y", {(10^6), "\!\(\*SubscriptBox[\(ϵ\), \(N, \\ y\)]\) (mm-mrad)"},
	"emit_s", {(10^6), "emit_s (?)"},
	"Dx", {1, "(Subscript[D, x]) (m)"}, (*Dispersion, need to test*)
	"By_ref", {1, "(Subscript[B, y]) reference (T)"},
	"energy",   {(10^(-6)), "Energy (MeV)"},
	"kinetic_energy",   {(10^(-6)), "Kinetic Energy (MeV)"},
	"dE", {(10^6), "Energy Spread (eV)"},
	_, Print["Bad variable: "<>PSvariable];Abort[]
]
];


OPALStepList[h5file_]:=Module[{datasetlist, steplist},

Check[datasetlist=Import[h5file, {"Datasets"}];, Abort[]];
steplist=ToExpression[(StringSplit[#, {"/Step#", "/"}]&/@Select[datasetlist,StringMatchQ [#,"*Step*"]&])[[All,1]]];
Intersection[steplist,steplist]

];

Clear[OPALTracks];
OPALTracks[xVar_, h5file_]:=Module[{xlist, steplist, datasetlist, nparticles,tlist, max},
Check[datasetlist=Import[h5file, {"Datasets"}];, Abort[]];
steplist=ToExpression[(StringSplit[#, {"/Step#", "/"}]&/@Select[datasetlist,StringMatchQ [#,"*Step*"]&])[[All,1]]];
maxstep=Max[steplist];
If[xVar=="t",
(*Make time list out of the steps*) 
nparticles=Length[Import[h5file, {"Datasets", "/Step#0/id"}]];
tlist = OPALstepsize *Table[i-1, {i,1, maxstep+1}];
xlist=Transpose[Table[tlist, {nparticles}]];,
(*else use this function*)
xlist = Table[OPALH5Trajectories[h5file, xVar, i], {i,0,maxstep}];
];
Transpose[xlist]
];


Clear[OPALPlot];
OPALPlot[yVar_String, xVar_String, opalh5file_,  whichparticles_:All,opts:OptionsPattern[]]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist},
{xscale,xlabel}=OPALGetScaleAndLabel[xVar];
{yscale,ylabel}=OPALGetScaleAndLabel[yVar];

xlist = xscale*OPALTracks[xVar, opalh5file];
ylist = yscale*OPALTracks[yVar, opalh5file];


ListPlot[Table[Transpose[{xlist[[i]], ylist[[i]]}], {i,1,Length[xlist]}][[whichparticles]],Evaluate[FilterRules[{opts}, Options[ListPlot]]],
 Joined->True, FrameLabel->{xlabel, ylabel}, PlotStyle->Red, PlotRange->All]

];

Clear[OPALPointPlot3D];
OPALPointPlot3D[xVar_String, yVar_String,zVar_String, opalh5file_,  whichparticles_:All]:=Module[{yscale, ylabel, xscale,xlabel,xlist,  ylist, zscale, zlist, zlabel},
{xscale,xlabel}=OPALGetScaleAndLabel[xVar];
{yscale,ylabel}=OPALGetScaleAndLabel[yVar];
{zscale,zlabel}=OPALGetScaleAndLabel[zVar];

xlist = xscale*OPALTracks[xVar, opalh5file];
ylist = yscale*OPALTracks[yVar, opalh5file];
zlist = zscale*OPALTracks[zVar, opalh5file];


ListPointPlot3D[Table[Transpose[{xlist[[i]], ylist[[i]],  zlist[[i]]}], {i,1,Length[xlist]}][[whichparticles]], PlotRange->All,AxesLabel->{xlabel,ylabel,zlabel}, BaseStyle->{FontSize->15} ]

];

Clear[OPALPhaseSpace];
OPALPhaseSpace[xVar_,yVar_, opalh5file_, step_:0, opts:OptionsPattern[]]:= Module[{xlist, ylist, yscale, ylabel, xscale,xlabel},
{xscale,xlabel}=OPALGetScaleAndLabel[xVar];
{yscale,ylabel}=OPALGetScaleAndLabel[yVar];


xlist =xscale OPALH5Trajectories[ opalh5file,xVar, step];
ylist =yscale  OPALH5Trajectories[ opalh5file,yVar, step];

ListPlot[Transpose[{xlist, ylist}],Evaluate[FilterRules[{opts}, Options[ListPlot]]],   Joined->False, FrameLabel->{xlabel, ylabel, "Step: "<>ToString[step]}, PlotStyle->Black, PlotRange->All]

];

(* Single phase space variable at a single step *)
OPALPhaseSpaceData[xVar_, opalh5file_, step_:0]:= Module[{xlist, xscale,xlabel},
{xscale,xlabel}=OPALGetScaleAndLabel[xVar];
xscale OPALH5Trajectories[ opalh5file,xVar, step]

];

OPALSigmaMatrixTerm[i_Integer,j_Integer,opalh5file_, step_:0]:=Module[{getvar, charge,wC,averagei,averagej,nparticles, iVar, jVar, ilist, jlist, pairs},
getvar[ii_Integer]:=Switch[ii, 
1, "x",
2, "px",
3, "y",
4, "py",
5, "z",
6, "pz",
_, Print["Bad index: "<>ToString[i]];Abort
];


iVar = getvar[i];
jVar = getvar[j];
ilist = OPALH5Trajectories[ opalh5file,iVar, step];
jlist = OPALH5Trajectories[ opalh5file,jVar, step];
pairs = Transpose[{ilist, jlist}];
averagei=Mean[ilist ];
averagej=Mean[jlist ];
(*Subtract off Centroid for Covariance*)
Mean[(#[[1]]-averagei) (#[[2]]-averagej)&/@pairs]
];
OPALSigmaMatrix[opalh5file_, step_:0]:=Module[{coords},
coords = Transpose[Table[ OPALH5Trajectories[ opalh5file,var, step], {var, {"x", "px", "y", "py", "z", "pz"}}]];
Covariance[coords,coords]
];

Line test

workdir = root<>"new_sbend/"

Bmad

bmadTracks=ImportBmadTimeParticles[workdir<>"lat.out"];

BmadPSList= {"x", "y", "z","cp_x", "cp_y", "cp_z"};
BmadEBList={"E_x", "E_y", "E_z", "B_x", "B_y", "B_z"};

Clear[BmadPlotAll];
BmadPlotAll[ bmadTracks_, xVar_String, yVarList_, opts:OptionsPattern[]]:=Module[{},
Table[BmadPlot[xVar, i, bmadTracks,opts], {i, yVarList}]
]

BmadPlot["z","cp_z", bmadTracks]

BmadPlot["z","x", bmadTracks]

BmadPSList

xVar="z";
gBmad=BmadPlotAll[bmadTracks, xVar, Join[BmadPSList]]

bmadEnd = Import[workdir<>"lat.end_opal_particles", "Table"][[2;;]];
ListPlot[{#[[5]], #[[6]]}&/@bmadEnd, PlotRange->All]

OPAL

OPALstepsize =100(* psdumpfreq *)1.0 (10^(-12))(*DT*); 

file="opal.h5";
h5 =workdir<>file

datasets=Import[h5];
maxstep=Max[OPALStepList[h5]]

params = Import[h5, {"Attributes", "/"}];
Grid[Table[{k, params[k]}, {k, Keys[params]}], Alignment->Left, Frame->True]

params = Import[h5, {"Attributes", "/Step#"<>ToString[maxstep]}];
Grid[Table[{k, params[k]}, {k, Keys[params]}], Alignment->Left, Frame->True]

OPALPhaseSpace["z","x",h5, 5]

OPALH5Trajectories[h5, "pz",maxstep]

OPALPlot["y","z", h5,All]

OPALPSList= {"x", "y", "z","px", "py", "pz"};
OPALEBList={"Ex", "Ey", "Ez", "Bx", "By", "Bz"};

Clear[OPALPlotAll];
OPALPlotAll[ opalh5_, xVar_String, yVarList_, opts:OptionsPattern[]]:=Module[{},
Table[OPALPlot[i,xVar, opalh5,All,opts], {i, yVarList}]
]

xVar="z";
gOPAL=OPALPlotAll[h5, xVar, Join[OPALPSList], PlotStyle->Red];

Comparison

Show/@Transpose[{ gOPAL,gBmad}]//Column

Sbend

workdir = root<>"sbend/"

Bend fringe fields

1DProfile1 6[entrance field order]  7[exit field order] 3.0[gap height]
 -6.0[beginning of fringe, relative to ELEMEDGE, -2.0[Origin of ENGE function] 2.0[end of fringe] 1000
 24.0 28.0 32.0 0
0.00000e+00 
4.36222e-06
 8.83270e-06
  + 9 lines 
  1.32490e-05 
  1.73710e-05 
  2.18598e-05

Figure B.5: A 1D field map describing the fringe field of an element using 7 Enge coefficients for the entrance fringe field and 8 Enge coefficients for the exit fringe field (polynomial order 6 and 7 respectively). The element has a gap height of 3.0 cm, the entrance fringe field is non-negligible from −6.0 cm and reaches the core strength at 2.0 cm relative to ELEMEDGE. The origin of the Enge function for the entrance fringe field is at −2.0 cm relative to ELEMEDGE. The exit fringe field is non-negligible up to 32.0 cm from ELEMEDGE and starts to deviate from the core strength after 24.0 cm from ELEMEDGE. The origin of the Enge function for the exit fringe field is at 28.0 cm from ELEMEDGE.

EngeF[s_, d_, coef_List]:=(1/(1+Exp[Sum[coef[[n]] (((s/d))^(n-1)), {n,1, Length[coef]}]]))

EngeF[s, 1, {1,2}]

Integrate[(1/(1+(E^(a+b s)))), s]

clist={1,10};
Plot[EngeF[s, 1,clist], {s, -5,5},PlotRange->All, AxesOrigin->{0,0}]

clist={-10,145};
Plot[EngeF[s, 1,clist], {s, -.1,.1},PlotRange->All, AxesOrigin->{0,0}]

d

Plot[EngeF[s, 1, clist](1-EngeF[s, 1, clist]), {s, -1,1},PlotRange->All, AxesOrigin->{0,0}]

clist= {a,b};
Integrate[EngeF[s,H, clist](1-EngeF[s,H, clist]), {s,-∞,∞}]

OPAL numbers: c0 = 0.478959 c1 = 1.911289 c2 = −1.185953 c3 = 1.630554 c4 = −1.082657 c5 = 0.318111

clist = {0.478959, 1.911289, −1.185953,1.630554,−1.082657,0.318111};
Plot[EngeF[s, 0.02,clist], {s, -.1,.1},PlotRange->All, AxesOrigin->{0,0}, Frame->True]

Bmad

(* Import Data *)
trackDat=Import[workdir<>"lat.out", "Table"][[3;;]];Dimensions[trackDat]

(* Find where each data set begins and ends *)
begin=Position[trackDat, "particle"][[All,1]]+numHeaderLines;
end=begin[[#+1]]-numHeaderLines-1&/@Range[Length[begin]-1];
end=Append[end,Length[trackDat]];

(* Make an array of all data sets and one of phases *)
myTrackDat=trackDat[[begin[[#]];;end[[#]]]]&/@Range[Length[begin]];

tableTrack=Table[{#[[ix["t"]]],#[[ix["x"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];

ListPlot[tableTrack, Joined->True, PlotRange->All]

gbmadx=ListPlot[Table[ {(10^9)#[[ix["t"]]],1000*#[[ix["x"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "x (mm)"}]

gbmady=ListPlot[Table[ {(10^9)#[[ix["t"]]],1000*#[[ix["y"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "y (mm)"}]

gbmads=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["s"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "s (m)"}]

gbmadpx=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["cpx"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, x]) (eV)"}]

gbmadpy=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["cpy"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, y]) (eV)"}]

gbmadps=ListPlot[Table[ {(10^9)#[[ix["t"]]],(10^(-6))#[[ix["cps"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, s]) (MeV)"}]

gbmadenergy=ListPlot[Table[ {(10^9)#[[ix["t"]]],(10^(-6))(Sqrt[(mec2N^2)+((#[[ix["cpx"]]])^2)+((#[[ix["cpy"]]])^2)+((#[[ix["cps"]]])^2)])}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "Energy (MeV)"}, GridLines->{Automatic, Table[i,{i,1,12}]}]

OPAL

h5file = workdir<>"opal.h5";
Run["cp "<>h5file <>" ~/Scratch"];
localh5file="~/Scratch/opal.h5";

(*datasets=Import[localh5file];
maxstep=Length[datasets]/9-1*)
maxstep=10

gopalx=OPALPlot["x",maxstep]

gopaly=OPALPlot["y",maxstep]

gopalz=OPALPlot["z",maxstep]

gopalpx=OPALPlot["px",maxstep]

gopalpy=OPALPlot["py",maxstep]

gopalpz=OPALPlot["pz",maxstep]

Comparison

Show[gbmadx, gopalx]

Show[gbmady, gopaly]

Show[gbmadpx, gopalpx]

Show[gbmadpy, gopalpy, PlotRange->All]

Show[gbmads, gopalz]

Show[ gopalpz,gbmadps]

Comparison

Linac

workdir = root<>"linac/"

Stat file

ParseOPALStat[workdir<>"opal.stat"]

Show[PlotOPALStat[1,5], AxesOrigin->{0,0}]

Bend fringe fields

1DProfile1 6[entrance field order]  7[exit field order] 3.0[gap height]
 -6.0[beginning of fringe, relative to ELEMEDGE, -2.0[Origin of ENGE function] 2.0[end of fringe] 1000
 24.0 28.0 32.0 0
0.00000e+00 
4.36222e-06
 8.83270e-06
  + 9 lines 
  1.32490e-05 
  1.73710e-05 
  2.18598e-05

Figure B.5: A 1D field map describing the fringe field of an element using 7 Enge coefficients for the entrance fringe field and 8 Enge coefficients for the exit fringe field (polynomial order 6 and 7 respectively). The element has a gap height of 3.0 cm, the entrance fringe field is non-negligible from −6.0 cm and reaches the core strength at 2.0 cm relative to ELEMEDGE. The origin of the Enge function for the entrance fringe field is at −2.0 cm relative to ELEMEDGE. The exit fringe field is non-negligible up to 32.0 cm from ELEMEDGE and starts to deviate from the core strength after 24.0 cm from ELEMEDGE. The origin of the Enge function for the exit fringe field is at 28.0 cm from ELEMEDGE.

EngeF[s_, d_, coef_List]:=(1/(1+Exp[Sum[coef[[n]] (((s/d))^(n-1)), {n,1, Length[coef]}]]))

EngeF[s, 1, {1,2}]

Integrate[(1/(1+(E^(a+b s)))), s]

clist={1,10};
Plot[EngeF[s, 1,clist], {s, -5,5},PlotRange->All, AxesOrigin->{0,0}]

Plot[EngeF[s, 1, clist](1-EngeF[s, 1, clist]), {s, -1,1},PlotRange->All, AxesOrigin->{0,0}]

clist= {a,b};
Integrate[EngeF[s, H, clist](1-EngeF[s, H, clist]), {s,-∞,∞}]

Bmad

(* Import Data *)
trackDat=Import[workdir<>"ReferenceParticles.out", "Table"][[3;;]];Dimensions[trackDat]

(* Find where each data set begins and ends *)
begin=Position[trackDat, "particle"][[All,1]]+numHeaderLines;
end=begin[[#+1]]-numHeaderLines-1&/@Range[Length[begin]-1];
end=Append[end,Length[trackDat]];

(* Make an array of all data sets and one of phases *)
myTrackDat=trackDat[[begin[[#]];;end[[#]]]]&/@Range[Length[begin]];

tableTrack=Table[{#[[ix["t"]]],#[[ix["x"]]]}&/@myTrackDat[[i]],{i,1,Length[myTrackDat]}];

ListPlot[tableTrack, Joined->True, PlotRange->All]

gbmadx=ListPlot[Table[ {(10^9)#[[ix["t"]]],1000*#[[ix["x"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "x (mm)"}]

gbmady=ListPlot[Table[ {(10^9)#[[ix["t"]]],1000*#[[ix["y"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "y (mm)"}]

gbmads=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["s"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "s (m)"}]

gbmadpx=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["cpx"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, x]) (eV)"}]

gbmadpy=ListPlot[Table[ {(10^9)#[[ix["t"]]],#[[ix["cpy"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, y]) (eV)"}]

gbmadps=ListPlot[Table[ {(10^9)#[[ix["t"]]],(10^(-6))#[[ix["cps"]]]}&/@myTrackDat[[i]], {i,1,Length[myTrackDat] }], PlotRange->All, Joined->True, BaseStyle->{FontSize->15}, Frame->True, FrameLabel->{"t (ns)", "c (Subscript[p, s]) (MeV)"}]

OPAL

h5file = workdir<>"opal.h5";
Run["cp "<>h5file <>" ~/Scratch"];
localh5file="~/Scratch/opal.h5";

(*datasets=Import[localh5file];
maxstep=Length[datasets]/9-1*)
maxstep=100

gopalx=OPALPlot["x",50]

gopaly=OPALPlot["y",maxstep]

gopalz=OPALPlot["z",maxstep]

gopalpx=OPALPlot["px",maxstep]

gopalpy=OPALPlot["py",maxstep]

gopalpz=OPALPlot["pz",maxstep]

Comparison

Show[gbmadx, gopalx]

Show[gbmady, gopaly]

Show[gbmads, gopalz]

Show[gbmadpx, gopalpx]

Show[gbmadpy, gopalpy, PlotRange->All]

Show[gopalpz]

Show[gbmadps, gopalpz, PlotRange->{1.8, 4}]