php极简wiki,科学网—CLASS极简教程 - 钱磊的博文

CLASS是gildas中的一个软件,之前对CLASS有所接触,但仅限于相关的文件格式转换(http://blog.sciencenet.cn/home.php?mod=space&uid=117333&do=blog&id=303542)。最近课题组需要用CLASS处理数据,所以我们在阅读CLASS手册(http://www.iram.fr/IRAMFR/GILDAS/doc/pdf/class.pdf)以及一个网上简易教程(http://www.mpifr-bonn.mpg.de/staff/bparise/gildas/2010-pety-class-demo-max-planck.pdf)、张智昱博士的网站(http://159.226.71.202/wiki/FAST_2013)一些网页(http://www.astro.umass.edu/~fcrao/library/manuals/map.html)以及德令哈观测站编写的CLASS使用手册(http://www.radioast.csdb.cn/tools.php)的同时做一些笔记,形成此CLASS极简教程。有一些命令是新版本的用fortran 90重写的CLASS里才有的,在这种情况下,命令前面的提示符用LAS90>表示。

本手册不定期更新,牛人无须看。

1. FITS文件(.fits)转换为gildas文件(.gdf, .30m)       天文中一些不同领域定义了不同的文件格式,不同的数据处理软件有时候也有自己的数据格式,比如说gildas相关软件用的就是自己的文件格式(后缀为.gdf)。不同的领域交流时通常使用标准的FITS文件,所以将FITS文件转换为gildas文件是使用CLASS的一个基础操作。.30m文件是CLASS文件默认的后缀名。

最简单的转换方法是使用vectorfits工具(vector包里的fits工具)。在shell里敲class,回车,打开class后,使用如下命令

LAS> vectorfits spectra.fits to spectra.gdf

就可以把FITS文件spectra.fits转换为gildas文件spectra.gdf。不过要注意,FITS文件需要满足一定要求,否则转换出来的gildas文件也不能为CLASS所用。比如,如果spectra.fits是数据块(data cube),那么转换出来的文件CLASS是无法直接使用的。

CLASS能接受的文件是光谱文件,如果有FITS格式的光谱文件(spectra),那么可以如上述将文件转换为gildas格式,也可以进行如下操作将fits文件转换为CLASS能使用的文件(.30m),其中spectra.fits是一个FITS格式的光谱文件。

LAS> set angle sec

LAS> lasfits read spectra.fits

LAS> sic delete spectra.30m

LAS> file out spectra.30m new

LAS> write 1

LAS> file in spectra.30m

LAS> get f

LAS> plot

注意到,读入fits文件用的是lasfits,和vectorfits是不同的。上述操作第一步是设置角度的单位,就文件转换而言是不必要的。第二步就是读入FITS格式的光谱文件。第三步是删除已有的CLASS文件,第四步是写入新的CLASS文件。之后就是读入新生成的CLASS文件,画图。

对于数据块fits文件(比如叫做qhz_cube.fits),可以使用lmv读入、转换。具体如下

LAS90> file out qhz.14m m

LAS90> lmv qhz_cube.fits

这样,数据块fits文件qhz_cube.fits就转换为calss文件qhz.14m了。

2. gildas文件(.30m)转换为FITS文件(.fits)

这个操作就是将前面一个操作倒过来,不过,在这个情况下,vectorfits就不顶用了,得用lasfits,将上面一个操作中的后一个方法倒过来。关于此操作,还可以参看(http://iram.fr/IRAMFR/GILDAS/doc/html/class-html/node15.html )。

LAS> set angle sec

LAS> file in spectra.30m

LAS> find

LAS> get f

LAS> set fits mode spectrum

LAS> lasfits write spectra.fits

其中get f是读出第一条光谱。

在CLASS90里可以将一系列光谱文件转换为FITS二进制表文件(生成FITS数据块文件还不知道如何实现)。

LAS90> file in spectra.30m

LAS90> find /range -200 200 -200 200

LAS90> fits write spectra.fits /mode index

3. gildas文件(.30m)转换为文本文件(.dat)

最简单的光谱文件就是在文本文件里写两列数,第一列是波长或者频率,第二列是强度。将CLASS文件转换为文本文件可以方便其他一些程序进行处理。具体操作如下,其中spectra.30m是输入的CLASS文件

LAS>file in spectra.30m

LAS>find

LAS>get f

LAS>sic output spectra.dat

LAS>for i 1 to channels

LAS: say ’rx[i]’ ’ry[i]’ /format g12.4 g12.4

LAS: next

LAS>sic output

4.读入文本文件(.dat)中的光谱并画图假设spectra.dat是含有光谱的文本文件。LAS> dev xl w

LAS> greg1column x 1 y 2 /file ’’spectra.dat’’

LAS> model y x /regular

LAS> plot

5. 减基线

减基线是谱线数据处理中经常需要进行的操作,不过光看手册的话会比较迷惑,参考对Effelsberg望远镜数据处理的简介(http://www.astro.uni-bonn.de/~rcbruens/tutorials/effberg/HI.html),我们总结出以下操作。

LAS> file in spectra.30m

LAS> find

LAS> get f

LAS> plot

LAS> set cursor on

LAS> set window

这个时候,当把鼠标放到画图窗口的时候会出现一条横线和一条竖线,交点是鼠标的位置。不看说明的话,一般都不知道怎么退出这种模式(答案是按“e”键)。这个模式是让你到画图窗口标明一下强发射线的位置,鼠标到某条发射线(吸收线)的左边,按空格键,到右边,按空格键,以此类推,对所有发射线(吸收线)进行此操作。这个过程你会感到心虚,因为你看不到任何变化。不要急,完成上面步骤之后,按“e”键退出。(也可以直接指定窗口的位置,比如LAS> set window -120 -100 0 20。其中的数就是各窗口左侧和右侧的坐标。)

LAS> draw window

此时就可以看到在下方出现一个矩形,标明发射线的位置(反射性所在的频率窗口(波长窗口))。然后就可以拟合基线了(比如用三个自由参数的,也就是二次多项式拟合,并画出基线)。

LAS> base 3 /plot

这时就看到拟合的基线了。然后

LAS> plot

就可以看到减完基线的谱线了。

结合下面提到的用CLASS90将光谱排好画出来

LAS90> find

LAS90> load

LAS90> plot /index

可以在二维上定义减基线的窗口

LAS90> set window /poly 1   ! Define 1 polygon

6. 构建数据块

成图观测最终能方便使用的数据时数据块,即三维数组,三维分别是赤经、赤纬和速度(频率)。但是成图观测的原始数据是一系列光谱,把这些光谱按坐标排好才能形成数据块。具体构建数据块的步骤如下(其中demo.30m(http://www.mpifr-bonn.mpg.de/staff/bparise/gildas/demo.30m)是相关数据)。

LAS90> file in demo.30m

LAS90> find

LAS90> list

这一步是看看文件中有哪些观测。

LAS90> find /tel 30M-V01-A100

寻找一批标号为30M-V01-A100的观测,然后检查一致性

LAS90> consistency

如果数据满足一致性,就可以进行以下操作

LAS90> let name thecube

LAS90> let type lmv

LAS90> table ’name’ new

LAS90> xy_map ’name’

至此就生成了数据块文件thecube.lmv,然后可以查看

LAS90> go view

LAS90> vectorfits outfile.fits from infile.lmv

7. 对多条谱线进行减基线操作

通常一个数据文件中有很多条谱线,

如果每条谱线都手动进行减基线操作将会非常费时费力。通常一次观测的谱线窗口相差不大,所以可以进行批处理。下面是从读入数据、设定窗口、到循环减基线的脚本。最后画出了所有谱线。

file in spectra.30m

set window 10 20

find /number 59 100

for i 59 to found

get i

base 3

next

stamp

将上面的语言写到一个叫baseline的文件中,在class程序里

LAS> @baseline

就可以执行了。

8. 谱线拟合

减完基线就可以进行谱线拟合了。读入数据画图的步骤和前面相同(qhz.14m是青海站观测的数据http://www.radioast.csdb.cn/tools.php)

LAS> file in qhz.14m

LAS> find

LAS> get f

LAS> plot

下面就是谱线拟合了,假定有四条谱线

LAS> line 4

此时就会出现和set window之后类似的情况,不过要简单一点。四条谱线共需要八个位置来定义,鼠标点八次(每条谱线的左右),然后自然就回到命令窗口。

LAS> min

就开始拟合,完了可以查看

LAS> vis

9. 存ps、eps文件画完图,如果不希望退出程序图就消失,那么就需要把图存在文件中。一般把图存成ps或eps格式,发表文章时也可以用。使用下面的命令就可以完成这个操作

LAS> hardcopy spectra.ps /dev ps fast

用来存ps文件。

LAS> hardcopy spectra.eps /dev eps color

用来存eps文件。

10. 画p-v图画p-v图需要取数据块中沿某个空间方向的一片数据,这可以通过strip命令实现(不过首先要通过find命令找到这一片数据),具体如下。

LAS> file in qhz.14m

读入数据,然后设定角度的单位为角分

LAS> set angle min

结下来用find命令寻找数据片

LAS> find /offset * 5

找y方向偏离5角分的一片数据(x-v二维数组),然后用strip命令将此数据存成一个文件

LAS> strip qhz.gdf

随后就可以将这片数据读入,画p-v图了。

LAS> define image a qhz.gdf read

设定画图的范围

LAS> limit 0 50 * * /rgdata a

然后画出数据,对数据进行重新的线性标度

LAS> greg2plot /scal linear 0 20

然后画出等值线

LAS> level 1 to 10

LAS> rgmap /abs 2

然后就是画轴和做标记

LAS> lab "Velocity (km/s)" /x

LAS> lab "gDDEC. (arcmin)" /y

LAS> axis xl /unit r

LAS> axis xu

LAS> axis yl /unit m

LAS> axis yr

将光谱排好画出来CLASS90还可以这么做

LAS90> find offset * 5

LAS90> load

LAS90> plot /index

11. 定义数组在CLASS中可以定义数组,并对数组进行运算,这对于数据处理是非常必要的。定义10个元素整数数组(数组的指标从1开始)可以用如下命令

LAS> define integer a[10]

如果要定义双精度数数组(我不知道单精度数怎么弄,囧rz……)可以用如下命令

LAS> define double b[10]

定义二维数组可以用如下命令

LAS> define double c[10,10]

多维数组也可以类似定义。经本人研究,最多能定义4维,比IDL的8维的上限要弱小一些。不过亲~,4维一般够用了吧。

12. 数组赋值、查看等操作       数组赋值可以用let命令

LAS> let a[1] 1

查看数组值可以用exam命令

LAS> exam a[1]

13.去除RFI(坏通道)

基本想法就是用上面提到的并排画光谱的操作将光谱排好,如果有RFI(坏通道),那么某几个通道的值就会很大,超过采样最大值。这个之后需要把这些通道用高斯噪声填满,这就是去RFI(坏通道)的基本方法。

LAS90> file out a100-fill single /over

LAS90> find

LAS90> for ient 1 to found

LAS90>    get next

LAS90>    fill -120 -100 60 75 /noise    ! File contaminated channels with Gaussian noise

LAS90>    write

LAS90> next ient

14. 思考题      如何画任意方向的p-v图?

spectra.30m

15. 画积分强度图

本段是向张智昱、潘之辰学习而来的。首先要对一个数据文件作积分(在CLASS里):

LAS90> file in co65_final.apex

LAS90> find

LAS90> get first

LAS90> set unit v

LAS90> plot ! 查看一下谱线的速度范围

LAS90> print moment 0 500 /out moment.dat

LAS90> exit ! 退出CLASS

打开GREG,读入之前得到的数据文件

GREG> column x 2 y 3 z 4 /file moment.dat

GREG> set box m

GREG> lim /rev x

GREG> random 100 /bl -1 !将数据转换为二维图像

GREG> plot !画积分强度图

GREG> box /abs !画坐标框

GREG> wed ! 画颜色棒

张智昱指出,在CLASS里是可以直接调用GREG里的函数的,在函数名前面加上greg就可以了。

16. 快速查看谱线常用命令

LAS90> file in xxx.cso                                           ! 后缀名一般是观测用的望远镜的名称,格式是一样的。

LAS90> find /source xxx /telescope "yyy"             ! 寻找用yyy仪器观测的xxx源的数据

LAS90> list

LAS90> get nnn                                                   ! 读入第nnn条光谱

LAS90> set nomatch                                            ! 不匹配坐标

LAS90> set u v                                                     ! 将横坐标的单位设为速度(而不是频率)

LAS90> average                                                   ! 把上面找到的所有谱线平均(对找弱谱线尤其有用)

LAS90> set mode x a b                                        ! 取一段x坐标a到b的数据

LAS90> set cursor on                                           ! 设置鼠标可以在窗口读数

LAS90> set window                                              ! 设置窗口,用鼠标读数(原来是靠左键点,现在靠按n,退出读数读取还是按e)

LAS90> plot

附录

一、脚本

1. fitscube2bur.class

FITS数据块文件转换为bur文件,用法(假设FITS文件为data.fits):

LAS90> @fitscube2bur data

=======================================

sic delete &1.bur

file out &1.bur m

lmv &1.fits

=======================================

2. bur2fitscube.class

将bur文件转换为FITS数据块(先转换为lmv文件,然后再转换为FITS),用法(假设bur文件为data.bur):

LAS90> @bur2fitscube data

=======================================

file in &1.bur

find

consistency

let name thecube

let type lmv

table ’name’ new

xy_map ’name’

vectorfits thecube.fits from thecube.lmv

=======================================

转载本文请联系原作者获取授权,同时请注明本文来自钱磊科学网博客。

链接地址:http://blog.sciencenet.cn/blog-117333-567008.html

上一篇:思考题(十三)赤道附近的望远镜有什么优势?

下一篇:科学家、科研技术员和科研工人

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值