""" Casa script for the manipulation of fits line cubes produced with ProDiMo. Does a beam convolution of the line cube with a beam and also does the continuum subtraction. Further the following observables are produced: - channel maps for the line (with and without continuum subtraction) - momements (zeroth and first moment) - a position-velocity diagram - A spectral line profile - an azimuthally averaged radial intensity profile. (requires task casairing see https://www.oso.nordic-alma.se/software-tools.php) The script works only for one single line. But it can be (easily) adapted to e.g. run for all ProDiMo line cubes of a single model. For some reason, I do not understand, the script cannot be run twice within one CASA session. I recommend to run the script directly from the command line (casa -c convsim.py) in the directory with the ProDiMo line cube. The script was only tested with CASA 4.7.2. The produced output files (and there names) are compatible to the plotting routines for ALMA simulations in the prodimopy package. (see https://bitbucket.org/cheesyog/prodimopy) There is no guarantee for correctness of this script and its results. If you find any bugs or have suggestions for improvement plese contact Christian Rab. """ """ Setting of the input parameters. These values need to be set so that the script can run. """ # what line to take from prodimo, and output directory and file/prefix pline="001" proj="L"+pline+"CONV" # what franction of the chanales should be used for the continuum e.g. 3 mean 1/3 at each side of the line nfrachan=3.5 # the beam used for the convolution beam=["0.1473167527972E-03deg","0.1222582832068E-03deg",str(0.8938728915312e2*0.01745329252)+"rad"] # position of the target (musst be in degree like in the fits file) # Provide some coordinates for the for the center of the image # FIXME: could use some abitrary values, put this way it is easier. RA="11h01m51.9067s" DEC="-34d42m17.032s" # the output image region bregion="centerbox[["+RA+","+DEC+" ], [10.0arcsec, 10.0arcsec]]" # length and velocity range for the PV diagramm pvlength="8.0arcsec" pvrange=("range=[-2.3km/s, 2.3km/s]") # rms for each channel of the cube, # if there is no real noise, set it to a very low to avoid numerical noise cuberms="1.e-6" """ Here starts everything. It should not be require to change anything here, but there might be cases where it is necessary. All the files in the output directory startgin with the proj prefix are deleted before anything is done! """ ################################################################################################### # Some Preparation project=proj fprefix=proj+"/"+proj fitsfile="LINE_3D_"+pline+".fits" # fits input file print "Output directory: ", project # wipe the directory any create it (if it does not exist) os.system("rm -rfv "+proj+"/"+proj+"*") os.system("mkdir "+proj) nchanfits=imhead(fitsfile, mode="get", hdkey="shape")[2] ################################################################################################### # Beam Convolution ia.fromfits(infile=fitsfile) imconv=ia.convolve2d(major=beam[0],minor=beam[1],pa=beam[2],outfile=fprefix+".conv") imconv.done() ia.close() # set the center coordinates imhead(imagename=fprefix+".conv",mode='put',hdkey='CRVAL1',hdvalue=RA) imhead(imagename=fprefix+".conv",mode='put',hdkey='CRVAL2',hdvalue=DEC) ################################################################################################### ## Continuum subtraction ## also outputs the continuum (for each channel) default("imcontsub") #nc=nchanfits/nsplit nc=nchanfits nrange=max([1,int(nc/nfrachan)]) fitspwstring="0~"+str(nrange-1)+","+str(nc-nrange)+"~"+str(nc-1) linefile=fprefix+".conv.contsub" contfile=fprefix+".conv.cont" os.system("rm -rf "+linefile) os.system("rm -rf "+contfile) imcontsub(imagename=fprefix+".conv", linefile=linefile, contfile=contfile, chans = fitspwstring) # average over all channels for the continuum (only one channel is needed) immoments(imagename=contfile,moments=[-1],outfile=contfile+".mean") ################################################################################################### # Image moments (zeroth and first) default(immoments) immoments(imagename=linefile,moments=[0],outfile=linefile+".integrated") # define a mask to make it look nicer, simply avoid very small values threesigmask="'"+linefile+"'"+" > ("+cuberms+"*3)" immoments(imagename=linefile,moments=[1],outfile=linefile+".mom1",mask=threesigmask) ################################################################################################### # Postition velocity Diagramm # read the PAobj from the fits file and use it # FIXME: the +90 is required, but I do not know really know why, # something with the orientation is strange pvpa=str(imhead(imagename=fitsfile,mode='get',hdkey='P_PAOBJ')+90.0)+"deg" impv(imagename=linefile, outfile=linefile+".pv", mode="length", center=[RA,DEC],length=pvlength,pa=pvpa, chans=pvrange,width=3,overwrite=True) ##################################################################################################### # Exports all images/cubes to fits files exportfits(imagename=fprefix+".conv", fitsimage=fprefix+".cube.lpc.fits",dropdeg=True,overwrite=True) exportfits(imagename=linefile, fitsimage=fprefix+".cube.fits",dropdeg=True,overwrite=True) exportfits(imagename=linefile+".integrated", fitsimage=fprefix+".cube.integrated.fits",dropdeg=True,overwrite=True) exportfits(imagename=linefile+".mom1", fitsimage=fprefix+".cube.mom1.fits",dropdeg=True,overwrite=True) exportfits(imagename=linefile+".pv", fitsimage=fprefix+".cube.pv.fits",dropdeg=True,overwrite=True) exportfits(imagename=contfile+".mean", fitsimage=fprefix+".cont.fits",dropdeg=True,overwrite=True) ################################################################################################### # Spectral line profile #specflux(imagename=fprefix+".conv", logfile=fprefix+".conv.specprof", overwrite=True) specflux(imagename=linefile, logfile=fprefix+".cube.specprof", overwrite=True) ################################################################################################### # azimuthally averaged profile casairing(image=fprefix+".cube.integrated.fits",resultfile=fprefix+".cube.integrated.radial",nrad=50) casairing(image=fprefix+".cont.fits",resultfile=fprefix+".cont.radial",nrad=50) # remove all the CASA stuff os.system("rm -rf "+proj+"/*/") os.system("rm "+proj+"/*.last")