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

Reading the data files

{t1,specgw1}=read1D[sim1,"power_GWs.dat"];
{t2,specgw2}=read1D[sim2,"power_GWs.dat"];
{t1h,specgw1h}=read1D[sim1h,"power_GWs.dat"];
{t2h,specgw2h}=read1D[sim2h,"power_GWs.dat"];

n1=10;(*factor to get the minimum wavenumber*)
n2=10;
tstep1=1;
tstep2=1;

Figure 5 (Comparison of  GW power spectrum)

plot=Show[{ListLogLogPlot[{Table[{((i-1)/n1),specgw1[[1,i]]},{i,1,512}],Table[{((i-1)/n1),specgw1[[72,i]]},{i,1,512}],Table[{((i-1)/n2),specgw2[[1,i]]},{i,1,512}],Table[{((i-1)/n2),specgw2[[72,i]]},{i,1,512}]},Joined->True,PlotStyle-> {Black,Blue,{Dotted,Red,Thick},{Dashed,Gray,Thickness[0.003]}},Frame->True,FrameStyle->Directive[Black,Bold,10],FrameLabel->{Style["k",Black,Bold,12],Style["Sp((OverDot[h]))",Black,Bold,12]},PlotLegends->Placed[LineLegend[{Style["t=1.5",10],Style["t=37",10],Style["t=1.5, with modulus of stress",10],Style["t=37, with modulus of stress",10]},LegendFunction->"Frame",LegendMarkerSize->15],{0.42,0.28}],Epilog->{Text["(b)",{-2,-18}]}]}]

ploth=Show[{ListLogLogPlot[{Table[{((i-1)/n1),specgw1h[[1,i]]},{i,1,512}],Table[{((i-1)/n1),specgw1h[[72,i]]},{i,1,512}],Table[{((i-1)/n2),specgw2h[[1,i]]},{i,1,512}],Table[{((i-1)/n2),specgw2h[[72,i]]},{i,1,512}]},Joined->True,PlotStyle-> {Black,Blue,{Red,Dashed},{Dashed,Gray}},Frame->True,FrameStyle->Directive[Black,Bold,10],FrameLabel->{Style["k",Black,Bold,12],Style["Sp((OverDot[h]))",Black,Bold,12]},PlotLegends->Placed[LineLegend[{Style["t=1.5",10],Style["t=37",10],Style["t=1.5, with modulus of stress",10],Style["t=37, with modulus of stress",10]},LegendFunction->"Frame",LegendMarkerSize->15],{0.42,0.28}],Epilog->{Text["(a)",{-2,-18}]}]}]
SetDirectory["/home/ramkishor/Desktop/GWspectra_from_mhd_turbulence (Version 1229)/mathematica/figures"];
Export["helical_with_and_wo_phase.eps",ploth];
Export["nonhelical_with_and_wo_phase.eps",plot];