SetDirectory["/home/ramkishor/pencil-code/mathematica"];
Needs["pc`"]
sim1=SetDirectory["/home/ramkishor/cvs_axel/sayan/GW/LowFreq1024sig1_k01_kf10b_rep2"];
sim2=SetDirectory["/home/ramkishor/cvs_axel/sayan/GW/LowFreq1024sig0_k01_kf10b2_rep2"];

Power Spectrum

{th,spechdot}=read1D[sim1,"power_GWs.dat"];
{tn,specnhdot}=read1D[sim2,"power_GWs.dat"];
{th,specstressh}=read1D[sim1,"power_Str.dat"];
{tn,specstress}=read1D[sim2,"power_Str.dat"];

Spectrum of (Overscript[h, (OverDot[~])])

tstep={1,50,100,200,349};
tstep2={1,50,75,100,197};
Show[{ListLogLogPlot[{spechdot[[tstep[[1]]]],spechdot[[tstep[[2]]]],spechdot[[tstep[[3]]]],spechdot[[tstep[[4]]]],spechdot[[tstep[[5]]]]},Joined->True,PlotStyle-> {Black,Black,Blue,Red,Orange},Frame->True,FrameStyle->Directive[Black,8],FrameLabel->{Style
["k",Black,8],Style["Sp((OverDot[h]))",Black,8]}],ListLogLogPlot[{specnhdot[[tstep2[[1]]]],specnhdot[[tstep2[[2]]]],specnhdot[[tstep2[[3]]]],specnhdot[[tstep2[[4]]]],specnhdot[[tstep2[[5]]]]},Joined->True,PlotStyle-> {Black,Black,Blue,Red,Orange},Frame->True,FrameStyle->Directive[Black,8],FrameLabel->{Style["k",Black,8],Style["Sp((OverDot[h]))",Black,8]}](*,ListLogLogPlot[{spechdot3[[tstep3[[1]]]],spechdot3[[tstep3[[2]]]],spechdot3[[tstep3[[3]]]],spechdot3[[tstep3[[4]]]],spechdot3[[tstep3[[5]]]]},Joined->True,PlotStyle-> {Black,Black,Blue,Red,Orange},Frame->True,FrameStyle->Directive[Black,8],FrameLabel->{Style["k",Black,8],Style["Sp((OverDot[h]))",Black,8]}]*)}];

Analytical Case

partn[t_,t1_,k_,β_]=Integrate[Cos[k(t-t1)](t1)^(β-1),t1];
partn2only[t_,t0_,k_,β_]=(partn[t,t0,k,β]-partn[t,1,k,β]);
partn2onlyh[t_,t0_,k_]:=NIntegrate[Cos[k(t-t1)](((t1^(-1))/(((1+(((t1-0.2))^20)))^(-(1/20)*(5/24))))),{t1,1,t0}];
partn2onlynh[t_,t0_,k_]:=NIntegrate[Cos[k(t-t1)](((t1^(-7/6))/(((1+(((t1-1))^10)))^((1/10)*(5/14))))),{t1,1,t0}];

FunctionEnvelope[f_,{t_,a_,b_},n_: 40]:=Module[{seeds,x,y,points,progress=0,tempf,union},seeds=Rescale[Range[0,1,1/n]+1/(2 n),{0,1},{a,b}];
points=Last[Last[Reap[Monitor[Do[Quiet@Check[progress++;{y,x}=FindMaximum[Abs[f],{t,x0}];
x=t/.x;If[a≤x≤b,Sow[{x,y}]],Null],{x0,seeds}],ProgressIndicator[progress,{0,Length[seeds]}]]]]];
union[]:=points=Union[points,SameTest->(Abs[First[#1]-First[#2]]/Replace[Max[Abs[First[#1]],Abs[First[#2]]],u_/;u==0:>1]<10^-6&)];
union[];
tempf=Interpolation[points];
points=Quiet[Join[{{a,tempf[a]}},points,{{b,tempf[b]}}]];
union[];
Interpolation[points]]

t=100;
β1=1/3;
β2=-1/4;

th[[349]]

Clear[β1,t,g1,f1,h1]
t=100;
β1=(1/3);
g1=FunctionEnvelope[(partn2only[10 t,10 t,k,β1])^2,{k,0.001,1},500];
f1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,1,30},200];
h1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,10,1000},200];

