tempo2 -gr plk -f xxx.par xxx.tim -epoch centre

tempo2 -gr glitch -f xxx.par xxx.tim -t xxx.dat -fitf1 -offset (起始时间) -gt (glitch时间) -combine -font 1.2 -foot 0.1 -head 0.1 -title “(图像标题)” -p 6 -p 4

%需要手动调整的数据有：
% 导入数据名称
% 跃变时刻tglitch -----line 9
% *F1画图数量级slj ----line 36
% figure title -------line 69
%保存图片前，先全屏画图窗口。可保障图形纵横比一致性

%import data
tglitch=57964
data1=readtable('glitch1.txt');
data1=table2array(data1(:,1:4)); %前四列（第2，3，4列）是我们需要的数据
data2=readtable('g1.dat');
data2=table2array(data2(:,1:5));   %前五列数据。分别为toa(MJD),F0,F0误差,F1,F1误差
% 如果导入的文件含有字符，需要table to array把数值部分变为数组
%===========================================================================
% time residual 
toa1=data1(:,4);
residual=data1(:,2);
residual_error=data1(:,3);
l1=length(toa1) ; %residual数据的长度
idx1 = find(toa1>tglitch,1) ; %查找从第几个数开始，toa大于tglitch(这里不要改为>=)
idx1=idx1-1; %toa的前idx个数据都是在跃变前的
meanpreR=mean(residual(1:idx1)) ; %跃变前残差的平均值
residual=residual-meanpreR;  %计时残差减去跃变前残差的平均值，以使得数据贴近0，画图好看
%===========================================================================
% delatF0,F1 or possible F2
toa2=data2(:,1);
F0=data2(:,2);
F0_error=data2(:,3);
F1=data2(:,4);
F1_error=data2(:,5);
F0=F0*10^(9) ;  %因为画图时标注的单位是nHz
F0_error=F0_error*10^(9) ;
l2=length(toa2);  %delta v,即F0,F1数据的长度
slj=real(floor(log10(F1(1)))); %获取F1的数量级，取F1(1),用log10再向下取整(F1数量级是10的负数次方)，再只取实部
slj=-18   %手动调一调吧
F1=F1*10^(-slj); F1_error=F1_error*10^(-slj); %把F1数据改为数量级为0
%这里我有一个问题，只能先修改F1的数量级再X_F1修正改为贴近0。如果顺序反过来，F1误差的值会比F1的值还大，画图很难看
idx2=find(toa2>tglitch,1);
idx2=idx2-1;
mean_pre_F1=mean(F1(1:idx2));
X_F1=F1-mean_pre_F1; %使F1跃变前的数据尽可能贴近0，修改后的F1名为X_F1
[a1,E1]=polyfit(toa2(1:idx2),F0(1:idx2),1);  %一阶多项式拟合F0跃变前
[a2,E2]=polyfit(toa2(1:idx2),X_F1(1:idx2),1);  %拟合F1跃变前
[a3,E3]=polyfit(toa2((idx2+1):l2),X_F1((idx2+1):l2),1);  %拟合F1跃变后
%拟合函数ployfit返回值，a是数组，依次是各个阶的多项式系数值
% E是误差结构体，
deltaF0=F0-(a1(1)*toa2+a1(2)); %deltaF0，即画图的-p 6
mean_pre_deltaF0=mean(deltaF0(1:idx2)); % mean(pre_deltaF0),即deltaF0跃变前的均值
XdeltaF0=deltaF0-mean_pre_deltaF0 ; % XdeltaF0是deltaF0减去跃变前均值后的，使数据画图时贴近0，即画图的-p 5
XdeltaF0(idx2+1:l2)=deltaF0(idx2+1:l2)-mean(deltaF0(idx2+1:l2));
dF0=deltaF0 ; %开启此句，是选择画-p 6，没有调整的图
%dF0=XdeltaF0 ;  %开启此句，是选择画-p 5，跃变后数据减去跃变后均值的图
prefitF1=a2(1)*toa2+a2(2); %拟合的跃变前F1
postfitF1=a3(1)*toa2+a3(2); %拟合的跃变后F1
%===========================================================================
%figure together
Toa_min=min(toa1(1),toa2(1)); %判断两组数据的x轴谁更大，以更大的那个显示，强行把他们x轴拉成一样长的
Toa_max=max(toa1(l1),toa2(l2));
%---------------------------------------------------------------------------
figure(1) 
h1=subplot(3,1,1); % time residual
plot(toa1,residual,'.','MarkerSize',15,'color','blue')
axis([Toa_min-80,Toa_max+100,min(residual)-0.2*(max(residual)-min(residual)),max(residual)+0.2*(max(residual)-min(residual))]) %设定坐标显示区间
ylim1=get(gca,'Ylim'); % gca代表此时的绘图区，'Ylim'代表X轴的范围'
line([tglitch tglitch],ylim1,'linestyle','--','color','black','LineWidth',1.5);
set(gca,'xticklabel',get(gca,'xtick')) % 取消x轴科学计数法显示
xlabel('Modified Julian Date')
ylabel('Residual(ms)')
title("PSR J1731-4744")
%我没有画time residual的误差，否则应该把上边改为errorbar函数绘图
%---------------------------------------------------------------------------
h2=subplot(3,1,2); % deltaF0
plot(toa2,dF0,'.','MarkerSize',15,'color','black')
axis([Toa_min-80,Toa_max+100,min(dF0)-0.2*(max(dF0)-min(dF0)),max(dF0)+0.2*(max(dF0)-min(dF0))]) %设定坐标显示区间
ylim2=get(gca,'Ylim');
line([tglitch tglitch],ylim2,'linestyle','--','color','black','LineWidth',1.0);
set(gca,'xticklabel',get(gca,'xtick')) % 取消x轴科学计数法显示
% xlabel('Modified Julian Date')
ylabel("\Delta\nu (nHz)")
%---------------------------------------------------------------------------
h3=subplot(3,1,3); % F1 with error
errorbar(toa2,X_F1,F1_error,'o','color','black','linewidth',0.8)
axis([Toa_min-80,Toa_max+100,min(X_F1)-0.2*(max(X_F1)-min(X_F1)),max(X_F1)+0.3*(max(X_F1)-min(X_F1))]) %设定坐标显示区间
ylim3=get(gca,'Ylim');
line([tglitch tglitch],ylim3,'linestyle','--','color','black','LineWidth',1.0);
%set(gca,'yticklabel',get(gca,'ytick')) % 取消y轴科学计数法显示
%++----++
hold on %这5行是给F1画拟合线的，不想要就把这5行注释掉
plot(toa2(1:idx2),prefitF1(1:idx2),'-.','color','blue') %跃变前的F1拟合曲线
hold on
plot(toa2((idx2+1):l2),postfitF1((idx2+1):l2),'-.','color','red') %跃变后的F1拟合曲线
hold off
%++----++
set(gca,'xticklabel',get(gca,'xtick')) % 取消x轴科学计数法显示
xlabel('Modified Julian Date')
ylabel (sprintf ("\nv^{.}+%g (10^{%g} s^{-2})",-mean_pre_F1,slj)) % !!我不会表达\dot{\nu}
pos1 = get(h1,'Position'); % 获取第一个子图的位置 (每个子图的位置都是四元素向量，代表图左下角的坐标、宽和高
pos2 = get(h2,'Position'); % 获取第二个子图的位置 
pos3 = get(h3,'Position'); % 获取第三个子图的位置
set(h2,'Position',[pos1(1) pos1(2)-pos2(4) pos2(3) pos2(4)]) % 设置第二个子图的位置
set(h3,'Position',[pos1(1) pos1(2)-pos2(4)-pos3(4) pos3(3) pos3(4)]) % 设置第三个子图的位置
%我没有画F0的误差，否则应该把上边改为errorbar函数绘图