LAB巡天全天的中性氢数据的可视化，对天球作了Mollweide投影，注意tvscl里Startx和Starty以及Xsize和Ysize对于图和坐标线的对齐有关键性的作用。目前程序表示的只是积分强度，尚没有转化为柱密度。


PRO processLAB

fitsname='lab.fit';default file name
deal,fitsname

END


PRO deal,fitsname


;Distance = 140.0; distance in units of parcec
;Distance = Distance*3.086d18; distance in units of cm

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

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

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');
bzero   =   fxpar(head,'BZERO');
bscale  =   fxpar(head,'BSCALE');

;x=(findgen(nx)-crpixx)*cdeltax+crvalx;
;y=(findgen(ny)-crpixy)*cdeltay+crvaly;
;z=(findgen(nz)-crpixz)*cdeltaz+crvalz;

;area = Distance^2*(abs(cdeltax)/57.3)*(abs(cdeltay)/57.3)
                                                   ; physical area per pixel
;print, area
;print, area*1.67d-24*weight/2d33


;a=ptr_new(/allocate_heap)
a=mrdfits(fitsname,0,/fscale); read the data cube to an array
vchannel=abs(cdeltaz)/1.0e3 ; km/s
;a=readfits(fitsname)

;b=total(a,3); co-add the third dimension
;????
;vchannel=bw/(nz*1.0d0)/freq*ckm; the velocity width corresponding to a channel
                              ; units in km/s
;b=b*vchannel; intensity integrated with velocity (with unit: Jy km/s)
;????
;b=b*vchannel; antenna temperature integrated with velocity (units K km/s)
;b=b*vchannel*1.93e3*nu^2*(1.0/3.0)/A; column density (units /cm^2)
                               ; N_l=1.93*10^3*(g_l/g_u)*(nu^2/A_ul)int T dv
;b=ptr_new(/allocate_heap)
;/for test
;*b=total(*a(1:100,1:100,*),3)*vchannel*1.93e3*nu^2*(1.0/3.0)/Alu*area*1.67d-24*weight/2d33
                               ; mass distribution in solar mass
;for test
b=total(a,3)
;b=a[*,*,500]
bad=where(finite(b) eq 0, count)
if(count gt 0) then b(bad)=0.0
;b=(b+abs(b))/2.0
;b=b*vchannel*1.93e3*1.42^2/4.65e-17*(4.0/3.0)
b=b*vchannel*(4.0/3.0)
b=b/1.0e3
print,max(a)
maxb=max(b)
print,maxb
print,(min(a)-bzero)/bscale
colors=200
clrtablelength=colors





;color_array=findgen(256L)
;clevs=findgen(256L)

LoadCT, 5, NColors=colors, Bottom=1, /Silent
device,decompose=0

;loadct,5
;tvlct,[70,255,0],[70,255,255],[70,0,0],1



latmin=-90
latmax=90
lonmin=-180
lonmax=180
;position = [0.1, 0.1, 0.9, 0.75]
position = [0.1, 0.90, 0.9, 0.95]

margin=0.12
;margin=0.5
wall=0.03
xsize=18.8
aa=xsize/8.8-(margin+wall)
bb=aa*2d/(1+sqrt(5))
ysize=(margin+bb+wall+bb+wall)*6.8

;================================================
set_plot,'PS'
filename=fitsname+'_TV.eps' ; set the file name of the output ps file
device,file=filename,/ENCAPSULATED,/COLOR, BITS=8;,xsize=xsize,ysize=ysize
;tvscl,c(1:1000,1:1000)
;tvscl,b
map_set,0,0,/MOLLWEIDE,/ISOTROPIC,/HORIZON,/GRID
result=map_image(b,Startx,Starty,Xsize,Ysize,compress=1,LATMIN=latmin,$
LONMIN=lonmin,LATMAX=latmax,LONMAX=lonmax,scale=0.1)
result=bytscl(result,Min=0.01)
tvscl,result,Startx,Starty,XSIZE=Xsize,YSIZE=Ysize
;tvscl,result
;print,Startx,Starty
;TVimage,BytScl(result, Top=99)
;lons=indgen(360/20+1)*20*(-1)+180
lons=indgen(360/45+1)*45-180
lonnames=strtrim(-lons)
;print,lonnames
;map_grid,latdel=10,londel=20,lons=lons,color=0.30*!d.n_colors,/LABEL,/HORIZON
map_grid,latdel=20,londel=20,lonnames=lonnames,lons=lons,$
color=0.80*!d.n_colors,charthick=2,glinethick=3,/LABEL,/HORIZON
;;colorbar,ncolors=256,POSITION=[0.85, 0.15, 0.95, 0.95]
;ticks=strtrim(sindgen((clrtablelength/5+1)*5)-(clrtablelength/2),2)
;Colorbar, Range=[0, Max(b)], Divisions=10, $
;Minor=5, NColors=colors, Bottom=1, $
;Position=position,Charsize=1,ticknames=tickes
timesstr=textoidl('\times 10^{4} K\cdot km/s')
Colorbar, Range=[0, Max(b)], Divisions=10, $
Minor=5, NColors=colors, Bottom=1, $
Position=position,Charsize=1,title=timesstr
;xyouts, 90, 180,timesstr,charsize=1.5,charthick=1.5,color=0;.80*!d.n_colors
;tlb=widget_base()
;labeltext='shit'
;label=widget_label(tlb,value=labeltext,ysize=40,units=0)
;widget_control,tlb,/realize

device,/CLOSE
;================================================


;spectra of each pixel
;print,"enter the coordinate of the pixel you need:"
;print,"ix(1 ~",nx12,")"," iy(1 ~",ny12,")";
;read,ix,iy;
;plot,z,a[ix,iy,*];

;find peak

;(pixel number - crpix)*cdelta+crval

END

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

上一篇：查看光谱的IDL程序
下一篇：用IDL生成某个分布的二维数组