NotebookEvaluate["/home/max/projects/mathematica/style.nb"]

(*component parameters*)
(* eγ, μ, eγμ, eγπ *)
ca={0,1.,1.,1.1};
cγ={0.982,0.963,0.949,0.947};
cn19={1.16*10^10,7.21*10^8,1.04*10^8,8.15*10^8};
crg={410,660,725,174};
cs={0.4,1.1,1.14,1.33};
cdx1={-600,-400,-400,-500};
cdxm0={140,0,154,45};
cdxm1={0,0.426,0.141,0};
cλ0={58,280,278,98};
cλ1={0.014,0.201,0.162,0.011};

dx[r_,ψ_,θ_,Xmax_]:=880/Cos[θ] Exp[-(r Sin[θ]Cos[ψ])/8600]-Xmax
mgh[dx_, dxref_, dxmax_, dx1_,λ_]:=((dx-dx1)/(dxref-dx1))^((dxmax-dx1)/λ)Exp[-(dx-dxref)/λ]
nkg[r_,rg_,s_,n19_]:=n19/(2π rg^2)*Gamma[9/2-s]/(Gamma[s]Gamma[9/2-2s])(r/rg)^(s-2)(1+r/rg)^(s-9/2)

(*time model is not expressed here for all components, just total, also shape σ is assumed const. *)
(*x1 = 50g/cm^2 *)
RX[dx_,x_,xmax_,θ_]:=8600./Cos[θ] Log[(dx+xmax)/x]
tq[dx_,xq_,xmax_,θ_,r_]:=Sqrt[RX[dx, xq, xmax,θ]^2+r^2]-RX[dx, xq, xmax,θ]
σ=0.6;
n[t_,dx_,xmax_,θ_,r_]:=1/(Sqrt[2π]σ (t-tq[dx,50,xmax,θ,r]))Exp[-1/(2σ^2)(Log[(t-tq[dx,50,xmax,θ,r])/(tq[dx,xmax(*40pct quantile*),xmax,θ,r]-tq[dx,50,xmax,θ,r])]+Sqrt[2]σ InverseErf[-0.2])^2]

rho[lgE0_, γ_,dx_,dxref_,dxmax_,dx1_,λ_,r_,rg_,s_,n19_]:=10^(γ*(lgE0-19)) * mgh[dx, dxref, dxmax, dx1,λ]* nkg[r,rg,s,n19]
rhocomp[lgE0_,Rmu_,dx_,r_,c_]:=(ca[[c]](Rmu-1)+1)rho[lgE0, cγ[[c]],dx,200,cdxm0[[c]]+cdxm1[[c]]*r,cdx1[[c]],cλ0[[c]]+cλ1[[c]],r,crg[[c]],cs[[c]],cn19[[c]]]
rhotot[lgE0_,Rmu_,dx_,r_]:=Sum[
(ca[[c]](Rmu-1)+1)
 rho[lgE0, cγ[[c]],dx,200,cdxm0[[c]]+cdxm1[[c]]*r,cdx1[[c]],cλ0[[c]]+cλ1[[c]],r,crg[[c]],cs[[c]],cn19[[c]]]
,{c,{1,2,3,4}}]

Manipulate[

Plot[
{rhocomp[lgE0,Rmu,dx,r,1],
 rhocomp[lgE0,Rmu,dx,r,2],
 rhocomp[lgE0,Rmu,dx,r,3],
 rhocomp[lgE0,Rmu,dx,r,4],

rhotot[lgE0,Rmu,dx,r]},
{dx,-200,1200},PlotStyle->{Red,Blue,Green,Orange,Black},Evaluate@customStyle,
FrameLabel->{T["ΔX / g (cm^(-2))"],T["ρ"]},PlotLegends->{T["eγ"],T["μ"],T["eγμ"],T["eγπ"],T["total"]}]

,{{lgE0,19},18,20},{{Rmu,1},0.5,2.5},{{r,500},200,2000}]

