Clear[intf,intfsqr,intfoverp,exppTsqr,exppzsqr]
intf[f_,lp0_,lpf_,κf_,pgoal_:20]:=(2/(((2π))^2))NIntegrate[Exp[3lp]f,{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal];
intfsqr[f_,lp0_,lpf_,κf_,pgoal_:20]:=(2/(((2π))^2))NIntegrate[Exp[3lp](f^2),{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal];
intfoverp[f_,lp0_,lpf_,κf_,pgoal_:20]:=(2/(((2π))^2))NIntegrate[Exp[2lp] f,{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal]
exppTsqr[f_,lp0_,lpf_,κf_,pgoal_:20]:=(2/(((2π))^2))NIntegrate[(1-(κ^2))Exp[5lp]f,{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal];
exppzsqr[f_,lp0_,lpf_,κf_,pgoal_:20]:=(2/(((2π))^2))NIntegrate[(κ^2)Exp[5lp]f,{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal];

Clear[lCb,mDsqr,makeLogτvals]
lCb[pmax_,pmin_]:=Log[(pmax/pmin)]
mDsqr[Ifoverp_,Nc_,g_]:=4 Nc (g^2)Ifoverp
makeLogτvals[begin_,end_,length_]:=(
((10^(#))&/@Array[#&,length,{((Log[begin])/(Log[10])),((Log[end])/(Log[10]))}])//N
)

Clear[computeIaIbL]
computeIaIbL[f_,g_,lp0_,lpf_,κf_,pgoal_:20]:=Module[{α,Intf,Intfsqr,Intfoverp,lCbval},
Intf=intf[f,lp0,lpf,κf,pgoal];
Intfsqr=intfsqr[f,lp0,lpf,κf,pgoal];
Intfoverp=intfoverp[f,lp0,lpf,κf,pgoal];
lCbval=lCb[(Sqrt[exppTsqr[f,lp0,lpf,κf,pgoal]/Intf]),(Sqrt[mDsqr[Intfoverp,Nc,g]])];
(*Print["lCb = "<>ToString[lCbval]];*)
{4π (α^2)(Nc^2)lCbval (Intf+Intfsqr),4π (α^2)(Nc^2)lCbval2 Intfoverp,lCbval}//.{α->((g^2)/(4π))}]

Clear[logfinitlogp]
logfinitlogp[lp_,κ_,g_]:=(Log[((σ0)/(g^2)) ]-((Exp[2lp](1-(1-(ξ^2))(κ^2)))/(Qs^2)))//.{ξ->2,Qs->1}
(*Clear[logftinitlogp]
logftinitlogp[lp_,κ_,g_]:=(2lp+Log[((σ0)/(g^2)) ]-((Exp[2lp](1-(1-(ξ^2))(κ^2)))/(Qs^2)))//.{ξ->2,Qs->1}*)

Clear[writeParamfile];
writeParamFile:=(
Export["paramfile.dat",{OutputForm[logfnow],"g="<>ToString[gnow],"lp0="<>ToString[lp0]<>", lpf="<>ToString[lpf],"κf="<>ToString[κf],"τ0="<>ToString[τ0]<>", τF="<>ToString[τF]<>", ntau="<>ToString[ntau],"pgoal="<>ToString[pgoal],"max cell="<>ToString[maxcell]}]
)

nvals={{0,0,1},{0,0,2},{0,0,1}};
mvals={{0,1,0},{0,2,0},{1,2,1}};
momentsToCompute=DeleteDuplicates[Table[{nvals[[i,j]],mvals[[i,j]]},{i,Length[nvals]},{j,Length[nvals[[1]]]}]~Flatten~1];
allMoments=Table[{nvals[[i,j]],mvals[[i,j]]},{i,Length[nvals]},{j,Length[nvals[[1]]]}]
momentPositions=Table[Position[momentsToCompute,allMoments[[i,j]]][[1,1]],{i,Length[nvals]},{j,Length[nvals[[1]]]}]

Clear[solveIaIb]
solveIaIb[logfstart_,g_,Iboff_:False,starti_:1]:=(
Nc=3;
(*pgoal=5;*)
pgoal=10;
maxcell=(*5 (10^(-4))*)(10^(-3));
(*maxcellmax=(10^(-5));*)
Print["running new code!"];
Clear[Ia,Ib,inicond,τini,κ0];
replIb=If[Iboff,{Ib->0},{}];
myeqn={Derivative[0,0,1][logf][lp,κ,τ]+((κ ((-1+(κ^2)) Derivative[0,1,0][logf][lp,κ,τ]-κ Derivative[1,0,0][logf][lp,κ,τ]))/τ)== ((E^(-lp)) Ib (Derivative[1,0,0][logf][lp,κ,τ]+2 (1+Exp[logf[lp,κ,τ]](1+Derivative[1,0,0][logf][lp,κ,τ])))+Ia (E^(-2 lp))(-2 κ Derivative[0,1,0][logf][lp,κ,τ]-(-1+(κ^2)) (Derivative[0,2,0][logf][lp,κ,τ]+((Derivative[0,1,0][logf][lp,κ,τ])^2))+Derivative[1,0,0][logf][lp,κ,τ]+Derivative[2,0,0][logf][lp,κ,τ]+((Derivative[1,0,0][logf][lp,κ,τ])^2))),logf[lp,κ,τini]==inicond}//.replIb;
logfnow=logfstart;
writeParamFile;
{timeIaIb,{Ianow,Ibnow,Lnow}}=Timing[computeIaIbL[Exp[logfnow],g,lp0,lpf,κf,pgoal]];
(*Print["first time IaIb = "<>ToString[timeIaIb]];*)
Print["L0="<>ToString[Lnow]];
(*c=(Exp[logfnow]//.{lp->lp0,κ->0})Exp[lp0];
fluxnow=If[Ibnow (c^2)>Ianow c,-Ibnow (c^2),0];*)
(*fluxnow=-Ibnow (c^2);*)
alln={};
allIa={};
savefs={};
energies={};
pls={};
ΛTs={};
ΛLs={};
(*fluxes={};*)
If[!Iboff,allIb={};];
Do[
(*Print[ToString[i]];*)
τi=τvals[[i]];
τf=τvals[[i+1]];
Clear[inifun];
(*inifun[lp_,κ_]:=Piecewise[{{logfnow,lp≤lpf&&lp≥lp0&&κ≥0&&κ≤κf}},0];*)
inifun[lp_,κ_]:=Piecewise[{{logfnow//.κ->Abs[κ],lp≤lpf&&lp≥lp0&&Abs[κ]≤1},{logfnow//.κ->Abs[κ]//.lp->lpf,lp>lpf&&Abs[κ]≤1},{logfnow//.κ->Abs[κ]//.lp->lp0,lp<lp0&&Abs[κ]≤1},{logfnow//.κ->1,lp≤lpf&&lp≥lp0&&Abs[κ]>1},{logfnow//.κ->1//.lp->lpf,lp>lpf&&Abs[κ]>1},{logfnow//.κ->1//.lp->lp0,lp<lp0&&Abs[κ]>1}},-20];
eqnow=myeqn//.{Ia->Ianow,Ib->Ibnow,τini->τi,inicond->inifun[lp,κ]};
{timesol,ndsolFP}=Timing[NDSolve[eqnow,logf,{lp,lp0,lpf},{κ,0,κf},{τ,τi,τf},Method->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->maxcell}}}]];
(*If[i<25,Print["time sol = "<>ToString[timesol]]];*)
(*While[(logf//.ndsolFP[[1]])[[1,3,2]]≠τf&&maxcell≤maxcellmax,
ndsolFP=NDSolve[eqnow,logf,{lp,lp0,lpf},{κ,0,κf},{τ,τi,τf},Method->{"MethodOfLines","SpatialDiscretization"->{"FiniteElement","MeshOptions"->{"MaxCellMeasure"->maxcell}}}];
maxcell=maxcell/2;
Print["max cell revised to "<>ToString[maxcell]];
];*)
logfnow=logf[lp,κ,τf]//.ndsolFP[[1]];
AppendTo[alln,Table[(2/(((2π))^2))NIntegrate[Exp[logfnow] (Exp[(3+m+n)lp](((1-(κ^2)))^(m/2))((Abs[κ])^n))//.{n->moment[[1]],m->moment[[2]]},{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal],{moment,momentsToCompute}]];
(*Print["time moments = "<>ToString[timemoments]];*)
AppendTo[energies,(2/(((2π))^2))NIntegrate[Exp[4lp] Exp[logfnow],{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal]];
AppendTo[pls,(2/(((2π))^2))NIntegrate[Exp[4lp]((Abs[κ])^2)Exp[logfnow],{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal]];
AppendTo[ΛLs,(2/(((2π))^2))NIntegrate[Exp[6lp]((Abs[κ])^2)Exp[logfnow],{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal]];
AppendTo[ΛTs,(2/(((2π))^2))NIntegrate[Exp[6lp](1-(κ^2))Exp[logfnow],{lp,lp0,lpf},{κ,0,κf},PrecisionGoal->pgoal]];
{timeIaIb,{Ianow,Ibnow,Lnow}}=Timing[computeIaIbL[Exp[logfnow],g,lp0,lpf,κf,pgoal]];
(*If[i<25,Print["time IaIb = "<>ToString[timeIaIb]]];*)
AppendTo[allIa,Ianow];
If[!Iboff,AppendTo[allIb,Ibnow];];
savePositions=Position[(#≥τi&&#<τf)&/@τsave,True];
If[savePositions≠{},AppendTo[savefs,logf[lp,κ,τsave[[savePositions[[1]]]]]//.ndsolFP[[1]]]];
If[Mod[i,25]==0,Put[If[!Iboff,{logfnow,τvals[[;;(i+1)]],ΛLs,ΛTs,energies,pls,savefs,alln,allIa,allIb},{logfnow,τvals[[;;(i+1)]],ΛLs,ΛTs,energies,pls,savefs,alln,allIa}],"intermed_output_"<>ToString[i]]];,{i,starti,Length[τvals]-1}];)

Clear[averageMoments]
averageMoments[τs_,ns_,di_:5]:=(
τavg={};
navg={};
Do[AppendTo[τavg,Mean[τs[[i;;i+di-1]]]];
AppendTo[navg,Mean[ns[[i;;i+di-1,;;]]]];,{i,1,Length[τs]-di,di}];
{τavg,navg}
)
Clear[solveExponents]
solveExponents[allmoments_]:=(
logn=Table[Interpolation[Thread[{τvals[[;;-2]],Log[allmoments[[;;,i]]]}](*,InterpolationOrder->1*)],{i,Length[momentsToCompute]}];
(*Print["len(logn)="<>ToString[Length[logn]]];*)
sols=Table[Solve[Thread[τ D[logn[[#]][τ],τ]&/@momentPositions[[i]]==(α-(m+2)β-(n+1)γ)//.{m->mvals[[i]],n->nvals[[i]]}],{α,β,γ}][[1]],{i,Length[momentPositions]}];
(*Print["len(sols)="<>ToString[Length[sols]]];*)
{α,β,γ}//.#&/@sols
)

Clear[fsc]
fsc[pT_,pz_,A_,B_,CC_]:=A Exp[-(1/2)(((pT^2)/(B^2))+((pz^2)/(CC^2)))]

(* compute initial conditions for A,B,C from matching to scaling form *)
Clear[numericalMatch];
numericalMatch[f_,lp0_,lpf_,κf_,pgoal_:20]:=Module[{niformula,pzsqrformula,pTsqrformula,ninum,pzsqrnum,pTsqrnum,sols,ispositive,possol},(
niformula=Assuming[{A>0,B>0,CC>0},Integrate[(1/(((2π))^3))2π pT fsc[pT,pz,A,B,CC],{pT,0,Infinity},{pz,-Infinity,Infinity}]];
pzsqrformula=Assuming[{A>0,B>0,CC>0},Integrate[(1/(((2π))^3))2π pT (pz^2) fsc[pT,pz,A,B,CC],{pT,0,Infinity},{pz,-Infinity,Infinity}]]/niformula;
pTsqrformula=Assuming[{A>0,B>0,CC>0},Integrate[(1/(((2π))^3))2π pT (pT^2) fsc[pT,pz,A,B,CC],{pT,0,Infinity},{pz,-Infinity,Infinity}]]/niformula;
ninum=intf[f,lp0,lpf,κf,pgoal];
pzsqrnum=exppzsqr[f,lp0,lpf,κf,pgoal]/ninum;
pTsqrnum=exppTsqr[f,lp0,lpf,κf,pgoal]/ninum;
sols=Solve[{niformula==ninum,pTsqrformula==pTsqrnum,pzsqrformula==pzsqrnum},{A,B,CC}];
ispositive=Position[Table[AllTrue[Thread[(({A,B,CC}//.thissol)>0)],#&],{thissol,sols}],True][[1,1]];
sols[[ispositive]]
)]