g1=FunctionEnvelope[(partn2only[2.5 t,2.5 t,k,β1])^2,{k,10^-5,1},500];
f1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,1,30},200];
h1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,10,1000},200];
g2=FunctionEnvelope[(partn2only[2.5 t,2.5 t,k,β2])^2,{k,10^-5,1},500];
f2=FunctionEnvelope[(partn2only[t,t,k,β2])^2,{k,1,30},200];
h2=FunctionEnvelope[(partn2only[t,t,k,β2])^2,{k,10,1000},200];
g7h=FunctionEnvelope[(partn2onlyh[t,t,k])^2,{k,10^-5,1},100];
f7h=FunctionEnvelope[(partn2onlyh[t,t,k])^2,{k,1,30},200];
h7h=FunctionEnvelope[(partn2onlyh[t,t,k])^2,{k,10,1000},200];
g7nh=FunctionEnvelope[(partn2onlynh[t,t,k])^2,{k,10^-5,1},100];
f7nh=FunctionEnvelope[(partn2onlynh[t,t,k])^2,{k,1,30},200];
h7nh=FunctionEnvelope[(partn2onlynh[t,t,k])^2,{k,10,1000},200];

Visual check for the goodness of the envelope

pg1=LogLogPlot[g1[k],{k,0.01,1},PlotStyle->{Black,DotDashed}];
pf1=LogLogPlot[f1[k],{k,1,20},PlotStyle->{Black,DotDashed}];
ph1=LogLogPlot[h1[k],{k,20,50},PlotStyle->{Black,DotDashed}];
Show[LogLogPlot[(partn2only[t,t,k,β1])^2,{k,0.1,30}],pg1,pf1,ph1]
pg2=LogLogPlot[g2[k],{k,0.01,1},PlotStyle->{Black,DotDashed}];
pf2=LogLogPlot[f2[k],{k,1,20},PlotStyle->{Black,DotDashed}];
ph2=LogLogPlot[h2[k],{k,20,50},PlotStyle->{Black,DotDashed}];
Show[LogLogPlot[(partn2only[t,t,k,β2])^2,{k,0.1,30}],pg2,pf2,ph2]

SetDirectory["/home/ramkishor/Desktop/GWspectra_from_mhd_turbulence (Version 1229)/mathematica/figures"]

Clear[t,kph]
kph[t_]=10(((t))^(-1/3));
n=1;
magspech[k_,t_]=((((k/kph[t]))^4)/(((1+(((k/kph[t]))^(n(4+5/3)))))^(1/n)));
stressh[k_,t_]:=NIntegrate[magspech[k1,t]magspech[√((k^2)+(k1^2)-2k k1 μ),t]((k^2)/((k^2)+(k1^2)-2k k1 μ))((1+(μ^2))(1+(((k-k1 μ))^2)/((k^2)+(k1^2)-2k k1 μ))+4μ((k-k1 μ)/((((k^2)+(k1^2)-2k k1 μ))^(1/2))))(((2π)/(4π))),{k1,0.001,(10^5)},{μ,-1,1}];
findLogDivisions[{xmin_,xmax_},n_Integer]:=10^FindDivisions[Log10@{xmin,xmax},n]
interval=findLogDivisions[{0.1,50},100];
dd=Dimensions[interval];
plotsh1=Table[{interval[[i]],stressh[interval[[i]],1]},{i,1,dd[[1]]}];
fitfunc=Interpolation[plotsh1];
bb=Interpolation[Table[{((i-1)/10),(((((th[[1]])/6)))^2)specstressh[[1,i]]},{i,1,512}]];
plotsh2=Table[{interval[[i]],bb[interval[[i]]]},{i,1,dd[[1]]}];
fitfunc1h=Interpolation[plotsh2];
tokh[x_]=(10^(-8))*fitfunc[x];
tokh1[x_]=0.43fitfunc1h[x];

