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

{t1,specgw1}=read1D[sim1,"power_GWs.dat"];
{tmag1,specmag1}=read1D[sim1,"power_mag.dat"];
{tkin1,speckin1}=read1D[sim1,"power_kin.dat"];
{tstr1,specstr1}=read1D[sim1,"power_Str.dat"];
{t2,specgw2}=read1D[sim2,"power_GWs.dat"];
{t2mag,specmag2}=read1D[sim2,"power_mag.dat"];
{t2kin,speckin2}=read1D[sim2,"power_kin.dat"];
{tstr2,specstr2}=read1D[sim2,"power_Str.dat"];

Figure 11 (Time evolution fo the stress)

Clear[step]

plot2=Show[{ListLogLogPlot[{Table[{t2[[step]],(Sqrt[(((((t1[[step]])/6)))^2)((specstr1[[step,4]])/(4π *(0.3^2)))])},{step,1,25}],Table[{t2[[step]],(Sqrt[(((((t1[[step]])/6)))^2)((specstr1[[step,6]])/(4π *(0.5^2)))])},{step,1,25}]},Joined->True,PlotStyle->{Blue,{Blue,Dotted}},Frame->True,FrameStyle->Directive[Black,Bold,10],FrameLabel->{Style["t",Black,Bold,12],Style["\!\(\*OverscriptBox[\(\(\\ \)\(T\)\), \(~\)]\)(k,t)",Black,Bold,12]},PlotLegends->Placed[LineLegend[{Style["k=0.3",10],Style["k=0.5",10]},LegendMarkerSize->15],{0.2,0.78}],Epilog->{Text["(a)",{Log[6.3],Log[5.7*(10^(-6))]}]}],LogLogPlot[(Sqrt[((0.056*(10^(-9)))/(4π *(0.5^2)))(1/(((1+(((tt-0.2))^20)))^(-(1/20)*(5/12))))]),{tt,1,100},PlotStyle->{Red,Dashed},PlotLegends->Placed[LineLegend[{Style["fitted curve",10],Style["k=0.5",10]},LegendMarkerSize->15],{0.24,0.58}]]}];
plot1=Show[{ListLogLogPlot[{Table[{t2[[step]],(Sqrt[(((((t2[[step]])/6)))^2)((specstr2[[step,4]])/(4π *(0.3^2)))])},{step,1,25}],Table[{t2[[step]],(Sqrt[(((((t2[[step]])/6)))^2)((specstr2[[step,6]])/(4π *(0.5^2)))])},{step,1,25}]},Joined->True,PlotStyle->{Blue,{Blue,Dotted}},Frame->True,FrameStyle->Directive[Black,Bold,10],FrameLabel->{Style["t",Black,Bold,12],Style["\!\(\*OverscriptBox[\(\(\\ \)\(T\)\), \(~\)]\)(k,t)",Black,Bold,12]},PlotLegends->Placed[LineLegend[{Style["k=0.3",10],Style["k=0.5",10]},LegendMarkerSize->15],{0.2,0.48}],Epilog->{Text["(b)",{Log[6],Log[7*(10^(-6))]}]}](*,LogLogPlot[6*(10^(-9))(tt^(-2/6)),{tt,1,2}],LogLogPlot[12*(10^(-9))(tt^(-3/2)),{tt,2,7}]*),LogLogPlot[(Sqrt[20500*((10^(-14))/(4π *(0.5^2)))((tt^(-1/3))/(((1+(((tt-1))^10)))^((5/7)*(1/10))))]),{tt,1,100},PlotStyle->{Red,Dashed},PlotLegends->Placed[LineLegend[{Style["fitted curve",10],Style["k=0.5",10]},LegendMarkerSize->15],{0.24,0.28}]]}]
SetDirectory["/home/ramkishor/Desktop/GWspectra_from_mhd_turbulence (Version 1229)/mathematica/figures"];
Export["t_evolution_nonhelical.eps",plot1];
Export["t_evolution_helical.eps",plot2];