由于b站专栏文章排版的问题，许多链接和文字格式、文字颜色等都出现了问题，勉强能看吧

脉冲星周期跃变Pulsar Glitch

0.    开篇之前

[1]本文

本文作者(我)为XMU,CPST的wza(1972019-)

感谢给我莫大帮助的刘鹏学长

感谢素未谋面的周世奇学长

[2]说明：

本文作为一个操作手册，对脉冲星数据处理的完全新手做指导。较为详尽的阐述了从安装系统、安装软件、下载数据、处理数据、得到结果与具体参数、制作论文需要的数据图与表格的全过程。本文中：

绿色是文献的名字，没有给出链接的都可以在知网查到或轻易找到

橙色是代码或终端命令

深棕色是批处理文件内容(.bat)

蓝底：不得不看的信息、链接内容

浅蓝色：我的批注

        本文应该有1个可能需要的附带文件夹是：来自刘鹏师兄的压缩包，里边有使用源代码安装软件的方法、软件文件以及需要安装的库的内容的word

[3]推荐书目：

教科书：

向守平.《天体物理概论》.中国科学技术大学出版社

吴鑫基，乔国俊，徐仁新.《脉冲星物理》.北京大学出版社

《Handbook of Pulsar Astronomy》| Cambridge University Press |D. R. Lorimer,  M. Kramer,

《Pulsar Astronomy》 | Cambridge University Press. | Andrew Lyne and Francis Graham-Smith

相关阅读资料：

Radio Pulsars: Introductory Reading -- http://hosting.astro.cornell.edu/~shami/psrintro/

1. 在处理数据之前，需要一些理论基础       

1.1 Pulsar Star

        什么是脉冲星？脉冲星是中子星的天文对应体(counterpart)，建议看一看向守平老师的《天体物理概论》相关章节。我们需要知道脉冲星是快速旋转的中子星。自转周期0.1-10s，而有的，被称为“毫秒脉冲星”，具有1.3-30ms脉冲周期。（这些值有一些相关限制和计算，可以看吴鑫基，乔国俊，徐仁新老师的《脉冲星物理》的相关章节）

        脉冲星在快速自转，但根据灯塔模型，它的转轴和磁轴不重合。脉冲星的磁轴每旋转一周就扫过地球一次（我们能看到的脉冲星），因此我们观测到的脉冲星脉冲周期即为它的自转周期。

        对其辐射的理论为磁偶极理论，分为极冠模型，狭长间隙模型等。这些模型说，脉冲星发出辐射，自转动能提供辐射束能量，会给自己制动，所以自转频率会慢慢降低。根据磁偶极辐射模型，我们可以根据脉冲星的自转周期及其一阶导数求出来表面磁场和该脉冲星的特征年龄。

（以上公式来自刘鹏学长《FAST_19波束脉冲星漂移扫描巡天模拟》，可以看看此论文）

         但是这些模型不太好用，这些模型不能很好解释多波段的观测结果。

        对于脉冲星脉冲时序的特性的研究，我们采用称为“脉冲星计时”的方法（pulsar timing, 来自天文观测的O-C方法），主要研究脉冲星到达时间（time of arrivals, 即TOA）。

        望远镜观测得到一些列脉冲。对于脉冲轮廓，我们将尽可能多的脉冲轮廓进行相位对齐后进行叠加，得到一个标准脉冲轮廓。对于到达时间，由望远镜观测数据得到单个脉冲的到达时间，之后用尽可能多的脉冲到达时间进行加权平均，建立模板配置文件，用这个模板配置文件与所有到达时间做相关，以得到一系列TOA.（此段解释无法详尽说明，随着下文中的处理过程就会比较好理解了）

TOA需要考虑如下修正：

        主要就是观测站钟差校正，色散延迟校正，太阳系内的Roemer延迟，Einstein延迟，Shapiro延迟以及可能的（脉冲星的）双星系统的这三项修正。

        如果脉冲星辐射稳定，且不考虑可能的吸积过程时，观测数据将会严格按照如下公式描述：（此为一截断的泰勒展开）

1.2 周期跃变

        TOA受到两种不规则扰动影响：时序噪声（timing noise）和周期跃变（puslar glitch）。

时序噪声timing noise，论文《正常脉冲星计时噪声的研究进展_高旭东》有简明易懂的描述。

周期跃变pulsar glitch指脉冲星周期突然变化的情况。这个名词在以前的翻译是“周期突快”。但现在改为“周期跃变”，大抵是因为发现了反跃变和慢跃变，这些不够“突”也不变“快”的情况。

        在1969年船帆座脉冲星首次发现了周期跃变现象，如图：

        可以看到当时研究还是讨论自转周期，但我们现在大多讨论它的自转频率（在脉冲星数据处理软件中称为F0）和自转频率的一阶导数（称为F1）、二阶导数（称为F2）。

        周期跃变分为四类：正常跃变(normal glitch)，慢跃变(slow gliltch)，自旋延迟跃变(glitch with delayed spin-ups)，反跃变(anti-glitch)。并且有一些统计特征，我们需要知道周期跃变幅度()的分布有两个峰值：和左右。按此分为大跃变和小跃变。

        更多内容，可以读周世奇的文章，他的硕士毕业论文很好：《年轻脉冲星的周期跃变分析》以及他2022.12发表的综述《Pulsar Glitch: A Review》（https://www.mdpi.com/2218-1997/8/12/641），讲的很详细，值得一看。

        对于周期跃变的理论解释，有：[1]星震模型crustquake ，可以参看来小禹，徐仁新老师的论文：《固态夸克星星震模型的研究_云朝昂，来小禹》[2]超流，基于超流的几种模型。[3]中子星与外界交换动量、角动量模型，比如可能的吸积模型、小质量天体撞击，以及星风模型。但是这些模型都不是很管用，很多解释不了的东西，并且没有一个模型能预测周期跃变的发生。

        周期跃变研究的意义：我们做数据处理，以期得到统计规律；对脉冲星辐射发射特性的研究具有重要意义（现在磁偶极模型并不完善！）；周期跃变这些时序研究是研究中子星内部的探针（强相互作用力非微扰不可解，我们还不能得到很好的中子星星体EOS）；对PTA(Pulsar Timing Array脉冲星计时阵列)的研究具有意义，包括其对引力波的研究（可以看这个资讯这颗星星，献给他 :“中国天眼”FAST望远镜首次发现脉冲星！_百科TA说 (https://baike.baidu.com/tashuo/browse/content?id=66051d8862035cd29adc9ad5)）。它有个图很好：

2.  处理脉冲星观测数据首先需要有数据

        有很多射电望远镜会观测脉冲星，比如英国Jodrell Bank observatory（JBO）天文台Lovell 76m , 澳大利亚ATNF的 parkes 64m，以及中国FAST 500m。新疆天文台25m望远镜。但是很多望远镜的观测数据并不开放，我们只得处理Parkes的数据。

2.1 下载观测数据

        澳大利亚ATNF天文台，在网站CSIRO Data Access Portal - ATNF Pulsar Observation Search(https://data.csiro.au/domain/atnf)处可以下载pulsar数据，包括原始数据.rf和预处理数据.FTp以及.f8Tp 见如下示例：

在source name处搜索星名，如J1731-4744（这个源信号很强！）

Band name选20cm波段（我们主要处理观测中心频率为1369MHz波段的数据）

点最下边的search。在结果页，Options处可以进一步筛选：

中心频率选1369MHz(其实是1368.75MHz)，源选J1731-4744，带_R的源是校准源，不需要。Obersvation type选预处理的数据，为.FTp或.f8Tp（没有预处理的raw为.rf数据）

和

及

其他不管

文件名字前带感叹号的，是暂时无下载许可的文件（ATNF有18个月数据保护期，或是其他不可公开访问项目），这些文件不能直接下载。（可以先点击下载，填写邮箱。数据库会稍后给你的邮箱发送邮件，包含许可证和下载链接，从该链接下载。但是数据多了会有一大堆邮件很麻烦，或许可以写一个自动抓取的代码*）

下载数据，可以排序查找，比如这个按日期排序

选择数据，在数据前边的框里打对号，或者页面全选

下载数据，可以逐个下载（直接下载为FTp）或全选后download，选择TAR或ZIP。

（如果选择的数据超过8GB，需要另一种下载方式WebDAV,但是预处理过后的数据很小，平均只有几十kb，不会有这么大的数据出现）

有问题也可以点击页面上方的help来查看帮助文档

        下载好的zip或tar压缩包，解压之后会有非常多的文件夹，文件在这些文件夹里，可能还包括名称为CSIRO（Parks所在机构）的许可证之类的文件，这个没用。

        这时候可以用批处理文件提取，比如：新建文本文档，添加以下内容for /r "C:\source_folder" %%a in (*.FTp) do copy "%%a" "C:\destination_folder"其中，C:\source_folder是包含要提取的文件的文件夹，*.FTp是要提取的文件类型，观测数据文件格式为.FTp以及可能有的.f8Tp两种。C:\destination_folder是要将提取的文件复制到的目标文件夹。保存文件，改尾缀为.bat（在windows下），运行，可以把观测数据统一提取到一个文件夹里）。

        在linux下，tar压缩包，在终端中对.tar文件解压缩，命令为tar -xvf （文件名）.tar ，打开预处理的文件夹，看到预处理过的FTp文件。

        文件名称一般为，eg: t160133_154435.FTp，前半部分名字意思为16年01月33日，154435或许是那天观测的具体时间（在观测时长之内）

2.2 Glitch数据库，星表Psrcat及其他信息

        [1]对于glitch事件，JBO和ATNF两个天文台分别有数据库记录。数据库 JBO :

 (www.jb.man.ac.uk/~pulsar/glitches/gTable.html）和ATNF（The ATNF Pulsar Catalogue | Glitch Database (https://www.atnf.csiro.au/people/pulsar/psrcat/glitchTbl.html)）可以查看已经发现的所有glitch。两个数据库的数据互为补充。

我一般查看ATNF数据库：

括号中为误差。Q为恢复系数，T_d为恢复时间(?)。如果这两项有值，则说明该次glitch有跃变后（指数）恢复过程。而JBO数据库有个好处，点击参考文献可以直接打开原文连接。

        [2]在ATNF的网站The ATNF Pulsar Catalogue (https://www.atnf.csiro.au/research/pulsar/psrcat/),即psrcat的网页以查看已知源的更详细信息。（这个视频简单叙述了怎么在这个网站查找源的详细信息，包括利用网页功能作图，链接：https://www.youtube.com/watch?v=rjMljTpbygY&list=PLq68qn0GOZVbkSHIqoA5-LAX8LpbMTQV7&index=4&t=3s）这个网站功能很强大，可以列表显示脉冲星信息。选择参数并打勾，点击‘table‘会得到脉冲星的参数信息表。同时还有plot画图显示的功能。可以从这里筛选选择脉冲星。

        举例，我选了一些参数，再点击上方“TABLE”按钮，可以查看数据库结果表格。

参数很多，具体每个参数是什么，可以查看它的说明文档：

我们需要知道，怎么从psrcat网页上得到.par文件，此为脉冲星自传模型文件。在该网页下面框中输入脉冲星名，点击“get ephemeris”，如图：

即可得到：

我已经鼠标选中了数据部分，复制，新建一个.txt文件，黏贴进去，改名为J1731.par（文件尾缀改为.par）

一项补充内容

儒略日： MJD

MJD 0对应于1858年11月17日午夜。2014年1月1日对应56658 MJD。不足一天的按小数计算

我找到了一个网站用来换算Modified Julian Day Converter (csgnetwork.com)

3. 数据处理得有数据处理软件

3.1安装软件之前得先有合适的系统

这些软件需要linux系统，Mac似乎也行(只能在类Unix下运行?)。那么如果你的电脑是windows系统是不行的。我推荐你安装ubuntu18.04（随着ubuntu系统更新，也可以选择20.04或22.04，但一般不要选择最新版系统）

那么你有如下选择：

[1]双系统

[2]虚拟机

[3]WSL

实际上后两者类似，也可以使用远程服务器等。WSL应该是会比虚拟机更好。

我推荐使用双系统，只要你系统安的好，就问题少。

这里注意，WSL和虚拟机，没有自己的图形界面，在作图时（即pgplot）需要调用主机的图形界面，此时是使用Xserver，要xhost +权限。而据刘鹏师兄说xhost +不可以在conda下运行。应该是要比较麻烦的设置，使用Xming等(?)，比较麻烦。

[1]双系统安装ubuntu18.04

首先下载系统镜像，ubuntu官方网站下载 https://releases.ubuntu.com/18.04/ 如图：

然后正如安装window系统一样，安装linux你也需要一个启动U盘，制作安装镜像，现在有很多工具做这个事情，比如rufus，UltraISO，以及最近的Ventoy，我只用过rufus，挺不错的我推荐这个。但也听说Ventoy很高级很好用。

具体安装流程可以在CSDN或知乎上搜一搜，教程大致都差不多，可以看一下这些参考：

(1)    Windows10安装或重装ubuntu18.04双系统教程（平民教程）_win刷机ubuntu_四处炼丹的博客-CSDN博客

(2)    Windows10安装ubuntu18.04双系统教程 - 不妨不妨，来日方长 - 博客园 (https://www.cnblogs.com/masbay/p/11627727.html)

但是注意，听我一句劝，linux系统设置分区的时候最好按照默认分区！特别是第一次装linux的朋友。

总体步骤大致为：

|| 在官方网站上下载ubuntu18.04镜像 ，这些镜像按版本不同大小在2-4G

|| 下载rufus(或者其他镜像刻录软件 ,但我更推荐rufus)以制作启动引导U盘。注意启动引导项是否为UEFI ； （此时在windows）电脑硬盘管理中划分新分区(留给linux用)

|| 插入U盘，重启电脑，进入bios，安装系统。注意，在这教程大多教你怎么自己分区。如果你电脑分区比较利索的话，我建议直接按默认分区，ubuntu和windows共存那个选项就好，省事，并且做的越少错的越少

|| 按照提示安装。安装好了。

注意安装好了后，系统会提示你更新，可以软件更新，但不要更新系统版本！！

[2]虚拟机

使用VM安装虚拟机。VM下载网址 https://www.vmware.com/cn/products/workstation-pro/workstation-pro-evaluation.html

然后百度搜一搜，找一个注册码破解使用。安装方法详见链接https://blog.csdn.net/weixin_45912291/article/details/108901106

        为了便于ubuntu系统使用，可以安装VM tools, VM软件会自动安装或者提示。换源意义不大，如果真要换，记得备份默认源。虚拟机不能按语言安装，这时系统是英文的。按提示升级软件包打开设置，设置language中添加中文，应用到全局；这里可能需要更新要等一会。可能需要重启，之后就成为中文系统了。

[3]WSL

现在有WSL，即windows subsystem for linux

听说功能强大，操作方便

Tips

Linux常用命令：
有一推送可以看看Ubuntu：你必须知道的常用命令 (https://mp.weixin.qq.com/s?__biz=MjM5NjY1NjYyMQ==&mid=2650152742&idx=4&sn=88455ea528cb7b339e063fd362a41ad8&chksm=bee75ba48990d2b26daf355e1e7991031498c6aea93b1fac2a696e4d4b7c591c843141973e5d&mpshare=1&scene=23&srcid=0416557bvZlhcCFQO2HhQ4F4&sharer_sharetime=1681718480687&sharer_shareid=7b2d09c7f095de1bf7463b817742e58f#rd)

sudo apt-get install (software)  以apt安装软件，sudo为加权限

sudo passwd root （设置root 密码）

su root （打开root）

ls 查看文件目录

ls -a 查看隐藏文件目录

cd  ./(file） 打开文件目录 

cd ..  返回上一级目录  

sudo gedit (documents.eg:  .bashrc/)  利用自带编译器gedit编辑文件

source ~/.bashrc   执行对 .bashrc 的修改‘

sudo apt-get update 更新software

sudo apt-get upgrade 升级软件包

touch new.txt 新建.txt文件

修改的.bashrc/内容是指向软件安装目录。这是系统环境变量

注意bashrc里，export path那段，等号左右不许有空格

Linux系统中，修改文件后缀并不影响文件格式。

可以下载安装QQ，baidu搜索qq for linux， 选择x86 deb那个，下载后点击安装。另外，建议安装腾讯会议。

3.2安装软件

你可以采取如下这些方法。如果顺利，使用[1]就能完成所有安装了，这会非常棒。

[1]conda

写在3.2-[1]节前边：我专门写了这一部分怎么操作的文档，看这个：

脉冲星软件的一种安装方法Pulsar software install :tempo2, psrchive, etc. - 哔哩哔哩 (bilibili.com)

先安装anaconda，再添加conda-forge这个channel(为conda包管理器启用conda-forge通道)

conda-forge是一个好用的channel，允许我们直接安装脉冲星软件。

之后你可以为pulsar software创建环境，并直接以conda install -c conda-forge (software name)安装。

[2]kernsuite

网址为kernsuite.info

这是通过添加ppa的方法，允许我们直接用apt-get install脉冲星软件。

按照首页的那个第一个命令框提示把那几条命令在终端运行

注意，把第二行kern-8改为kern-5，因为我用的是ubuntu18.04，你打开网站顶部的packages看kern-5里，查看对应ubuntu版本(应该可以向下兼容)和对应的软件

之后就可以sudo apt-get install psrcat 等等软件

同时可能需要安装一些依赖和库，参看刘鹏师兄那个word里边那些dependence

这种方法仍然需要添加tempo2等软件到系统目录bashrc

[3]糟糕的古早方法

如非必要，请不要采取此法

安装软件，主要就是psrchive和tempo2

主要有以下几步：

(i)安装psrchive, tempo2, 编译运行安装；

(ii)安装他们需要的依赖，一些其他的包和库，比如(python)anaconda, pgplot5, fftw, cfitsio… (iii)把这两个软件（和某些所依赖的软件）添加到系统路径里，使得软件运行的时候能找到自己需要的文件（~/.bashrc）

(1)   你可能需要使用python，可以直接安装anaconda,比如：

https://zhuanlan.zhihu.com/p/459607806 注意文件名和路径，最后那个bashrc需要手动添加到路径，然后在终端中输入conda –version ; python –version可检验是否已安装

（2）对于tempo2, 推荐按照这个教程安装：

https://www.bilibili.com/read/cv12011917?spm_id_from=333.999.0.0

（这里也说了安装cfistio, gsl , pgplot5, fftw3 fftw3-dev pkg-config这些，他是直接下的，但是我感觉pgplot5这种还是需要添加kernsuite.info这个ppa的吧）

值得一提的是，这个git clone下来，git慢的很，慢得都连接超时下载失败了。也可以直接在tempo2的官方网站里download界面下.tar.gz再自己解压，指路：https://bitbucket.org/psrsoft/tempo2/src/master/）

其中，./configure --prefix=/your path ,这里=后边是编译到的路径，不写prefix的话就默认路径安装（可能把可执行文件放在/share里，把库放在/lib里等等，乱）

对于我，我是./configure --prefix=/usr/share/tempo2，这样后边添加path也就清晰：

export TEMPO2=/usr/share/tempo2/T2runtime

alias tempo2=/usr/share/tempo2/bin/tempo2 (我不知道这一行有没有用啊，好像注释掉也可以正常用tempo2。Tempo2官方教程里没有alias这行，但是这个B站教程有)

在安装TEMPO2之前需要安装pgplot、FFTW、CFITSIO等数学库，这些东西有问题的话，这一步会报错。

（3）安装其他一些库和依赖，参考刘鹏学长的word里边的依赖目录

一些软件安装资讯

下面的信息不是必看的,且由于b站专栏的原因，这些链接全部失效了

软件工具 | 国家天文科学数据中心 | NADC (china-vo.org)

如果你已经正常安装，那么下面的安装教程请忽略

脉冲星计时软件tempo2的安装 - 哔哩哔哩 (bilibili.com)

下载并处理脉冲星射电数据 (qq.com)

安装pulsar sfotware：https://mp.weixin.qq.com/s/i5EFJimKFGvoP8rsa_GCQQ

安装ubuntu链接：https://mp.weixin.qq.com/s/Qmuhd5AnuKM7rPvy6BgeEg

脉冲星射电数据处理入门 (qq.com)

Introduction to Pulsar Timing (qq.com) (处理原理介绍)

PTA软件的安装的一些坑 (qq.com) （这里他用的WSL安装ubuntu，需要用Xming处理Xserver调用主机图像界面，但他也没说的很详细）

 

Tips

[1]Ubuntu系统截图：

Alt+PrtSc 截图当前窗口，图片保存在 目录/图像

Shift+PrtSc 框选区域截图，图片保存在 目录/图像

Ctrl+Shift+PrtSc框选区域截图，保存在剪贴板

[2]没太大用的tips

对于数据下载和挑选常用的四个网站，可以用批处理文件快速打开，如下

在记事本中输入：

@echo off

start /min msedge.exe "https://www.atnf.csiro.au/research/pulsar/psrcat/"

ping n 0.25 127.0.0.1>>nul

start /min msedge.exe "https://www.atnf.csiro.au/people/pulsar/psrcat/glitchTbl.html"

ping n 0.25 127.0.0.1>>nul

start /min msedge.exe "http://www.jb.man.ac.uk/~pulsar/glitches/gTable.html"

ping n 0.25 127.0.0.1>>nul

start /min msedge.exe "https://data.csiro.au/domain/atnf/"

ping n 0.25 127.0.0.1>>nul

#taskkill /f /im msedge.exe

exit

保存，更改文件后缀为.bat，要打开网站时双击即可(在windows系统下)

4.  数据处理

我们已经在2.1节中从ATNF下载了观测数据（以J1731-4744为例，我们下载了其从2007-2021年间可获取的观测数据，共计128个预处理过的.FTp文件，带宽256MHz，中心频率1369Hz）

4.1数据筛选与得到计时残差图

        打开数据文件所在的文件夹，打开终端，执行sudo chmod 777* (注意*前边有空格。*是通配符，意思是所有文件。777是linux下文件权限的某种代码)，获取对文件的访问编辑权限（在命令行中该文件夹下输入ls，文件名变为绿色就是有权限）（有些文件目录下直接就有操作权限，不用此命令）

        可以创建三个文件来分门别类放置数据来方便查看，比如‘goodSN’ , ‘fit’ , ‘notgoodprofile’ ,第一个文件夹放置稍后得到的轮廓最好的模板文件，第二个文件夹放置处理得到的自传参数等文件，第三个文件夹放置轮廓不好（意味着观测不好和观测流量低）的数据。创建文件目录可以使用命令，如 mkdir goodSN

        本节主要讲述了对数据进行筛选，得到由脉冲星观测数据得出的观测计时模型和ATNF星表得出的自转模型，并以此得到脉冲残差图的过程，流程如下：

4.1.1轮廓筛选

        (1)首先我们需要查看数据文件bin的数量，bin是数据图的横坐标数(?)。

使用命令vap -c nbin *.FTp 来查看文件bin的数量（其他类型文件的话就把.FTp改成其他文件后缀），如果不一致，比如有的是1024，有的是512，那就需要改。修改bin数只能改小不能改大（只能降采样，不能进行虚假的插值），使用命令：pam –setbin 512 *.FTp -m

        (2)之后我们需要查看脉冲轮廓以进行筛选，看轮廓形状是否噪声太多。查看脉冲轮廓的图像，可以：

[1]使用psrchive的pav命令，在终端中输入pav -CD *.FTp（pav是绘图命令，CD是centre display, *.FTp指显示所有.FTp文件。如果有.f8Tp文件就再pav -CD *.f8Tp显示。可以看到如下图像：

        在终端中按‘Enter‘键转到下一张图。用👁👁眼观察数据，筛选掉坏的数据。选择标准：S/N信噪比大于5为适宜，小于5噪声太高，考虑不采用。或者一些明显噪声高，脉冲难以分辨的数据，也不采用

(如果要单独查看某一图像，如a070722_091632.FTp的轮廓，使用命令pav -CD a070722_091632.FTp）

[2]使用psrchive的psrplot命令，把.FTp批量转化为图片文件以进行查看

        使用命令psrplot -p flux *.FTp ，终端提示graphic device/typr(? to see list,default /NULL): 这是提示输入文件名和文件类型以供它生成图片文件。如果这时输入 ? 它会提示可以选的图形文件, 有ps或png等。可以输入all.ps/ps （斜杠前是包括尾缀的文件名，斜杠后是文件类型）可以得到一个包含所有数据图信息的.ps文件（在linux下非常方便查看和筛选轮廓）。或者得到.png等文件以方便查看。但显示的图像只有FTp文件名，没有信噪比数据。查看观测数据信噪比的方法见下文。

        不采用的数据，可以直接手动移动到notgoodprofile文件夹里，也可以用命令rm filename notgoodprofile （没有权限的时候需要使用命令sudo rm filename notgoodprofile）

4.1.2 FTp文件相关信息查看

4.1.2节不是操作流程，而是tips说明。

Psrchive有很多的工具包，这些插件实现了很多功能，详细内容可查看该软件文档https://psrchive.sourceforge.net/manuals/，一些工具如下

psredit filename 查看该文件信息，包括中心频率，观测历元，观测站信息等，以及psrstat工具

psrplot -p flux *.FTp 然后输入all.ps/ps 可以得到一个包含所有数据图信息的ps文件（非常方便查看和筛选轮廓）

psrplot -P可以查看psrplot支持的图形类别，其中flux是流量图，适合FTp

psrtxt filename.FTp > txt_filename.txt 可以把FTp文件导出为txt文件，数据为横坐标与流量，支持批量导出，如psrtxt *.FTp > J1731-4744.txt，将所有FTp导出为一系列txt文件

可以psrtxt -h 以查看更多用法。可以自己导出txt后，使用matlab或python以进行更多操作，如绘图等。

pas -r reference_file filename 以reference_file(.std或.FTp)为模板，进行相位对齐。filename可以为*.FTp

psradd -FTP *.FTp -o new.FTp -v轮廓叠加，把数个pulsar profile叠加在一起，以得到高信噪比、脉冲轮廓平滑的数据，一般是把glitch前后的数据以此叠加，得到轮廓形状。

4.1.3 得到计时模型.tim

接4.1.1节

（1）选择最好的脉冲轮廓作为标准模板

        可以通过如下命令查看FTp文件信噪比pdv -f *.FTp | sort -k 10 -g 终端显示的最后一行是是信噪比最高的FTp（sort是排序）。在终端输出的结果中，S(mJy)是望远镜接受流量，单位毫央(斯基)；W10是平均脉冲峰值强度10%的地方的宽度，W50是平均脉冲峰值强度50%的地方的宽度（半高全宽）(详见《脉冲星物理》3.1节)；最后一列S/N即为信噪比。初处理中一般选择信噪比最高的轮廓(.FTp)作为标准脉冲轮廓（如果该轮廓信噪比最高但轮廓形状很差，也可依信噪比大小选其他的）

        可以将该文件复制到goodSN文件夹中以方便查看，直接复制或cp filename goodSN

（2）计时模型获取

        得到标准模板.FTp后，将其复制为模板文件：cp filename std_filename.std (比如处理J1731-4744数据时，我将其命名为J1731.std)

        使用psrchive的pat，命令为pat -s J1731.std -f “tempo2 -i” *.FTp > tim_filename.tim得到到达时间文件.tim（如J1731-4744我命名为J1731.tim）（命令中，大于号>意思是输出为；以tempo2的文件格式： -f``tempo2 -i`` ; 用*.FTp文件与1731.std文件做相关，注意各处空格，注意”是英文双引号！）

        可以用emacs来查看、编辑并排序该.tim文件。也可以简单地用ubuntu自带的gedit编辑器打开gedit filename.tim（提示没有权限时，在命令前加sudo）。也可以用cat命令打开.tim文件，还可以使用vim编辑，在命令行中vim filename.tim即可打开。输入：冒号之后可以输入命令，比如:q是退出，而:wq是保存退出

.tim文件内容如上。第一列是文件名称，第二列是中心频率，第三列是到达时间（MJD）,第四列是误差(什么的误差?），pks是parkes望远镜简写。 -i是后端存放数据的模组。 PDFB4是最新的数据存放模组。

        gedit没有排序功能，需要在终端中使用ubuntu的sort命令，如sort -n -k 3 filename.tim是对该文件第三(3)列(即到达时间）(-k）按数值大小(-n)排序。 注意这是把结果打印在终端里的，得到新文件需要改为命令：sort -n -k 3 filename.tim > newfile.tim 输出为新的tim文件。我们只需要排序后的新.tim (注意，tim文件的注释符是一个单独的字母C ）

4.1.4得到自转模型.par

        使用psrcat软件可以直接得到自传模型.par，其中包括了对该脉冲星已有的自转模型的参数，如脉冲星名(PSR)、脉冲星赤经(RA)和赤纬(DEC)、自转频率(F0)、自转频率的一阶导数(F1)，以及得到这些参数的历元(Epoch)和其他信息。稍后的4.1.5节中讲述了提供这些信息给tempo2，让其根据.tim的到达时间建立计时模型。

命令 psrcat -e J1357-6429 这将直接在终端里显示该脉冲星的参数

命令 psrcat -e J1357-6429 > 1357.par 导出为.par文件

        现在通过conda-forge安装脉冲星软件还不支持psrcat，可以使用其的网页版：The ATNF Pulsar Catalogue (csiro.au)  网页下拉，在这里输入源的名字：

点get ephemeris可以看到信息。我们要得到.par文件，我是这样的：新建一个.txt，然后把网站上得到的那段数据复制进去，再改尾缀为.par

以下是.par的文件内容：J1731-4744的例子：

天球坐标：RAJ是赤经，DECJ是赤纬

DM指色散量，射电信号在星际传播中造成的时间延迟（cm-3 pc）

PEPOCH是参考时刻，即脉冲星周期和周期一阶导数的参考时刻（MJD）

RM 自旋测量(rad m-2)

DMEPOCH 星际传播的色散延迟

PMRA 赤经运动（mas/yr）PMDEC 赤纬运动(mas/yr)

EPHVER 脉冲星时空参考系的版本号（这个更新了，需要把2改成5）

UNITS 时间单位； TDB 即Barycentric Dynamical Time（这个更新了，需要把TDB改成TCB）

还有个这个.par文件没有的EPHEM DE405是什么太阳星历（这也更新了，更新到436了，最少改为430）

F0 质心自转频率（Hz）

F1 ：F0对时间的一阶导(s-2)

F2 ：F0对时间的二阶导(s-3)

P0 ：质心自转周期（s）

P1 ：P0对时间的一阶导（无量纲数）

4.1.5得到计时残差图

        拟合(详见《脉冲星物理》4.2.4节，p91),psrchive的pat包采用的方法原理为把脉冲轮廓与标准轮廓做相关，整体考虑脉冲的轮廓，得到拟合参数(即计时残差)及其误差。

        其中有些参数我们要拟合，有些不拟合。这时需要打开.par文件查看。我们直接在第二行数据后边添加0或1数字，数字1代表拟合这个数据，0代表不拟合。（另外，从tempo2图里保存new par时，保存的文件START和FINISH后边的数字可能是2，需要改成1.处理含有多个望远镜的数据时会报错）

之后使用命令：tempo2 -gr plk -f 之前得到的par（如：1357.par）之前得到的tim（如：J1731.tim）-epoch center

其中： gr是图形界面，残差分析的GUI ； plk残差分析插件的作图，调用这个插件； -f 是打印information ；-epoch centre是以数据中心为图像中心显示

这里注意，它可能会报几个warning，但不是error，没关系

        如果终端中一直报这个：

这是由于，它找不到图形界面输出pgplot的图，需要xhost+程序调用图形界面。如果是在本机安装的话是不用的，如果是虚拟机、VNC、远程桌面这些就需要调用本机的图形界面软件。 另外，刘鹏学长说conda环境下无法使用xhost +

        终端输入命令 xhost +（注意+前边有个空格）。并确认一下pgplot和psrchive都安装正确。这里虚拟机总是会出各种问题，需要在主机上安装Xming这种Xserver,然后再配置，WSL应该同理。总之我没搞明白，我最终用的还是双系统。

4.2拟合计时残差

4.2.1找到周期跃变(1)

        执行tempo2 -gr plk -f xxx.par xxx.tim -epoch centre命令后会出图，按J以线条连接数据点，如下：

        可从头开始，或在图像中选择线段部分进行拟合。理论上对拟合模型，我们可以先拟合2个，再拟合4个，再6个，再8个。两个两个地拟合。如果8个以上都拟合得很好，那后边的都按照这个拟合就可以了。对于复杂的残差图，我们做的时候，先拟合三个，再一个一个地往上加点。

        点击键盘” z “(放大，zoom），再用鼠标左键框选想要放大的部分。点击键盘“ U ”键，会转回显示所有数据点的界面。

        对着数据点点击鼠标中键，会出现一个小圆圈，在另一个图形窗口会显示这个数据点的.FTp文件的图像以供快捷查看（仅当该数据.FTp文件和.tim、.par文件在同一个文件夹里时才可显示）

        对着数据点点击鼠标左键，会在命令行里显示该数据点的详细信息，比如文件名、MJD、到达时间误差等

        按’ x ‘键，对时间残差进行拟合。此操作是调整timing model的相位，以得到残差最小的到达时间模型。多次点击进行拟合，使得‘rms’值降低至小值(注意，选的点少的时候，按’ x ‘拟合，可能把rms弄成Nan，这时按’ U ‘回去，接着按x，它还可以继续拟合)

        按键盘“Q”键，即可退出该图形窗口。

        对于图中明显不符合趋势的数据点，考虑可能是信噪比太低等原因导致拟合异常，可以鼠标右键删除该点。如要彻底删除，需要在.tim文件中注释掉该点。.tim文件的注释符是大写字母C，可以注释该行。

-----先选第一段直线结构(按z，框选，图像会放大到如下所示：)

反复按x键拟合，直至rms降到小值（如降到300或500以下），如图所示：

按U回到全部数据点视图，看到第一段曲线结构：

此时点击F2，可以拟合F2。选择第一段曲线结构：(注意只选曲线的点，不要选偏折后的点)

反复按X拟合，等到残差和rms又变好了：

这时按U回去看，发现第一段已经被修正的成近似直线：

这里发现前后两段直线中发生了偏折，这是明显的小跃变的现象。直线最后一组数据和折线第一组数据：s110618_124943.FTp与s110717_105323.FTp之间发生了第一次周期跃变，是小跃变。

4.2.2跃变数据保存

在找到第一处周期跃变后，我们需要保存周期跃变的数据。

跃变前后的两个数据点，即图中偏折前后的两个数据点，点击左键，在命令行中会出现他们的具体信息，如MJD和文件名。记下这些信息。记录表格如下：

        其中，glitch-epoch即为跃变发生的历元，取跃变前后两个数据点到达时间的中间值。这里的跃变历元误差为前后两个数据点到达时间间隔的一半。

        选中跃变前的直线，多按几次X，使得rms尽量小。然后选中该段，点击图形窗口中的“new par”按钮（或按shift+P），在命令行里会提示输入文件名称，可以输入”pre1.par”,即可保存第一段跃变前的参数为新的pre1.par文件。此为框选部分拟合得到的自转参数，数据对应的epoch是框选部分的时间中点。下次再处理这一段时就无须从头开始了。拟合跃变后的直线的数据，框选跃变后数据，像前一段一样重复按X拟合，直到跃变后的数据也拟合为水平直线，如图：

框选跃变后部分，拟合，shift+P保存为post1.par。跃变后的post.par使用尽可能多的数据点。

        （也可以在图形窗口中把除了跃变前后以外的数据点都删掉（挨个按鼠标右键），直到只剩跃变前后的数据点这是小跃变，如果可以，前后各取20个点，但注意不要包括其他跃变的数据点），然后点击图形窗口的”new tim”,在终端中输入文件名，如glitch1.tim，这样下次作图时使用这个tim就不会被其他无关数据点干扰。但我无此习惯）（对于保存.tim，也可以复制之前的.tim文件，然后把无关数据点删除。注意保留.tim文件的表头FORMAT 1）

        框选跃变前后的数据点，按ctrl+j保存为glitch1.txt文件，该文件保存的内容为这些跃变前后的数据点的到达时间toa和计时残差residual以及误error差，用于稍后画timing residual残差图。

        那么现在我们已经找到了第一个glitch，记录了glitch-epoch并保存了pre1.par和glitch1.txt

4.2.2找到周期跃变(2)

继续观察第一次跃变后的图形，看到一条一条平行的直线，这些折现实际上是一条直线，因为相位跳变才出现了这种结构，可以手动消除相位跳变，也可以在.par文件中添加一行，内容为TRACK  50，意思为把数据点前后50天的数据捆绑（改成TRACK 100，那就是100天的），以此消除相位跳变。

        手动消除相位跳变，对于相位跳变的点，左键点击该数据点左下角，图中的点相位升高，故再点击键盘“-“减号键（若相位降低，则点击”+“加号键），会出现如左侧所示黄色虚线，折线会变为直线。如下图为示意：

J1731-4744的数据很好，只需要选跃变后的第一段折线进行拟合即可得到结果。

之后框选该段直线，只选F1，反复点击“X”拟合直至为直线。若出现曲线结构则加上F2拟合曲线结构。应该会得到如下图像：

记录下该次跃变前后的两个数据点的数据，shift+P保存第二次跃变前的参数pre2.par，跃变后的参数post2.par，第二次跃变的残差和toa数据glitch2.txt。该段数据点恰好在跃变时缺失，故此次跃变epoch的误差很大。

        对于后边的数据依次按此处理，或许可以找到更多的跃变事件。

4.2.3找到周期跃变(3)

以上的举例都是小跃变，接下来我们看看大跃变是怎么处理的。对于源J1731-4744，在MJD~57978附近发生了一次大跃变。按上边几节的方法可以处理得到大跃变前的计时模型，如下：

放大第四处周期跃变，可以看到大跃变的特征：

跃变后的计时残差呈散点状，此时无论怎么拟合跃变后的数据点，都得不到直线的拟合结果。这是由于大跃变幅度很大，跃变前的计时模型(.par)无法拟合跃变后的情况。按照前边几节所说的保存pre4.par，glitch4.txt和跃变前后数据点的epoch。更多数据处理见下面章节。

4.2.4注意事项与总结

[1]跃变的判定

大跃变后的数据呈散点状。小跃变要求跃变前后的数据点是明显的直线偏折。下图这种带弧线的情况不是小跃变：

对于我们采取的Parkes数据来说，它一般间隔两周观测同一个源，间隔太长，不容易看到小跃变。而在这方面新疆南山的优点就体现出来了。

[2] 对于缺失数据，像这种：

可能的跃变后的关键数据缺失（青色圈里部分），所以不能确定是不是发生了glitch。

想进一步确定的话，可以去ATNF上看看这段时间内的rf文件可不可以补充这段观测缺失部分, 要是有的话可以自己去对rf文件进行消干扰，然后生成FTp文件，补充这部分。

[3]标识glitch-epoch

有个小窍门：对于ATNF, JBO已经记录的glitch，可以在.par文件中添加

例如源J1357-6429的在我观测数据范围内的几次跃变，添加以下内容到.par里：

GLEP_ 55576

GLEP_2 57795

GLEP_3 58148

在tempo2画残差图时会在图像里在对应时刻显示紫色虚线，方便识别。

[4]不同频率的观测数据

有的源有不同中心频率的观测数据，在tempo2的图中会以不同颜色来表示。按‘ J ‘连线，它只在不同频率观测值之间连线。但我们简单处理数据不考虑频率是否相同。我们需要消除不同频率观测数据间的相位差异，在.par文件里加jump内容（在tempo2残差图里用“+”或者“-”消除相位后，保存新的.par文件，会自动在新.par文件里添加该相位跳变信息）

[5]关于拟合数据保存

shift+P保存数据时注意, 先框选数据点，再拟合，再保存。不要拟合后又重新框选，会出问题。

[6].tim文件的异常误差

有的.tim文件里边，脉冲到达时间的误差是0.000，这不行，手动给他改成，比如0.001

[7]更好的数据点数据

对于下图：青色为脉冲轮廓。红色是其他干扰。

想要更好的处理，就自己去下载.rf(是raw)，没有预处理的文件，把那些尖锐的消掉

我们这是粗略处理，那些不好的轮廓也能用。但是要精准数据的话，就把那一段内的所有.FTp进行叠加，得到标准轮廓，以此建立模板.std。使用psrchive包含的如下命令：

pas -r (referrence file) filename 相位对齐

psradd -FTP *.FTp -o new.FTp -v轮廓叠加

[8]小跃变拟合

以上举例中，这个源J1731-4744的源轮廓非常好，然而好多数据的轮廓是不好的。拟合的时候，先选3个点，拟合出来后再加一个点，一个一个地往上加。小跃变不能一次拟合太多点！

拟合不出来也可能是.par文件的F0,F1数据不对，可以在JBO的数据库的网址里找到参考文献，看前边发生的glitch后的F0的具体数据，修正自己的.par文件的F0。

拟合的时候，尽量只拟合一阶导。如果一阶拟合不出来，再去拟合二阶导。有的小跃变可能一阶导有结果，二阶导F2就看不出来了。我们需要一段一段地拟合F1，每一段拟合的点不要太多！把图放大，看到明显的直线和后边的偏折。F1拟合是直线偏折，F2是曲线结构开始和中断。比如

 

下图为F2拟合glitch点的变化，曲线结构起点

 

这个图很平滑，所以那个曲线结构不会再有其他跃变了。如果没这么平滑，还是要一段一段的去拟合去试，看到底有没有。

再看以下这个例子：

青色处是有相位跳变，故而不连续

这个很明显，我以第一次glitch后第二次glitch前之一段的.par做的图，第二次跃变后的部分也能很好拟合，而第一次跃变前的部分只有散点，说明了：第二次跃变是小跃变，第一次跃变较大（频率变化有100左右，单位10-9Hz）

[9]精确测定glitch epoch：需要拟合计算

周世奇视频，00：51：00后讲了他的画图程序。再往后讲了怎么更精确地确定glitch发生的时刻：通过拟合频率及其导数的变化。

[10]总结

（1）     通过psrcat得到.par文件，通过.std和所有.FTp得到.tim文件，注意.par需要改的内容和.tim排序。我习惯以源的名字命名这两个文件，如J1731.par 和J1731.tim

（2）     之后通过命令画图tempo2 -gr plk -f xxx.par xxx.tim -epoch centre 逐段拟合，找到glitch

（3）     记录跃变前最后一次数据点的MJD和跃变后第一次数据点的MJD，简单地认为这两个值的平均值就是跃变发生的时刻

（4）     选择包含跃变前后的一段数据，ctrl+J保存.txt用于画残差图 命名为glitch1.txt

（5）     只保留跃变前的数据段，shift+P保存为.par文件，用于后边tempo2 glitch插件画图，也可以用于再次做残差图时不用从头再来。 我习惯命名：如第一次跃变前pre1.par 跃变后post1.par（大跃变可能拟合不出来，无法得到post.par，处理方法见后续章节）

那么就有了初始的.tim（需要）和初始的.par文件（不再需要）。对于比较复杂需要重复拟合处理的某次跃变，可以把.tim文件只截取该次跃变前后的部分。

以及每一次跃变前后的.par和每一次跃变的.txt。把以上这些文件放到fit文件夹中

好了，到这里我们就知道怎么拟合找到周期跃变了。那么跃变数据怎么得到，图怎么画呢？

 

5. 从跃变数据得到跃变参数

第四章中，我们已经找到了一个源的跃变现象，可能有小跃变或大跃变。那么具体的跃变参数应该怎么获得呢？

5.1外推法Extrapolated——小跃变

小跃变外推法只要把post.par中的F0,F1,F2与pre.par中的做差即可得到，非常方便。

但外推法不如拟合法精确，且不能得到跃变后恢复过程的参数。

5.2外推法Extrapolated——大跃变

大跃变无法得到post.par。新建一个空的post.par文件，把pre.post的数据复制进来，只保留重要的数据，删除以下内容：

（保存异常数据，是由于conda安装tempo2异常引起的，但在这里不影响，只要删除就好）

只保留：

注意F2，这里是把F2注释了。但是大跃变可能需要F2的值。可以把F2的值记为0，参数改为0，意为不拟合该项，误差值删掉：

对于大跃变，如果该次跃变已经被ATNF等数据库记录，那就把他们的变化值记下来，修改post,par里跃变前的F1,F0为跃变后的值，保存post.par。也可以去找该跃变链接的原论文内容， 会给出跃变幅度等参数。如果数据库里没记录，那就猜一个，多猜几次，看运气，可能得猜100次，拿修改后的去tempo2残差拟合，直到能拟合出来跃变后的自传模型。能拟合出来之后，就可以把跃变前后的F0和F1作差得到ΔF0和ΔF1了，这就是外推法的结果。

5.3拟合法Fitted

现在我们已经拥有了pre.par和post.par，用tempo2 plk做post.par的残差图。再次拟合得到跃变后的数据使其拥有更好的拟合效果。选取尽量多的数据点，把跃变后的数据整段拟合（大跃变把数据拟合到一起就行，可能不是直线，而是曲线或弧线，重点是要连续的线，不能是散点），保存为新的.par，覆盖post.par就行，直接仍然命名为post.par

复制pre.par为pre-x1.par作为拟合法的参考自传模型，在文件最后添加如下信息：

其中第二行，为跃变相变，默认给它0值。对于跃变幅度，如果采用ATNF等数据库的值，数据库中显示的为和值，记得乘以跃变前F0和F1才是ΔF0和ΔF1

对于跃变后恢复，一般大跃变才有恢复。看ATNF数据库：

这两列是恢复指数和恢复时间(?)，这两列没值的就是没恢复。没有恢复行为的，不用加恢复信息那两行。如果不确定有没有，先不加。

        把修改后的信息post-x1.par用tempo2画图，软件会考虑到跃变，使得跃变前后的数据在残差图中会显示为直线。拟合该段数据至得到直线（数据段包括跃变前，跃变历元，跃变后的数据，一般长度为几百天。比如说选跃变前后各300天的数据），保存为fit.par。在这个.par里，可以看到跃变的ΔF0和ΔF1被程序修正了，它的值比post-x1.par的值更精确。具体的信息自己查看。那么拟合法得到的跃变幅度就有了。

        如果post-x1.par不能被tempo2画出来跃变前后都被拟合的图，那说明很可能是有恢复过程没被考虑到，可以在该.par里添加恢复过程那两项，先添加一条，如果仍不行，判断是不是有更多的恢复过程，添加两条，可能有多条，直至能完全拟合为止。一般不会出现这种情况。恢复过程具体参数需要自己通过残差图猜测。

        如图为自动修正后的fit.par：

5.4结果处理，误差分析：误差与误差传递

我们得到的pre.par和post.par里有F0与F1值，只需要在相应文件夹下终端输入

grep -E “F0”  *.par和grep -E “F1” *.par就可以在终端里显示所有.par文件里F0和F1，及其误差的值。示例：

制作结果表格时，要把pre.par, post.par, fit.par所有拟合后的.par文件的数据抓取出来，需要：

Range: START-FINISH

(par) epoch : PEOCH

rms : TRES

No. of TOAs ; NTOA

F0, F1 (pre.par和post.par里的值)

和 用fit.par里的GLEP_F0乘以F0和GLEP_F1乘以F1

有个小办法：在.par文件所在目录的终端中输入命令：

 grep “F0\|F1\|START\|FINISH\|PEOCH\|TRES\|NTOA” *.par > J1731.txt

抓取这一堆参数，保存在一个文本文档里，如命名为J1731.txt，用excel打开时转化为表格，记得选取空格和冒号：都为分隔符！

要考虑误差，怎么通过误差传递公式得到结果呢？

        如果是拟合法，那Fit.par中有和的误差。如果是外推法，就是preF0和preF1、postF0和postF1及其误差，再利用公式(1)得到和

            使用excel处理，以F0为例，有pre.par的preF0和post.par的postF0,记为x和y。有preF0_error和postF0_error，记为Ux和Uy。则

                                                     

利用误差传播公式，可由F0和F1的误差得到跃变幅度的误差值。误差公式如下：

由上式可得误差，但值为相对误差，需要乘以结果值，即x的误差为

得到：

选择这两个数据单元格，右键，设置单元格格式，选其他，输入0.00(00)即为两位小数的以括号表示误差的办法(?)

到此，我们已经通过外推法或拟合法得到了跃变幅度

6.  结果与结果图

6.1写论文需要哪些数据结果和图表？

对于一个小幅度的周期跃变现象，在论文中我们通常需要给出如下的数据图：

        该图需要包含残差图(a)，自转频率变化量在跃变前后的变化(b),和自转频率一阶导的变化。图中红色是跃变发生的历元(glitch epoch)

        对于大幅度的周期跃变，通常需要以下的数据图：

大跃变前的自转模型不能拟合大跃变后的自转，因此不画计时残差图。需要画自转频率变化量在跃变前后的变化图(a)和它的放大图(b),以及自转频率一阶导数的变化图。

        以上数据图中，有时我们还需要画自转频率二阶导数F2的变化情况图。

同时完整的glitch描述还应该包括以下三个表：（以下例表来自于周世奇的工作[M]）

        (a)    所采用的源的参数和数据宽度‘

各项内容为脉冲星名(PSR)、赤经、赤纬、该表中各项参数测量对应的历元、自转周期、自转周期一阶导数、色散量、特征年龄、表面磁场、数据宽度。除了数据宽度以外的所有值都可以在ATNF的psrcat星表的数据库中查找到。数据宽度是指你处理该源的周期跃变所采用的数据的时间跨度，单位为年。

        (b)   周期跃变发生前后的计时参数

Int为数据计数。第一行1 - 2为第一次跃变后到第二次跃变前的自传模型(文中所述该源第一次跃变已在ATNF中记录，故只列出了第一次跃变后的数据)，即为pre2.par中的数据。包括自转周期，自转周期一阶导，二阶导，该.par的历元epoch(见.par文件中的Pepoch值)，拟合得到该自传模型.par的数据的MJD跨度(见.par文件中的start和finish的MJD值)，拟合得到该自传模型的数据的数量(见.par文件中的NTOA值)，以及拟合时的rms(应该为.par文件中的tres值,.如果保存的不对(比如拟合完了又重新框选)，这个值会有问题。也可以这样：在保存.par文件的时候，最后拟合，tempo2界面上显示的rms)

同理2 – 3为post2.par的数据（注意post2.par与pre3.par是否为同一自传模型，可以为同一个自传模型，只要拟合良好）

        (c)跃变参数

        表中数据为脉冲星名、跃变次序、跃变发生的历元。“New?”指的是该次跃变是否为他的工作新发现的跃变。我们可以忽略此项，因为我们写的不是硕士毕业论文。Extrapolated指的是外推法，内容为我们第五章得到的跃变幅度。Fitted指的是拟合法，其中还包含了跃变后的恢复过程的参数Q和，以及使用拟合法时采用的数据点数量和跨度以及拟合res。

        我们也可以简单表述，使用如下表格：

只提供一种方法的处理结果，并声明是哪种方法。

6.2使用插件gnuplot画图预览

在tempo2 plk图形界面ctrl+j保存的glitch.txt可以在命令行里简单用cat filename来查看，也可以使用gedit来编辑，文件格式如下：

未知数据  || 残差  || 残差误差 || toa   || 未知数据 ||中心频率||观测数据名||文件保存格式

可以使用gnuplot命令画图查看, 没有安装的话安装：sudo apt-get install gnuplot 然后在命令行中输入gnuplot进入插件。之后，输入plot “filename.txt” u 4:2 意思是按第4列和第2列作图：

纵坐标可能随数据值大小而有不一样的数据单位，可能是ms或s 。图中为ms

6.3作图

6.3.1 tempo2 glitch插件作图

Tempo2 glitch插件是tempo2自带的作图插件，需要使用.dat文件作图。

        在tempo2 -gr plk -f xxx.par xxx.tim -epoch centre命令后，按Ctrl+R可以在图像中选点。实线是起始时刻，虚线是终止时刻。鼠标左键选点，右键终止。 我们画图的话，就在这一段数据中取几段，然后得知的变化来作图。此时使用的.par文件为pre.par（使用外推法得到跃变幅度）或fit.par（使用拟合法得到跃变幅度）

        一段一段的选点，比如第一段选十个点，然然第二段和第一段5个点重叠，即1-10，5-15这样来取点。Ctrl+R选取数据时，每段数据小于3个点时无法拟合，每段数据只有3个点时得到的结果只有F0，大于等于四个点时才能得到F0和F1（以及它们的误差）。即最少取5个点。由周世奇的报告PPT：（然后拟合，根据卡方检验，要拟合的参数个数，要求数据量需要大于参数个数+1，即拟合F0,F1,F2,就要求数据点数量最少5次观测（4+1=5,因为考虑自转模型当中有(arbitrary phase)，这也是一个参数个数），小于5次观测则拟合结果是错误的。）（到底需要几个，我有些迷惑。你看regions.dat里会显示，如果数据点太少它会报错）

如下为选点（示意图）：

选完后，它会自动在文件目录下生成一个regions.dat文件。我习惯命名为 R_g1.dat 意思是ctrl+r处理glicth-1得到的.dat文件

文件长这样：

        作图的完整命令为 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 该命令会生成result.dat文件，同时得到图像

        命令tempo2 -gr glitch -h会提示此命令的参数，（命令中没有括号）。glitch插件画图标注开始的时间是为了让画图的时候横坐标初始值为零。标注glitch时间，有两个目的：1、是画出glitch epoch 这条虚线。2、是为了确定pre-glitch 部分，然后该画图程序计算出此部分的自旋模型，然后跃变后的频率减去这个模型值，就得到图中的Δ𝜈（画Δ𝜈图时，需要用一阶多项式拟合跃变前数据，然后用所有数据减去这个值，才是Δ𝜈

        其中 font 1.2及后边那些是画图参数， p6 p4是plot参数，在help里可以看到, p4是跃变前的Δ𝜈。大跃变加上p5。p5的策略是，跃变后的Δ𝜈减去跃变后Δ𝜈的均值

.par用跃变前的.par数据

帮助命令如下：

我们的结果图需要画类似如下的图。如果是大跃变，再加上减去平均值的（相当于放大的）Δ𝜈

终端里输入命令开始画图，它有好几种保存文件的格式，按问号’ ? ‘ ，它会显示你可以保存的文件格式：

这时候我需要输入 filename.yyy/yyy   yyy即为要保存的文件格式，比如all.ps/ps那么我就保存了文件名为all.ps的.ps文件

Linux可以用命令evince xxx.yyy来查看图片，如：evince all.ps

但是这样画图会存在一些误差，也不方便对图形的调整。因此我们可以使用matlab来作图。

6.3.2 matlab作图

tempo2 glitch命令画图之后，会生成result.dat文件。文件中的数据点就是我们图里所对应的数据。如下：

MJD     || F0          || F0误差      || F1        || F1误差 || 后两列数据标号，无用

        (自己画图的时候，本身tempo2 glitch出来的数据比较平滑，但是你画出来的可能不很平滑，这时候你就需要这么干：

打开regions.dat看一下，第一个点的起止时间，复制它的起始时间，把画图用的参考.par里start改成这个，再运行tempo2 -gr plk -f xxx.par xxx.tim -epoch centre 拟合后，shift+P保存分段得到许多段的.par（就像ctrl+r那样分段），然后通过命令grep -E “F0|F1” 5 *.par 来抓取数据这一些.par里的F0，F1（？以及他们的残差）数据值)

        现在我们有了g1.dat（自动保存时命名为result.dat）。或者是提取一系列.par文件得到的F0，F1值，也按照result.dat里数据的顺序制成表格，再画图。画残差图、ΔF0和F1的变化图。

        使用的代码如下：matlab代码尝试绘制脉冲星周期跃变数据图 - 哔哩哔哩 (bilibili.com)

画图的方法是：

        对于，需要拟合，先挑选出glitch发生之前的那一段的数据F0，然后我们手动对跃变前的那些F0数据进行线性拟合，得到相应自转模型。然后再用现在的F0减去该模型下的F0得到

        这个拟合算法的策略是，选择周期跃变前（如果有多次的画，选第一次的）的数据，然后用一阶多项式拟合他们。用真实的F0与拟合曲线上对应到达时间点的（预测）F0值做差，这个差就是Δv

        画计时残差图时要注意：残差数据和F0数据范围不要相差太多。数据只取跃变附近的就行了，不用取太多。

        F1也需要做拟合曲线。分段选择跃变前与跃变后的数据，分别用一阶多项式拟合他们，并画在结果图中。F1的拟合曲线，看F1是否有线性趋势，没有的话就不要画拟合曲线了

        大跃变，无法拟合跃变后数据，不需要画残差图，但是deltaF0要有减去均值的贴近0的数据图（放大图），跃变前的各个值减去跃变前均值，跃变后的减去跃变后均值。

6.3.3总结

[1]得到J1731.tim J1731.par，见4.2.4节

tempo2 -gr plk -f xxx.par xxx.tim -epoch centre画图后：

shift+P得到的pre1.par和post.par（外推法），以及可能有的fit.par（拟合法）

crtl+J得到的glitch1.txt

crtl+R得到的R_g1.dat （生成时自动命名为regions.dat）

[2]之后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 得到的g1.dat（生成时自动命名为result.dat）和自动生成图g1.ps

[3]画残差图

如5.1节中的gnuplot画残差图，也可以自己画

[4]画和频率一阶导变化图

按照我的matlab代码画图

 

Supplement

最后，一些其他资讯，可以随便看看

由于b站专栏排版，这些超链接都失效了，我也懒得放链接了

[i]FAST开机测试：

这颗星星，献给他 :“中国天眼”FAST望远镜首次发现脉冲星！_百科TA说 (baidu.com)

来自新疆天文台的一些观测资讯：

[ii]脉冲星观测系统-中国科学院新疆天文台 (xao.ac.cn)

[iii] 25米射电望远镜-中国科学院新疆天文台 (xao.ac.cn)

以及该页面索引处的几个介绍网页

[iv]首页 | 天文学名词 | Astronomical Terms (china-vo.org)

[v]以及该网站首页首页 | 国家天文科学数据中心 | NADC (china-vo.org)

[vi]来自公众号的一些文章，看看其中的原理介绍即可，不要按他们的流程来，容易混乱：

下载并处理脉冲星射电数据 (qq.com)

脉冲星射电数据处理入门 (qq.com)

Introduction to Pulsar Timing (qq.com)

tempo2介绍脉冲星计时软件包TEMPO2的使用手册 - 知乎 (https://zhuanlan.zhihu.com/p/537042917)

这篇文章介绍了tempo2的使用方法，较为详细。以及：

使用Pulsar计时软件包Tempo2 - 哔哩哔哩 (bilibili.com)

Tempo2 和psrcat的一些help 我已经放在了word文件‘tempo2 psrcat-h’里

补充图片

1.psrcat的参数

2.tempo2操作键位

 Ps.最近听了戴子高老师的报告，反跃变对于FRB来说还是挺重要的