β1=1/3;
β2=-1/4;
pgh=LogLogPlot[36tokh[k]g7h[k],{k,0.1,1},PlotStyle->{Red,Dotted},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from Eq.(A5)",10]},LegendMarkerSize->15],{0.6,0.48}]];
pfh=LogLogPlot[36tokh[k]f7h[k],{k,1,20},PlotStyle->{Red,Dotted}];
phh=LogLogPlot[36tokh[k] h7h[k],{k,20,50},PlotStyle->{Red,Dotted}];
pgh1=LogLogPlot[36tokh1[k]g7h[k],{k,0.1,1},PlotStyle->{Green,Dashed},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from simulation",10]},LegendMarkerSize->15],{0.615,0.38}]];
pfh1=LogLogPlot[36tokh1[k]f7h[k],{k,1,20},PlotStyle->{Green,Dashed}];
phh1=LogLogPlot[36tokh1[k] h7h[k],{k,20,50},PlotStyle->{Green,Dashed}];
pg1=LogLogPlot[36tokh[k]g1[k],{k,0.1,1},PlotStyle->{Black,DotDashed},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from Eq.(A5)",10]},LegendMarkerSize->15],{0.6,0.28}]];
pf1=LogLogPlot[36tokh[k]f1[k],{k,1,20},PlotStyle->{Black,DotDashed}];
ph1=LogLogPlot[36tokh[k] h1[k],{k,20,50},PlotStyle->{Black,DotDashed}];
plot1=Show[{ListLogLogPlot[Table[{((i-1)/10),spechdot[[349,i]]},{i,1,512,1}],Joined->True,PlotStyle-> {Blue},Frame->True,FrameStyle->Directive[Black,Bold,12],FrameLabel->{Style["k",Black,Bold,12],Style["Sp((OverDot[h]))",Black,Bold,12]}],pgh,pfh,phh,pgh1,pfh1,phh1,pg1,pf1,ph1},PlotRange->All,Epilog->{Text["(a)",{-2,-18.5}],Arrow@Line[{{-2.3,-20},{-1.8,-20.9}},VertexColors->{Directive[Green,Opacity@1],Green}],Text["from sim.",{-2.1,-19.8}],Arrow@Line[{{-1.3,-22.2},{-2.1,-21.3}},VertexColors->{Directive[Black,Opacity@1],Black}],Text["single power law",{-1.3,-22.4}],Arrow@Line[{{-0.8,-21.4},{-1.1,-20.7}},VertexColors->{Directive[Red,Opacity@1],Red}],Text["double power law",{-0.6,-21.7}]}];
Export["LowFreq1024sig1_comparison_broken.eps",plot1];

Clear[t,kp,plotsh1]
kp[t_]=10(((t))^(-1/4));
n=1;
magspec[k_,t_]=((((k/kp[t]))^4)/(((1+(((k/kp[t]))^(n(4+5/3)))))^(1/n)));
stress[k_,t_]:=NIntegrate[magspec[k1,t]magspec[√((k^2)+(k1^2)-2k k1 μ),t]((k^2)/((k^2)+(k1^2)-2k k1 μ))(1+(μ^2))(1+(((k-k1 μ))^2)/((k^2)+(k1^2)-2k k1 μ))(((2π)/(4π))),{k1,0.001,(10^5)},{μ,-1,1}];
findLogDivisions[{xmin_,xmax_},n_Integer]:=10^FindDivisions[Log10@{xmin,xmax},n]
interval=findLogDivisions[{(10^(-4)),5000},200];
dd=Dimensions[interval];
plotsh1=Table[{interval[[i]],stress[interval[[i]],1]},{i,1,dd[[1]]}];
fitfunc1=Interpolation[plotsh1];

aa=Interpolation[Table[{((i-1)/10),(((((tn[[1]])/6)))^2)specstress[[1,i]]},{i,1,512}]];
plots1=Table[{interval[[i]],aa[interval[[i]]]},{i,1,dd[[1]]}];
fitfunc2=Interpolation[plots1];
tok[x_]=(10^(-8))fitfunc1[x];
tok1[x_]=0.58fitfunc2[x];

Show[{(*ListLogLogPlot[Table[{(i/10),spechdot[[tstep[[3]],i]]},{i,0,511,1}],Joined->True,PlotStyle-> {Black,Black,Blue,Red,Orange},Frame->True,FrameStyle->Directive[Black,8],FrameLabel->{Style["k",Black,8],Style["Sp((OverDot[h]))",Black,8]}],*)ListLogLogPlot[Table[{((i-1)/10),specstress[[1,i]]},{i,1,512,1}],Joined->True,PlotStyle-> {Blue},Frame->True,FrameStyle->Directive[Black,8],FrameLabel->{Style["k",Black,8],Style["Sp((OverDot[h]))",Black,8]}],LogLogPlot[c3[k],{k,0.1,5000},PlotStyle->{Black}]},PlotRange->All]

β1=1/3;
β2=-1/4;
pgnh=LogLogPlot[36 tok[k] g7nh[k],{k,0.1,1},PlotStyle->{Red,Dotted},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from Eq.(A3)",10]},LegendMarkerSize->15],{0.6,0.58}]];
pfnh=LogLogPlot[36tok[k] f7nh[k],{k,1,20},PlotStyle->{Red,Dotted}];
phnh=LogLogPlot[36tok[k]h7nh[k],{k,20,50},PlotStyle->{Red,Dotted}];
pgnh1=LogLogPlot[36tok1[k] g7nh[k],{k,0.1,1},PlotStyle->{Green,Dashed},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from simulation",10]},LegendMarkerSize->15],{0.615,0.48}]];
pfnh1=LogLogPlot[36tok1[k] f7nh[k],{k,1,20},PlotStyle->{Green,Dashed}];
phnh1=LogLogPlot[36tok1[k]h7nh[k],{k,20,44},PlotStyle->{Green,Dashed}];
pg2=LogLogPlot[36tok[k] g2[k],{k,0.1,1},PlotStyle->{Black,DotDashed},PlotLegends->Placed[LineLegend[{Style["|(Subscript[(OverTilde[T]), 0])(k)| from Eq.(A3)",10]},LegendMarkerSize->15],{0.6,0.38}]];
pf2=LogLogPlot[36tok[k] f2[k],{k,1,20},PlotStyle->{Black,DotDashed}];
ph2=LogLogPlot[36tok[k]h2[k],{k,20,50},PlotStyle->{Black,DotDashed}];
plot1=Show[{ListLogLogPlot[Table[{((i-1)/10),specnhdot[[197,i]]},{i,1,512,1}],Joined->True,PlotRange->{{0.06,70},{(10^(-10)),2*(10^(-8))}},PlotStyle-> {Blue},Frame->True,FrameStyle->Directive[Black,Bold,12],FrameLabel->{Style["k",Black,Bold,12],Style["Sp((OverDot[h]))",Black,Bold,12]},Epilog->{Text["(b)",{-2,-18.5}],Arrow@Line[{{-2.3,-20},{-2,-21.1}},VertexColors->{Directive[Green,Opacity@1],Green}],Text["from sim.",{-2.1,-19.8}],Arrow@Line[{{-1.3,-22.4},{-2.1,-21.4}},VertexColors->{Directive[Black,Opacity@1],Black}],Text["single power law",{-1.3,-22.6}],Arrow@Line[{{-0.8,-21.4},{-1.3,-20.4}},VertexColors->{Directive[Red,Opacity@1],Red}],Text["double power law",{-0.6,-21.7}]}],pgnh,pfnh,phnh,pgnh1,pfnh1,phnh1,pg2,pf2,ph2}];
Export["LowFreq1024sig0_comparison_broken.eps",plot1];