跳转到内容

脉冲星和中子星/脉冲星性质/Ppdot 图形脚本

来自维基教科书,开放的书籍,开放的世界

制作 p-pdot 图的脚本

[编辑 | 编辑源代码]

此图是用 PSRCAT 和 Julia 制作的。首先,要制作数据文件,请使用以下 PSRCAT 命令


psrcat -c "p0 p1" -o short -l "type(radio)  && exist(p0) && p1 > 0" -nohead -nonumber > all.dat
psrcat -c "p0 p1" -o short -l "type(radio)  && exist(p0) && p1 > 0 && pb > 0" -nohead -nonumber > binary.dat

然后使用以下 Julia 命令

#
using Winston;

all = readdlm("all.dat");
binary = readdlm("binary.dat");

all = sortrows(all,by=x->(x[1]));
binary = sortrows(binary,by=x->(x[1]));

# Use Chen, Ruderman 1993 death line

dline = 10^(2.0*78.0/7.0)*all[:,1].^(19.0/7.0)/3.2e19^2;
b10_14line = (1e14/3.2e19)^2./all[:,1];
b10_12line = (1e12/3.2e19)^2./all[:,1];
b10_10line = (1e10/3.2e19)^2./all[:,1];
b10_8line = (1e8/3.2e19)^2./all[:,1];

age_1000yr = all[:,1]/(1e3*365.25*86400.0*2);
age_100000yr = all[:,1]/(100e3*365.25*86400.0*2);

closefig()

# WARNING: color("red") is deprecated, use colorant"red" or parse(Colorant, "red")
# should be     add(p, Curve(x, c, color=colorant"red")), but this doesn't work
#
p = FramedPlot(xlabel="Period (s)",ylabel="Period derivative (s/s)",xlog=true,ylog=true,yrange=(1e-22,1e-9));
pts_all = Points(all[:,1],all[:,2],"symbolsize",0.5,"color","black");
pts_binary = Points(binary[:,1],binary[:,2],"color","red");
setattr(pts_binary,label="Binary");

pts_dline = Curve(all[:,1],dline[:,1],"type","dotted","color","green");
plot_b10_14line = Curve(all[:,1],b10_14line[:,1],"type","dashed","color","blue");
plot_b10_12line = Curve(all[:,1],b10_12line[:,1],"type","dashed","color","blue");
plot_b10_10line = Curve(all[:,1],b10_10line[:,1],"type","dashed","color","blue");
plot_b10_8line = Curve(all[:,1],b10_8line[:,1],"type","dashed","color","blue");

plot_age_1000yr = Curve(all[:,1],age_1000yr[:,1],"type","dashed");
plot_age_100000yr = Curve(all[:,1],age_100000yr[:,1],"type","dashed");

l = Legend(.1, .9, {pts_binary})
text1 = PlotLabel(0.1,0.63,"1 kyr");
text2 = PlotLabel(0.1,0.47,"1 Myr");
text3 = PlotLabel(0.65,0.9,"10^{14} G","color","blue");
text4 = PlotLabel(0.15,0.77,"10^{12} G","color","blue");
text5 = PlotLabel(0.85,0.24,"10^{10} G","color","blue");

add(p,pts_all,pts_binary,pts_dline,plot_b10_14line,plot_b10_12line,plot_b10_10line,plot_b10_8line,plot_age_1000yr,plot_age_100000yr,l,text1,text2,text3,text4,text5)

savefig(p,"ppdot.png")
华夏公益教科书