Manipulate[

LogPlot[{
rhocomp[lgE0,Rmu,dx[r,ψ,θ,Xmax],r,1],
 rhocomp[lgE0,Rmu,dx[r,ψ,θ,Xmax],r,2],
 rhocomp[lgE0,Rmu,dx[r,ψ,θ,Xmax],r,3],
 rhocomp[lgE0,Rmu,dx[r,ψ,θ,Xmax],r,4],
rhotot[lgE0,Rmu,dx[r,ψ,θ,Xmax],r]}
,{r,200,1800},PlotStyle->{Red,Blue,Green,Orange,Black},Evaluate@customStyle,PlotRange->{{0,2000},{1,Automatic}},FrameLabel->{T["r / m"],T["ρ"]},PlotLegends->{T["eγ"],T["μ"],T["eγμ"],T["eγπ"],T["total"]}]

,{{θ,12Degree},0,50Degree},{{lgE0,19},18,20},{{Xmax,700},650,900},{{Rmu,1},0.5,2.5},{{ψ,π/2},0,π}]

Manipulate[

LogPlot[
{rhotot[lgE0,Rmu,dx[r,ψ,θ Degree,{Xmax-dxmax,Xmax+dxmax}],r],
rhotot[lgE0,{Rmu-drmu,Rmu+drmu},dx[r,ψ,θ Degree,Xmax],r]
},{r,200,1800},PlotStyle->{Orange,Black},Evaluate@customStyle,PlotRange->{{0,2000},{1,Automatic}},FrameLabel->{T["r / m"],T["ρ"]},PlotLegends->{T["δ(Subscript[X, max])"],T["(Subscript[δR, μ])"]}]

,{{θ,12},0,50},{{lgE0,19},18,20},{{Xmax,700},650,900},{{dxmax,20},0,70},{{Rmu,1},0.5,2.5},{{drmu,0.1},0,0.3},{{ψ,π/2},0,π}]

Manipulate[

Plot[
{rhotot[lgE0,Rmu,dx[r,ψ,θ Degree,{Xmax-dxmax,Xmax+dxmax}],r]/rhotot[lgE0,Rmu,dx[r,ψ,θ Degree,Xmax],r],
rhotot[lgE0,{Rmu-drmu,Rmu+drmu},dx[r,ψ,θ Degree,Xmax],r]/rhotot[lgE0,Rmu,dx[r,ψ,θ Degree,Xmax],r]
},{r,200,1800},PlotStyle->{Orange,Black},Evaluate@customStyle,PlotRange->{{0,2000},{0.5,1.5}},FrameLabel->{T["r / m"],T["δρ/ρ"]},PlotLegends->{T["δ(Subscript[X, max])"],T["δ(Subscript[R, μ])"]}]

,{{θ,12},0,50},{{lgE0,19},18,20},{{Xmax,700},650,900},{{dxmax,20},0,60},{{Rmu,1},0.5,2.5},{{drmu,0.1},0,0.3},{{ψ,π/2},0,π}]

Manipulate[

Plot[n[t,dx[r,ψ,θ Degree,Xmax],Xmax,θ Degree ,r],{t,1,3000},
PlotStyle->{Orange},Evaluate@customStyle,PlotRange->{{0,3000},{0,0.005}},FrameLabel->{T["t / ns"],T["n"]}
]
,{{θ,12},0,50},{{r,1000},300,2000},{{Xmax,700},650,900},{{ψ,π/2},0,π}]

rgr[x_,y_,θ_]:=Sqrt[x^2+(y/Cos[θ])^2]
tgr[t_,x_,y_,θ_]:=t+rgr[x,y,θ]Tan[θ]y/Sqrt[x^2+y^2]

animation=Evaluate[With[{θ=45,lgE0=19,Xmax=700,Rmu=1},
showergif=Animate[
Quiet[Plot3D[Log10[rhocomp[lgE0,Rmu,dx[rgr[x,y,θ Degree],ArcCos[y/Sqrt[x^2+y^2]],θ Degree,Xmax],rgr[x,y,θ Degree],1]+10]n[tgr[t,x,y,θ  Degree],dx[rgr[x,y,θ],ArcCos[y/Sqrt[x^2+y^2]],θ Degree,Xmax],Xmax,θ Degree ,rgr[x,y,θ Degree]],
{y,-4000,4000},{x,-4000,4000},
ImageSize->500,PlotRange->{{-5000,5000},{-5000,5000},{0,0.005}},
PerformanceGoal->"Quality",WorkingPrecision->MachinePrecision,RegionFunction->Function[{x,y},Sqrt[x^2+y^2]>100]]
],{t,-1000,4000 },AnimationRate->75,RefreshRate->2,AnimationDirection->Forward]
]];

(*animation*)