#!/bin/csh -f # # ccombine.csh - CARMA STING track combination script (Tony Wong, UIUC) # # History # 2009aug10: Initial version # 2010aug12: Allow specification of either vlast or nmaps; support vmom0 # # 20100815: Copied from Tony # log files will be renamed with prefix # default showplot changed to # # 20100816: calculation more robust # do beamfit before MOSSDI2 # # 20110320: add new keywords related to MOSMEM # : choose the deconvolution task (default: mossdi2) # : extra keyword for mosmem # e.g.: exmem="model=model.sdi flux=10" # <- the model map is an intial guessing, # in unit jy/pix # exmem="default=default.map options=doflux,dofactor" # : sd map in jy/beam for joint deconvolution # (the script will help you regrid) # : sd beam for jointdeconvolution (from imgen) # (don't need to match the size of dirty beam, # position independent) # 2011jul13: rewritten to allow gain normalization & noisemap calculation #-------------------------------------------------------------------------- # Basic params set files= set source= set interact=y set showplot=n set purge=n # Mapping and deconvolution params set vfirst= # velocity of first channel set delv= # channel spacing in km/s set vlast= # velocity of last channel, overrides nmaps set nmaps= # number of maps to generate set imsize=513 # for INVERT set cell=.5 # for INVERT set robust=0.5 # for INVERT set exinvert= # extra parameters for INVERT set rmsfac=3 # rmsfac for mossdi set clnbox= # deconvolution region, senmsk if unset set maxsen=3 # sensitivity contour level for senmsk set msklev=2 # sensitivity contour level for final cube set bmaj= # beam major axis to override default, arcsec set bmin= # beam minor axis to override default, arcsec set bpa= # beam position angle to override default, deg # Moment params set fwhm=9 # FWHM of smoothed map for masking method set mom0fun=sqr # Use n^2 (sqr) or log contours from mom-0 plots set mom0lev= # First contour level, default is 3-sigma set loglevs=0.2 # logarithmic contour level spacing in dex set mom1del=20 # Spacing between velocity contours, in km/s set fluxclip=4 # min allowed flux (Jy/bm km/s) for gaufits set velclip=20 # max allowed velocity error (km/s) for gaufits set nspan=3 # number of channels from peak to integrate for wmom0 #set mrms= # cheat srms in moment.csh set psbox= # custom box for plotting set decon='mossdi2' # deconvolution task, alternative: mosmem set exmem= # extra keyword for mosmem set sdmap= # sd map for joint deconvolution set sdbeam= # sd beam for jointdeconvolution ### LOCAL SETTINGS alias date 'date "+%m/%d/%y %T %Z"' set version=2011jul13 set log=ccombine.log ### OVERRIDE ABOVE PARAMETERS IF GIVEN ON COMMAND LINE foreach par ( $* ) set check=`echo $par | awk -F= '{print NF}'` if ( "$check" >= 2 ) set $par end ### CHECK THAT CUBE IS FULLY SPECIFIED if ($vfirst == "" || $delv == "") then echo "### FATAL ERROR: both vfirst and delv must be specified" exit 0 endif if ($nmaps == "" && $vlast == "") then echo "### FATAL ERROR: either nmaps or vlast must be specified" exit 0 endif if ($vlast != "") set nmaps=`calc -i "($vlast-($vfirst))/$delv+1"` #============================================================================ # OPEN LOGFILE IF NEEDED echo "" if (-e $log && $interact != 'n') then echo -n "$log exists. Overwrite (o) or append (a) [default=a]? " set resp = $< if ($resp == "o") rm -f $log echo "" endif # BEGIN EXECUTION goto continue # MOVE 'continue:' around to where you want to work. # Shouldn't be necessary in standard operation. continue: # MAIN MENU if !($?start) then echo "Where would you like to start?" echo " 1. Make dirty maps and plot" echo " 2. Deconvolve maps and plot" echo " 3. Make moment maps" echo " 4. Make Gaussian fit velocity field" echo " 5. Clean up" echo "" echo -n "? " set start = $< endif switch ($start) case 2: goto deconv case 3: goto domom case 4: goto dogauss case 5: goto dorm endsw # Output some interesting information echo "" echo "ccombine.csh $version executed by $user on `date`" | tee $log echo "" | tee -a $log #-------------------------------------------------------------------------- # 1. Make dirty map: #-------------------------------------------------------------------------- # Make the map rm -rf $source.map $source.beam echo "Included files: $files" | tee -a $log echo "Making $nmaps maps beginning at $vfirst and spaced by $delv km/s" \ | tee -a $log echo "" | tee -a $log invert vis="$files" map=$source.map beam=$source.beam cell=$cell \ imsize=$imsize line=velocity,$nmaps,$vfirst,$delv,$delv \ options=systemp,double,mosaic slop=1 robust=$robust $exinvert # Calculate rms # set nmin1 = `expr $nmaps - 1` # echo "Dirty map, Channels 1-2, $nmin1-$nmaps, Mean:" | tee -a $log # set exec="histo in=$source.map" # set box="box(-60,-60,60,60)" # set crmsi=`$exec region="a,$box(1,2)" |grep Rms |awk '{print $4}'` # set crmsf=`$exec region="a,$box($nmin1,$nmaps)" |grep Rms |awk '{print $4}'` # set crms = `calc -f f7.5 "($crmsi+$crmsf)/2"` # echo "RMS in central 2 x 2 arcm: $crmsi $crmsf Mean: $crms" | tee -a $log # set frmsi=`$exec region="image(1,2)" | grep Rms | awk '{print $4}'` # set frmsf=`$exec region="image($nmin1,$nmaps)" | grep Rms | awk '{print $4}'` # set frms = `calc -f f7.5 "($frmsi+$frmsf)/2"` # echo "RMS in entire field: $frmsi $frmsf Mean: $frms" | tee -a $log # Get the (flattened) noise level. rm -rf $source.sen $source.gain mossen in=$source.map sen=$source.sen gain=$source.gain rm -rf $source.mapnorm set minval=`histo in=$source.sen region='image(1)'|grep Min|awk '{print $3}'` echo "Normalizing $source.sen by $minval" | tee -a $log maths exp="$source.map*$minval/$source.sen" out=$source.mapnorm set clip=`calc "$maxsen*$minval"` rm -rf $source.senmsk maths exp=$source.sen out=$source.senmsk \ mask="<$source.sen>/<$source.gain>.lt.$clip" set nmin1 = `expr $nmaps - 1` set exec="sigest in=$source.mapnorm" set drmsi=`$exec region="image(1,2)" |grep Est|awk '{print $4}'` set drmsf=`$exec region="image($nmin1,$nmaps)" |grep Est|awk '{print $4}'` set drms = `calc -f f7.5 "($drmsi+$drmsf)/2"` echo "Estimated RMS noise, dirty map: $drms" | tee -a $log set scale=`calc -f f7.3 "$drms/$minval"` echo "Scaling factor for rmsfac: $scale" | tee -a $log puthd in=$source.mapnorm/rms value=$drms type=real puthd in=$source.mapnorm/scale value=$scale type=real echo " " | tee -a $log # Make the plot. cgdisp in=$source.mapnorm device=$source.dmaps.ps/ps \ region='arc,box(-90,-90,90,90)' \ options=full,3value labtyp=arcsec nxy=4,3 type=contour \ slev=a,$drms levs1=2,3,4,5,6,8,10,12,14 csize=0.6,0.8 if ($showplot != 'n') then gv $source.dmaps.ps & endif if ($interact != 'n') then echo -n "Continue with deconvolution [n]? " set resp=$< if ($resp != 'y') goto end endif #-------------------------------------------------------------------------- # 2. Clean the map: #-------------------------------------------------------------------------- deconv: # Deconvolve using MOSSDI set scale=`gethd in=$source.mapnorm/scale` set cutoff=`calc "$scale*$rmsfac"` echo "${decon}: rmsfac=$rmsfac, cutoff=$cutoff" | tee -a $log rm -rf $source.psf mospsf beam=$source.beam out=$source.psf if ($bmaj == '') then imfit in=$source.psf object=beam region="a,box(-10,-10,10,10)" \ | tee tmp.imfit set bmaj=`grep 'Major axis' tmp.imfit | sed 's/.*://' | awk '{print $1}'` set bmin=`grep 'Minor axis' tmp.imfit | sed 's/.*://' | awk '{print $1}'` set bpa=`grep 'Position angle' tmp.imfit | sed 's/.*://' | awk '{print $1}'` rm -f tmp.imfit else if ($bmin == '') set bmin=$bmaj if ($bpa == '') set bpa=0 endif rm -rf $source.cln if ($decon == 'mossdi2') then if ($clnbox == '') then mossdi2 map=$source.map beam=$source.beam out=$source.cln \ cutoff=$cutoff niters=20000 region="mask($source.senmsk)" else mossdi2 map=$source.map beam=$source.beam out=$source.cln \ cutoff=$cutoff niters=20000 region=$clnbox endif endif if ($decon == 'mosmem') then rm -rf $sdmap.rg regrid in=$sdmap out=$sdmap.rg tin=$source.map if ($clnbox == '') then mosmem map=$source.map,$sdmap.rg beam=$source.beam,$sdbeam \ out=$source.cln rmsfac=$scale \ niters=20000 region="mask($source.senmsk)" $exmem else mosmem map=$source.map,$sdmap.rg beam=$source.beam,$sdbeam \ out=$source.cln rmsfac=$scale \ niters=20000 region=$clnbox $exmem endif endif # Produced cleaned maps and noise-flattened image cube rm -rf $source.cm restor map=$source.map beam=$source.beam model=$source.cln out=$source.cm \ fwhm=$bmaj,$bmin pa=$bpa copyhd in=$source.map out=$source.cm items=mask rm -rf $source.cmnorm maths exp="$source.cm*$minval/$source.sen" out=$source.cmnorm set nmin1 = `expr $nmaps - 1` set exec="sigest in=$source.cmnorm" set crmsi=`$exec region="image(1,2)" |grep Est|awk '{print $4}'` set crmsf=`$exec region="image($nmin1,$nmaps)" |grep Est|awk '{print $4}'` set crms = `calc -f f7.5 "($crmsi+$crmsf)/2"` echo "Estimated RMS noise, cleaned map: $crms" | tee -a $log puthd in=$source.cmnorm/rms value=$crms type=real # Produce gain-corrected image cube and noise cube set minval=`histo in=$source.sen region='image(1)'|grep Min|awk '{print $3}'` set clip=`calc "$msklev*$minval"` rm -rf $source.nse maths exp="<$source.sen>*$crms/(<$source.gain>*$minval)" out=$source.nse \ mask="<$source.sen>/<$source.gain>.lt.$clip" rm -rf $source.cmmsk maths exp="<$source.cm>/<$source.gain>" region="mask($source.nse)" \ out=$source.cmmsk rm -rf $source.cmnse maths exp="<$source.nse>" region="mask($source.nse)" out=$source.cmnse # Produced noise-flattened residual cube rm -rf $source.res restor map=$source.map beam=$source.beam model=$source.cln out=$source.res \ fwhm=$bmaj,$bmin pa=$bpa mode=residual rm -rf $source.resnorm maths exp="$source.res*$minval/$source.sen" out=$source.resnorm copyhd in=$source.cm out=$source.resnorm items=bmaj,bmin,bpa set rflux=`histo in=$source.resnorm | grep Flux | awk '{print $6}'` set rrms=`sigest in=$source.resnorm | grep Est | awk -v c=$crms '{print $4/c}'` set rintflux=`calc -f f8.2 "$rflux*$delv"` set cflux=`histo in=$source.cln | grep Flux | awk '{print $6}'` set cintflux=`calc -f f8.2 "$cflux*$delv"` echo "Used restoring beam of $bmaj x $bmin arcsec, pa $bpa" | tee -a $log echo "Flux in residuals map: $rintflux Jy km/s" | tee -a $log echo "Flux in deconv. model: $cintflux Jy km/s" | tee -a $log echo "RMS of residual in sigma units: $rrms" | tee -a $log echo " " | tee -a $log set maxval=`histo in=$source.cmnorm | grep Max | awk '{print $3}'` cgdisp in=$source.cmnorm,$source.cmnorm device=$source.cmaps.ps/ps \ region='arc,box(-90,-90,90,90)' options=full,3value,mirr \ type=c,p labtyp=arcsec nxy=4,3 lines=1,1,2 slev=a,$crms \ levs1=2,3,4,5,6,8,10,12,14 csize=0.6,0.8 range=0,$maxval,lin,1 if ($showplot != 'n') then gv $source.cmaps.ps & endif if ($interact != 'n') then echo -n "Continue with moment maps [n]? " set resp=$< if ($resp != 'y') goto end endif #-------------------------------------------------------------------------- # 3. Make moment maps: #-------------------------------------------------------------------------- domom: # Make moment maps moment.csh file=$source.cmmsk noisemap=$source.cmnse start=straight \ nspan=$nspan interact=$interact moment.csh file=$source.cmmsk noisemap=$source.cmnse start=mask \ fwhm=$fwhm interact=$interact moment.csh file=$source.cmmsk noisemap=$source.cmnse start=gaufit \ fluxclip=$fluxclip velclip=$velclip moment.csh file=$source.cmmsk start=cleanup #rm -rf $source.senmsk1 # set minval=`histo in=$source.sen region='image(1)'|grep Min|awk '{print $3}'` # echo "Normalizing $source.sen by $minval" | tee -a $log # set clip=`calc "$maxsen*$minval"` # rm -rf $source.senmsk # maths exp=$source.sen out=$source.senmsk mask="<$source.sen>.lt.$clip" # # imsub in=$source.senmsk out=$source.senmsk1 region='image(1)' # foreach mom (dmom0 mmom0 vmom0) # immask in=$source.$mom region="mask($source.senmsk1)" logic=and # immask in=$source.$mom.nse region="mask($source.senmsk1)" logic=and # end # Plot the mom-0 and peak brightness maps if ($mom0fun == 'sqr') then set levs=1,4,9,16,25,36,49,81,100 else set levs=`awk -v del=$loglevs 'BEGIN {for (i=0;i<10;i++) printf "%3.1f,",10**(i*del); printf "%3.1f\n",10**(10*del)}'` endif set maxval=`histo in=$source.mmom0 | grep Max | awk '{print $3}'` set mom0rms=`histo in=$source.dmom0.nse | grep Min | awk '{print $3}'` echo "The rms of the dmom0 image is $mom0rms" if ($mom0lev == '') then set mom0lev=`calc "3*$mom0rms"` endif mv $source.mmom0/mask $source.mmom0/mask1 cgdisp in=$source.mmom0,$source.mmom0 type=p,c \ region=$psbox device=$source.mmom0.ps/cps \ range=0,$maxval,lin,-8 slev=a,$mom0lev levs1=$levs \ labtyp=hms,dms nxy=1,1 options=black,wedge,full cols1=1 csize=1 \ lines=2,1 beamtyp=b,l if ($showplot != 'n') then gv $source.mmom0.ps & endif mv $source.mmom0/mask1 $source.mmom0/mask set crms=`gethd in=$source.cmnorm/rms` set maxval=`histo in=$source.peak | grep Max | awk '{print $3}'` cgdisp in=$source.peak,$source.peak type=p,c \ region=$psbox device=$source.peak.ps/cps \ range=0,$maxval,lin,-8 slev=a,$crms levs1=$levs \ labtyp=hms,dms nxy=1,1 options=black,solneg2,wedge,full cols1=1 csize=1 \ lines=2,1 beamtyp=b,l if ($showplot != 'n') then gv $source.peak.ps & endif # cgdisp in=$source.vmom0,$source.vmom0,$source.nse type=p,c,c \ # region=$psbox device=$source.vmom0.ps/cps \ # range=0,$maxval,lin,-8 slev=a,$mom0lev,a,1 levs1=$levs levs2=2 \ # labtyp=hms,dms nxy=1,1 options=black,solneg2,wedge,full cols1=1 csize=1 \ # lines=2,1,2 beamtyp=b,l # if ($showplot != 'n') then # gv $source.vmom0.ps & # endif # Plot the mom-1 maps set vlast=`calc "$vfirst+($nmaps-1)*$delv"` set nlevs=`calc -i "($vlast-$vfirst)/$mom1del"` awk -v v0=$vfirst -v dv=$mom1del -v nl=$nlevs \ 'BEGIN { for (i = 1; i < nl; i++) {printf "%d,",v0+i*dv} printf "%d\n",v0+nl*dv}' \ > vlevs.txt cgdisp in=$source.vpeak,$source.vpeak type=p,c \ region=$psbox device=$source.vpeak.ps/cps \ range=$vfirst,$vlast,lin,3 slev=a,1 levs1=@vlevs.txt \ labtyp=hms,dms nxy=1,1 options=black,wedge,full cols1=1 csize=1 lines=2,1 if ($showplot != 'n') then gv $source.vpeak.ps & endif cgdisp in=$source.gmom1,$source.gmom1 type=p,c \ region=$psbox device=$source.gvel.ps/cps \ range=$vfirst,$vlast,lin,3 slev=a,1 levs1=@vlevs.txt \ labtyp=hms,dms nxy=1,1 options=black,wedge,full cols1=1 csize=1 lines=2,1 if ($showplot != 'n') then gv $source.gvel.ps & endif rm -f vlevs.txt mv moment.log $source'_moment.log' mv mossdi2_iteration.log $source'_iteration.log' if ($interact != 'n') then echo -n "Continue with cleanup [n]? " set resp=$< if ($resp != 'y') goto end else if ($purge != 'y') then goto end endif #-------------------------------------------------------------------------- # 5. Clean up files #-------------------------------------------------------------------------- echo "-----dorm------" dorm: # Clean up: rm -rf $source.res $source.resnorm $source.psf rm -rf $source.nse $source.senmsk $source.gau* gzip -f $source.*.ps end: