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');
;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
;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,*];