ClearAll["Global`*"]
Needs["pc`"]
pcInitialize[]

dir="/Users/hzhou/Work/pencilCode/runs/projects/ssdy/";
runs=dir<>#&/@{"dp-2dq1"};

sim=runs[[1]];
{t,Em}=read1D[sim,"power_mag.dat"];
{t,Ek}=read1D[sim,"power_kin.dat"];
t
t//Length
0.4*4*2Ek[[1]]//Total//Sqrt

Clear[readSpecFluxPQ,readSpecFlux]
readSpecFluxPQ[sim_]:=Module[{dp,dq,nk,nlk,grid},
dp=readParamNml[sim,"run.in","SPECFLUX_DP"];
dq=readParamNml[sim,"run.in","SPECFLUX_DQ"];
nk=getParam[sim,"Nx"]/2;

nlk[d_]:=Floor@If[d>0,(nk-1)/d+1,Log[-d,nk-1]+1];
grid[d_]:=If[d>0,d*Range[0,nlk[d]-1],(-d)^Range[0,nlk[d]-1]];

Round@{grid[dp],grid[dq]}
]
readSpecFlux[sim_,file_]:=Module[{pq,grid,t,tra},
pq=readSpecFluxPQ[sim];
grid=Table[{pp,qq},{qq,pq[[2]]},{pp,pq[[1]]}]//Flatten[#,1]&;

{t,tra}=read1D[sim,file,Times@@(Length/@pq)];
tra=Flatten/@Transpose[{grid,#}]&/@tra;

Return[{t,tra}]
]

{t,tra}=readSpecFlux[sim,"power_mag_E_transfer_magE.dat"];
pq=readSpecFluxPQ[sim]
Length[pq[[1]]]

k=pq[[1]];
knu=461;
style={ImageSize->250,ImagePadding->{{80,15},{50,5}}};
labelk={{0.001,"0.001"},{0.01,"0.010"},{0.1,"0.100"},{1,"1.000"}};
pq=readSpecFluxPQ[sim];
Table[ListLogLinearPlot[Table[
With[{data=Apply[Association,tra[[i]]/.{x_,y_,z_}:>Round[{x,y}]->z]},{1/knu,1}*#&/@Table[{qq,10^12*data[Round@{kk,qq}]},{qq,pq[[2]]}]]
,{i,Length[t]}],
PlotLegends->Placed[LineLegend[{"k="<>ToString[kk]<>"="<>ToString[kk/knu//N//SetPrecision[#,2]&]<>"(Subscript[k, ν])"}],{Right,Top}],
PlotStyle->(Hue/@Rescale[t]),
FrameLabel->{"p/(Subscript[k, ν])","(10^12)\[Cross](Subscript[T, kp])"},
GridLines->{{4/knu},None},
FrameTicks->{{Automatic,Automatic},{labelk,Automatic}},
style
],{kk,k}]//Partition[#,UpTo@4]&//Grid

k=pq[[1]];
knu=461;
style={ImageSize->250,ImagePadding->{{80,15},{50,5}}};
labelk={{0.001,"0.001"},{0.01,"0.010"},{0.1,"0.100"},{1,"1.000"}};
pq=readSpecFluxPQ[sim];
Table[ListPlot[Table[
With[{data=Apply[Association,tra[[i]]/.{x_,y_,z_}:>Round[{x,y}]->z]},
Accumulate@Table[10^12*data[Round@{kk,qq}],{qq,pq[[2]]}]]
,{i,Length[t]}],
DataRange->{1/knu,511/knu},
PlotLegends->Placed[LineLegend[{"k="<>ToString[kk]<>"="<>ToString[kk/knu//N//SetPrecision[#,2]&]<>"(Subscript[k, ν])"}],{Right,Bottom}],
PlotStyle->(Hue/@Rescale[t]),
FrameLabel->{"p/(Subscript[k, ν])","(10^12)\[Cross]∑(Subscript[ , p'=1]^p)(Subscript[T, kp'])Δp'"},
GridLines->{{4/knu},None},
FrameTicks->{{Automatic,Automatic},{Automatic,Automatic}},
style
],{kk,k}]//Partition[#,UpTo@4]&//Grid

tra={Mean[tra]};
k={43,110,281};
knu=461;
i=1;
plot=With[{specK=Ek[[i]],spec=Em[[i]],data=Apply[Association,tra[[i]]/.{x_,y_,z_}:>Round[{x,y}]->z]},
style={ImageSize->350};
labelk={{0.001,"0.001"},{0.01,"0.010"},{0.1,"0.100"},{1,"1.000"}};
Grid[{{

Show[ListLogLogPlot[{Rest[specK],Rest[spec]},
DataRange->{0,512/knu},
PlotStyle->{Black,Red},
FrameLabel->{"k/(Subscript[k, ν])","Energy spectra"},
FrameTicks->{{pcTicks["10^i"][[1;;-1;;4]],Automatic},{labelk,Automatic}},
PlotLegends->Placed[LineLegend[{"(Subscript[E, "K"])","(Subscript[E, "M"])"}],{Right,Top}],
GridLines->{{4/knu,(kpeak=Position[spec,Max[spec]][[1,1]]-1)/knu},None},
Epilog->{Inset[Style["(a)",pcLabelStyle],Scaled@{0.9,0.1}],
Inset[Style["∝(k^(-5/3))",pcLabelStyle],Scaled@{0.5,0.7}],
Inset[Style["∝(k^(3/2))",pcLabelStyle/.Black->Red],Scaled@{0.5,0.32}]}
],
LogLogPlot[5*10^(-5)(kk/10*knu)^(-5/3),{kk,10/knu,90/knu},PlotStyle->{Black,Dashed}],
LogLogPlot[10^(-13)(kk/10*knu)^(3/2),{kk,5/knu,60/knu},PlotStyle->{Red,Dashed}],
style
],

Print["kpeak=",kpeak];
pq=readSpecFluxPQ[sim];
ListLogLinearPlot[Table[Table[{qq/knu,10^If[kk<280,12,11]*data[Round@{kk,qq}]},{qq,pq[[2]]}],{kk,k}],
PlotLegends->Placed[LineLegend[{"(10^12)\[Cross](Subscript[T, 43p])","(10^12)\[Cross](Subscript[T, 110p])","(10^11)\[Cross](Subscript[T, 281p])"},Spacings->0.15],{Right,Top}],
PlotStyle->{Black,Red,Blue},
FrameLabel->{"p/(Subscript[k, ν])","Normalized (Subscript[T, kp])"},
GridLines->{{4/knu},None},
Epilog->{Inset[Style["(b)",pcLabelStyle],Scaled@{0.9,0.1}]},
FrameTicks->{{Automatic,Automatic},{labelk,Automatic}},
style
]},

{ListPlot[Table[10^11*Accumulate[data[Round@{kk,#}]&/@pq[[2]]],{kk,k}],
DataRange->{1,511}/knu,
PlotLegends->Placed[LineLegend["k/(Subscript[k, 1])="<>ToString[#]&/@k,Spacings->0.15],{Left,Top}],
PlotRange->{{Automatic,250/knu},All},
PlotStyle->{Black,Red,Blue},
FrameLabel->{"p/(Subscript[k, ν])","(10^11)\[Cross]∑(Subscript[ , p'=1]^p)(Subscript[T, kp'])Δp'"},
GridLines->{{4/knu},None},
Epilog->{Inset[Style["(c)",pcLabelStyle],Scaled@{0.9,0.1}]},
style
],

Show[ListLogLogPlot[Table[{kk/knu,10^10*Total[data[Round@{kk,#}]&/@pq[[2]]]},{kk,pq[[1]]}],
FrameLabel->{"k/(Subscript[k, ν])","(10^10)\[Cross]∑(Subscript[ , p])(Subscript[T, kp])Δp"},
Epilog->{Inset[Style["(d)",pcLabelStyle],Scaled@{0.9,0.1}],
Inset[Style["∝(k^(3/2))",pcLabelStyle],Scaled@{0.5,0.72}]},
GridLines->{{4/knu,kpeak/knu},None},
FrameTicks->{{pcTicks["10^i"],Automatic},{labelk,Automatic}}
],
LogLogPlot[0.01(kk/5*knu)^1.5,{kk,4/knu,kpeak/knu},PlotStyle->{Dashed,Black}],
style
]}},Spacings->{0.5,-1.5}
]
]

If[True,Export["/Users/hzhou/Dropbox/Apps/Overleaf/Kazantsev-Subinertial/fig/tra.eps",plot],"No output"]