在13CO数据里找出来的团块速度都很小，总觉得有点不对，所以想看看每个波束的光谱的形状。Gaia用起来有些不方便，还是IDL比较顺手，只是很多功能不是现成的。

这个程序基本就是把三维数据从数据文件中读出来，用户给一个坐标，程序把对应的光谱取处理作高斯拟合，然后画出数据和拟合结果，并给出线宽。

PRO viewfits

fitsname='xxx.fits'
a=mrdfits(fitsname)

head=headfits(fitsname);read the header of the fits file to a vector

bw      =   fxpar(head,'BW'); band width (MHz)
freq    =   fxpar(head,'LINEFREQ'); central frequency (GHz)

nx      =   fxpar(head,'NAXIS1'); number of elements in the first dimension
ny      =   fxpar(head,'NAXIS2');
nz      =   fxpar(head,'NAXIS3');
crvalx  =   fxpar(head,'CRVAL1'); reference value of the first dimension
cdeltax =   fxpar(head,'CDELT1'); increasement of the first dimension
                                ; in units of degree, when calculate physical
                ; scale, must changed to arcdegree
crpixx  =   fxpar(head,'CRPIX1'); reference position of the first dimension
crvaly  =   fxpar(head,'CRVAL2');
cdeltay =   fxpar(head,'CDELT2');
crpixy  =   fxpar(head,'CRPIX2');
crvalz  =   fxpar(head,'CRVAL3');
cdeltaz =   fxpar(head,'CDELT3');
crpixz  =   fxpar(head,'CRPIX3');
f = (dindgen(nz)+0.5-crpixz)*cdeltaz+crvalz; m/s
f = f/1.0e3; km/s
width = cdeltaz/1.0e3

expr = 'P[0]+GAUSS1(X,P[1:3],/PEAK)'

status = 'init'
while (status ne 'exit') do begin
print,'enter x y coordinates:'
read,x,y
c=a(x,y,*)
maxc = max(c)
index = where(c eq maxc)
start = [0.1,f(index),width,maxc]
result = MPFITEXPR(expr,f,c,c*0.0+1.0,start)
sigma = result[2]
plot,f,c,psym=4,symsize=1.0
oplot,f,result[0]+gauss1(f,result[1:3],/PEAK),color=100,thick=3
strsigma = textoidl('sigma =')+strtrim(sigma)
xyouts, f[index],c[index]*0.8, strsigma,charsize=1.5,charthick=1.5
print,'enter command:'
read,status
endwhile
END

转载本文请联系原作者获取授权，同时请注明本文来自钱磊科学网博客。
链接地址：https://blog.sciencenet.cn/blog-117333-377218.html

上一篇：小记IDL读FITS文件（数据标定）
下一篇：IDL画LAB巡天全天中性氢分布图