plotopts = {Frame->True, BaseStyle->{FontSize->15}, ImageSize->600};
SetOptions[{ParametricPlot,Plot,ListPlot},plotopts];

Import

file = "csr_wake.dat";
rdat = Import[NotebookDirectory[]<>file, "Table"]

smatch = "!#-----------------------------";

plist = Position[rdat, smatch][[All,1]];
nwake = Length[plist]

wakei[i_] := rdat[[plist[[i]]+4;;plist[[i]]+4 + 38]];
{ListPlot[{#[[1]], #[[3]]}&/@wakei[nwake]], ListPlot[{#[[1]], #[[2]]}&/@wakei[nwake]]}

ListPlot[{#[[1]], #[[5]]}&/@wakei[nwake]]

wmax = Max[Flatten[Table[Abs[wakei[i][[All,3]]], {i,nwake}]]]

densitymax = Max[Flatten[Table[Abs[wakei[i][[All,2]]], {i,nwake}]]]

si[i_] := rdat[[plist[[i]]+1]][[3]];
si[nwake]

Waverage[wakei_]:= Module[{dtot = Total[wakei[[All,2]]]},Total[(1/dtot) #[[2]] #[[3]]&/@wakei]];
Waverage[wakei[nwake]]

Wsigma[wakei_]:= Module[{dtot = Total[wakei[[All,2]]]},Sqrt[Total[(1/dtot) #[[2]] ((#[[3]])^2)&/@wakei]-((Waverage[wakei])^2)]];
Wsigma[wakei[nwake]]

scale = (10^6);
waverages = Table[{si[i],scale Waverage[wakei[i]]}, {i,nwake}];
wsigmas = Table[{si[i],scale Wsigma[wakei[i]]}, {i,nwake}];

waveragesDict[file] = waverages;
wsigmasDict[file] = wsigmas ;
ListPlot[{waverages, wsigmas}, Joined->True, PlotRange->All, FrameLabel->{"s (m)", "Average, sigma wake"}]

cumulativeloss =Total[waverages[[All,2]]]

Manipulate[ListPlot[{#[[1]], #[[3]]}&/@wakei[i], FrameLabel->{"z", "wake", "S = "<>ToString[si[i]]}, PlotRange->3 (10^(-7)){-1,1}], {i,1,nwake,1}]

Comparison

ListPlot[{waveragesDict["bcsr_wake.dat"], wsigmasDict["bcsr_wake.dat"],waveragesDict["qcsr_wake.dat"], wsigmasDict["qcsr_wake.dat"]}, Joined->True, PlotRange->{{0,.6},All}, PlotStyle->{Red,Blue,{Dashed,Red}, {Dashed, Blue}}]

Energy Loss

plotopts = {Frame->True, BaseStyle->{FontSize->15}, ImageSize->600};
SetOptions[{Plot,ListPlot}, plotopts];

c = 299792458.
rc = 2.817 (10^(-15));
σ = 0.000899377;
hchamber = 24 (10^(-3));

echarge = 1.6 (10^(-19));

qbunch = 1 (10^(-9));
np = qbunch/echarge;
energy = 1000 (10^6);
mc2 = .511 (10^6);
γ = energy/mc2;
β = (Sqrt[1-1/(γ^2)]);

rawedat = {#[[2]], (10^(-3 ))#[[3]]}&/@Import[NotebookDirectory[]<>"beampz.dat", "Table"];
edat ={ #[[1]], energy#[[2]]}&/@rawedat;
ListPlot[edat]

κf[bfield_]:= bfield (c/energy);

bdat = {#[[2]],#[[3]]}&/@Import[NotebookDirectory[]<>"b_field.dat", "Table"];
{ListPlot[bdat], ListPlot[{#[[1]], κf[#[[2]]]}&/@bdat]}

ipad= {{150,50}, {50,10}};
gout = Column[
{ ListPlot[{#[[1]], κf[#[[2]]]}&/@bdat, AspectRatio->1/4, FrameLabel->{"s (m)", "κ (1/m)"}, PlotStyle->Blue, Joined->True, ImagePadding->ipad],
ListPlot[bdat, AspectRatio->1/4, FrameLabel->{"", "(Subscript[B, y]) (T)"}, PlotStyle->Blue, Joined->True, ImagePadding->ipad]}]

Export[NotebookDirectory[]<>"curvature.pdf", gout]

ListPlot[bdat, AspectRatio->1/4, FrameLabel->{"s (m)", "(Subscript[B, y]) (T)"}, PlotStyle->Blue, Joined->True]

Δ = Mean[Differences[bdat[[All,1]]]];

wdat = {#[[1]], np rc mc2 (((Abs[κf[#[[2]]]])^(2/3))/(σ^(4/3)))}&/@bdat;
ListPlot[wdat]

shieldingfactordat = {#[[1]], hchamber  (((Abs[κf[#[[2]]]])^(1/3))/(σ^(2/3)))}&/@bdat;
ListPlot[shieldingfactordat]

val=0; 
westimatedat = Table[val = val -.35 Δ d[[2]]; {d[[1]],val}, {d, wdat}];
ListPlot[westimatedat]

gout2 = 
 Column[{
ListPlot[{Legended[edat, "Beam"], Legended[westimatedat, "Analytic"]}, FrameLabel->{"s (m)", "ΔE (eV)"}, Joined->True, ImagePadding->ipad, AspectRatio->1/3],
ListPlot[bdat, AspectRatio->1/4, FrameLabel->{"", "(Subscript[B, y]) (T)"}, PlotStyle->Blue, Joined->True, ImagePadding->ipad]}]

Export[NotebookDirectory[]<>"csr_loss.pdf", gout2]

Column[

Export[NotebookDirectory[]

model = a x;
pars = {a};
var = x;
a/.FindFit[edat, model, pars, var]

a/.FindFit[westimatedat, model, pars, var]