Analytical Case

kpn[t_,α_]=k0((((t/1)))^α);
partn[t_,t1_,k_,β_]=Integrate[Cos[k(t-t1)](t1)^(β-1),t1];
partn2only[t_,t0_,k_,β_]=(partn[t,t0,k,β]-partn[t,1,k,β]);
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;
norm=10^-5;
β1=1;
β2=2/3;
β3=1/3;
β4=1/6;
β5=-7/6;
β6=-3/4;
β7=-1/4
kp1=20;
kp=200;
klower=0.05;
kupper=1000;

g1=FunctionEnvelope[(partn2only[2.5 t,2.5 t,k,β1])^2,{k,10^-5,1},500];
g2=FunctionEnvelope[(partn2only[2.5 t,2.5 t,k,β2])^2,{k,10^-5,1},500];
g3=FunctionEnvelope[(partn2only[2.5 t,2.5 t,k,β3])^2,{k,10^-5,1},100];
g4=FunctionEnvelope[(partn2only[t,t,k,β4])^2,{k,10^-5,1},100];
g5=FunctionEnvelope[(partn2only[t,t,k,β5])^2,{k,10^-5,1},100];
g6=FunctionEnvelope[(partn2only[t,t,k,β6])^2,{k,10^-5,1},100];
g7=FunctionEnvelope[(partn2only[t,t,k,β7])^2,{k,10^-5,1},100];
f1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,1,30},200];
f2=FunctionEnvelope[(partn2only[t,t,k,β2])^2,{k,1,30},200];
f3=FunctionEnvelope[(partn2only[t,t,k,β3])^2,{k,1,30},200];
f4=FunctionEnvelope[(partn2only[t,t,k,β4])^2,{k,1,30},200];
f5=FunctionEnvelope[(partn2only[t,t,k,β5])^2,{k,1,30},200];
f6=FunctionEnvelope[(partn2only[t,t,k,β6])^2,{k,1,30},200];
f7=FunctionEnvelope[(partn2only[t,t,k,β7])^2,{k,1,30},200];
h1=FunctionEnvelope[(partn2only[t,t,k,β1])^2,{k,10,1000},200];
h2=FunctionEnvelope[(partn2only[t,t,k,β2])^2,{k,10,1000},200];
h3=FunctionEnvelope[(partn2only[t,t,k,β3])^2,{k,10,1000},200];
h4=FunctionEnvelope[(partn2only[t,t,k,β4])^2,{k,10,1000},200];
h5=FunctionEnvelope[(partn2only[t,t,k,β5])^2,{k,10,1000},200];
h6=FunctionEnvelope[(partn2only[t,t,k,β6])^2,{k,10,1000},200];
h7=FunctionEnvelope[(partn2only[t,t,k,β7])^2,{k,10,1000},200];

Clear[t,kp]
kp[t_]=100(((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 μ))(1/(((2π))^3)),{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];
c3[x_]=(6/100000)*1.5fitfunc1[x];

t=100;
kp=200;
pg1=LogLogPlot[c3[k] g1[k],{k,10^-4,1},PlotRange->{{5*10^-5,2000},{10^-17,10^-7}},PlotStyle->{Blue},Epilog->{{Arrow[{{-2,-36},{0,-38}}],Text["1/(Subscript[t, 0])",{-2.5,-36}],Arrow[{{-6,-36},{-4.6,-38}}],Text["1/t",{-6.5,-35}]},Arrow[{{6.5,-36},{5.3,-38}}],Text["2 (Subscript[k, p])",{6.5,-35}]},GridLines->{{1/t,1,kp},None},GridLinesStyle->{{Gray,Thick},None},Frame->True,FrameStyle->Directive[Black,15,Bold],FrameTicks->{{{1,kp},None},{{0.01,1,kp},None}},FrameLabel->{Style["k",Black,Bold,15],Style["Sp((OverDot[h]))",Black,Bold,15]},PlotLegends->Placed[LineLegend[{Style["p=1",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg2=LogLogPlot[c3[k] g2[k],{k,10^-4,1},PlotStyle->{Dashed,Blue},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=2/3",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg3=LogLogPlot[c3[k]g3[k],{k,10^-4,1},PlotStyle->{Black},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=1/3",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg4=LogLogPlot[c3[k] g4[k],{k,10^-4,1},PlotStyle->{Dashed,Red},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=0",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg5=LogLogPlot[c3[k]g5[k],{k,10^-4,1},PlotStyle->{Dashed,Green},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=-7/6",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg6=LogLogPlot[c3[k]g6[k],{k,10^-4,1},PlotStyle->{Green},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=-3/4",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pg7=LogLogPlot[c3[k]g7[k],{k,10^-4,1},PlotStyle->{Orange},GridLines->{None,None},PlotLegends->Placed[LineLegend[{Style["p=-1/4",8],Style["t=4",8],Style["t=10",8]},LegendMarkerSize->12],{0.68,0.5}]];
pf1=LogLogPlot[c3[k] f1[k],{k,1,20},PlotStyle->{Blue}];
pf2=LogLogPlot[c3[k] f2[k],{k,1,20},PlotStyle->{Dashed,Blue}];
pf3=LogLogPlot[c3[k] f3[k],{k,1,20},PlotStyle->{Black}];
pf4=LogLogPlot[c3[k]f4[k],{k,1,20},PlotStyle->{Dashed,Red}];
pf5=LogLogPlot[c3[k]f5[k],{k,1,20},PlotStyle->{Dashed,Green}];
pf6=LogLogPlot[c3[k] f6[k],{k,1,20},PlotStyle->{Green}];
pf7=LogLogPlot[c3[k] f7[k],{k,1,20},PlotStyle->{Orange}];
ph1=LogLogPlot[c3[k] h1[k],{k,20,1000},PlotStyle->{Blue}];
ph2=LogLogPlot[c3[k] h2[k],{k,20,1000},PlotStyle->{Dashed,Blue}];
ph3=LogLogPlot[c3[k] h3[k],{k,20,1000},PlotStyle->{Black}];
ph4=LogLogPlot[c3[k] h4[k],{k,20,1000},PlotStyle->{Dashed,Red}];
ph5=LogLogPlot[c3[k]h5[k],{k,20,1000},PlotStyle->{Dashed,Green}];
ph6=LogLogPlot[c3[k] h6[k],{k,20,1000},PlotStyle->{Green}];
ph7=LogLogPlot[c3[k] h7[k],{k,20,1000},PlotStyle->{Orange}];

Plot6=LogLogPlot[0.4*10^-8 k^0.5,{k,1.5*10^-2,0.6},PlotStyle->{Gray},PlotLabels->Placed[{"(k^(1/2))"},{Scaled[1],Above}]];
Plot7=LogLogPlot[0.4*10^-8 k,{k,10^-2,0.3},PlotStyle->{Gray},PlotLabels->Placed[{"k"},{Scaled[1],{5,0.2}}]];
Plot8=LogLogPlot[0.5*10^-8 k^2,{k,10^-2,0.3},PlotStyle->{Gray},PlotLabels->Placed[{"(k^2)"},{Scaled[1],{2,0.5}}]];

plott=Show[{pg1,pf1,ph1,pg2,pf2,ph2,pg3,pf3,ph3,pg4,pf4,ph4,pg7,pf7,ph7,pg6,pf6,ph6,pg5,pf5,ph5,Plot6,Plot7,Plot8}];

SetDirectory["/home/ramkishor/Desktop/GWspectra_from_mhd_turbulence (Version 1229)/mathematica/figures"];
Export["sp_hdot_for_different_p.eps",plott];