From dshepher@aoc.nrao.edu Mon Mar 5 10:12:13 2001 Date: Fri, 2 Mar 2001 17:46:35 -0700 (MST) From: Debra Shepherd To: Steven Myers Cc: Ed Fomalont , Greg Taylor , John Benson , John Hibbard Subject: aips++ script for qband 3c84 data ---------- Forwarded message ---------- Date: Fri, 2 Mar 2001 17:37:51 -0700 (MST) From: Debra Shepherd To: Tim Cornwell , George Moellenbrock Cc: Debra Shepherd Subject: aips++ script for qband 3c84 data 02mar01 *** This script is not ready for outside users -- still too much embedded testing. - aips++ script finished on 2mar01, most testing comments have been deleted. Times to run each task during calibration and imaging are given. Final comparison with aips-processed data is given at the end of this script. AIPS++ Script produced by Debra Shepherd ================================================================== source /aips++/stable/aipsinit.csh /aips++/stable/linux_gnu/bin/aips++ dowait := T; 22dec00 data -- VLA Q band data for 3c84, no fast-switching Information on the DATA: program number DSTST - image details: VLA A array synth beam Q: 0.05'' K: 0.1'' Primary beam: 60'' 120'' - Q band theory rms approx. 0.11 mJy/beam = 0.00011 Jy/bm - K band theory rms approx. 0.14 mJy/beam = 0.00014 Jy/bm AIPS (without selfcal) gets to - Q: AIPS 1: Mean= 9.2047E-03 rms= 2.3293E-02 JY/BEAM over 3731. pixels - K: AIPS 1: Mean= 6.8843E-04 rms= 2.2115E-02 JY/BEAM over 3731. pixels ** Visibility Tape Information ** Tape # File # Time of Final Record: N13920 1 src = '''' (3c84) gcal = '0319+415''' (3c84) 5s scans, no fast switching, just sit and stare at 3c84 for ~1 hr at Q ******* ******* ******* Q-band continuum in subarray 1 ******* ******* ******* 0319+415 J2000 A 03h19m48.1601s 41d30'42.106" 3C84 0316+413 B1950 A 03h16m29.5673s 41d19'51.919" NGC1275 ----------------------------------------------------- BAND A B C D FLUX(Jy) UVMIN(kL) UVMAX(kL) ===================================================== 90cm P S X X X 8 13 20cm L P P X X 23.9 12 6cm C P P P P 23.3 3.7cm X P P P P 21.70 visplot 2cm U P P P P 20.70 visplot 1.3cm K X S S S 16.4 0.7cm Q X S S S 9.00 1800 visplot This script fills, edits, calibrates, and images. The image is ok but suffers from a few bad visibility points. At the end of this, run visplot and edit on the residual visibilities. # # Fill from tape drive (a linux box) # Manually configure the Linux hardware driver to use variable length blocks. mt -f /dev/nst0 setblk 0 copy DN tape data with unix command onto disk: advance past the first record with housekeeping data. mt -f /dev/nst0 setblk 0 mt -f /dev/tape rewind mt -f /dev/tape fsf 1 dd if=/dev/tape bs=26624 of=22dec00.vla -- 6256+21 records in -- 6256+21 records out include 'vlafiller.g'; myvlafiller:=vlafiller(host=[i_am_unset="i_am_unset" ], forcenewserver=F); ok:=myvlafiller.diskinput(filename="22dec00.vla" ); ok:=myvlafiller.output(msname="q3c84.ms" , overwrite=F); ok:=myvlafiller.selectproject(project="DSTST" ); ok:=myvlafiller.selectfrequency(centerfrequency=[value=43.3, unit="GHz" ], bandwidth=[value=500.0, unit="MHz" ]); ok:=myvlafiller.fill(verbose=T, async=F); # # RESULTS: Writing to measurement set q3c84.ms 2 new entries in the spectral window sub-table Row 0 1 channel. Ref. frequency: 43.315 GHz IF: 0 Row 1 1 channel. Ref. frequency: 43.365 GHz IF: 1 1 new entry in the polarization sub-table Row 0 correlation products RR, LL, RL, LR 1 new entry(s) in the field sub-table Row 0: Field direction (J2000): (03:19:48.160, +041.30.42.106) The measurement set contains 421848 rows. The antenna sub-table contains 27 entries The field sub-table contains 1 entries The polarization sub-table contains 1 entries The spectral window sub-table contains 2 entries myvlafiller.done(); ================================================================================ MeasurementSet Name: /export/home/sola/dss/aips++.tsting/22dec00/q3c84.ms MS Version 2 ================================================================================ Observer: unavailable Project: DSTST Observation: VLA Fields: 1 ID Name Right Ascension Declination Type 1 0319+415 03:19:48.16 +41.30.42.11 J2000 Data records: 421848 Total integration time = 2830 seconds Observed from 22-Dec-2000/03:44:58 to 22-Dec-2000/04:32:08 Spectral Windows: 2 Ref.Freq #Chans Resolution TotalBW 43314.9MHz 1 60500 kHz 50000 kHz 43364.9MHz 1 60500 kHz 50000 kHz Polarization setups: 1 Correlations RR LL RL LR Feeds = 27: printing first row only Antenna Spectral Window # Receptors Polarizations 1 -1 2 [ R, L] Antennas: 27 1-9: VLA:W48 VLA:W40 VLA:W16 VLA:W8 VLA:W24 VLA:W32 VLA:W64 VLA:W72 VLA:W56 10-18: VLA:E48 VLA:E40 VLA:E16 VLA:E8 VLA:E24 VLA:E32 VLA:E64 VLA:E72 VLA:E56 19-27: VLA:N48 VLA:N40 VLA:N8 VLA:N32 VLA:N16 VLA:N24 VLA:N64 VLA:N72 VLA:N56 -------------------------------------------------------------------------------- Tables: MAIN 421848 rows ANTENNA 27 rows DATA_DESCRIPTION 2 rows DOPPLER FEED 27 rows FIELD 1 row FLAG_CMD FREQ_OFFSET HISTORY OBSERVATION 1 row POINTING POLARIZATION 1 row PROCESSOR SOURCE (see FIELD) SPECTRAL_WINDOW 2 rows STATE SYSCAL WEATHER ================================================================================ NOTE: 27 ants, no antennas were deleted from the subarray but there should be no data for a number of antennas. - Get an antenna plot: include 'visplot.g'; visplot:=visplot(msfile="q3c84.ms" , nrows=16384, edit=T, flagfile=[i_am_unset="i_am_unset" ], displaywidth=600, displayheight=350); select region "antenna", click on spanner and then "from MS" the Y config plot appears, choose file - print. labels are jumbled in the center of the array, need to see antenna numbers too. include 'flagger.g'; flag:=flagger(msfile="q3c84.ms" ); # flag autocorrelations for easier plotting (can't use them here). ok:=flag.auto(trial=F); # Flagging 30132 rows flag.done(); # # Set the flux for 3c84 -- the only source observed here. # include 'imagerwizard.g'; myimager:=imager(filename="q3c84.ms" ); ok:=myimager.setjy(fieldid=1, spwid=-1, fluxdensity=9.0); # 3c84 (expect 9 Jy at Q band): - 0319+415 spwid= 1 [I=9, Q=0, U=0, V=0] Jy, (user-specified) - 0319+415 spwid= 2 [I=9, Q=0, U=0, V=0] Jy, (user-specified) - timing: This takes approximately 30s for each spwid myimager.done(); # # Solve for the gains and apply # include 'calibrater.g'; calib:=calibrater(filename="q3c84.ms" ); ok:=calib.setdata(mode="none" , nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], uvrange=0, msselect=''); ok:=calib.setsolve(type="G" , t=5, preavg=0.0, phaseonly=F, refant=21, table="g3c84.gcal2" , append=F); ######### choose refant=21 (aips++ antenna index 21=vla ant 8) to be identical to VLA obs. include 'calibrater.g'; calib:=calibrater(filename="q3c84.ms" ); ok:=calib.setdata(mode="none" , nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], uvrange=0, msselect=''); ok:=calib.setsolve(type="G" , t=5, preavg=0.0, phaseonly=F, refant=21, table="g3c84.gcal2" , append=F); - timing: For interval of 5 seconds, found 1116 slots Solving for G G Jones Slot=1, 0319+415, spw=1: 22-Dec-2000/03:44:55 to 22-Dec-2000/03:44:59 G Jones Sum of weights is zero G Jones Slot=2, 0319+415, spw=1: 22-Dec-2000/03:45:00 to 22-Dec-2000/03:45:04 G Jones Sum of weights is zero G Jones Slot=3, 0319+415, spw=1: 22-Dec-2000/03:45:05 to 22-Dec-2000/03:45:09 G Jones Initial fit per unit weight = 0.528916 Jy, sum of weights = 376 G Jones Final fit per unit weight = 0.0297041 Jy after 8 iterations G Jones Slot=4, 0319+415, spw=1: 22-Dec-2000/03:45:10 to 22-Dec-2000/03:45:14 G Jones Initial fit per unit weight = 0.526907 Jy, sum of weights = 454 G Jones Final fit per unit weight = 0.0298886 Jy after 8 iterations G Jones Slot=5, 0319+415, spw=1: 22-Dec-2000/03:45:15 to 22-Dec-2000/03:45:19 G Jones Initial fit per unit weight = 0.527103 Jy, sum of weights = 456 G Jones Final fit per unit weight = 0.0307524 Jy after 8 iterations ... 8 iterations/slot, some slots have zero weight (flagged) each slot takes 3-4s after 67 minutes, 1116 slots processed: Finished calibrater::solve 4025.53 real 3720.57 user 41.56 system - Comparison: AIPS took 2 minutes 4 seconds to process the same data - with exactly the same inputs. AIPS CPU time was only 51.2 seconds: - CALIB1: Found 43500 good solutions - CALIB1: 8 solutions had insufficient data - CALIB1: Appears to have ended successfully - CALIB1: sola 31DEC01 TST: Cpu= 51.2 Real= 124 --- AIPS++ solver took 32 times longer than AIPS. # Don't need to bootstrap flux density scale from flux cal to gcal # because I only have 3c84 obs here. # Apply the gcal: # select all source uv data to be corrected: ok:=calib.setdata(mode="none" , nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], uvrange=0, msselect=''); ok:=calib.setapply(type="G" , t=0.0, table="g3c84.gcal2" , select=''); ok:=calib.correct(); - timing: Finished calibrater::correct 146.93 real 131.76 user 6.28 system Time to correct data = 2.4 minutes. # check the calibration... ok:=calib.plotcal(plottype="AMP" , tablename="g3c84.gcal2" , antennas=[], polarization=1, spwids=[], timeslot=1); calib.done(); --------------------------------------------------------------- --------------------------------------------------------------- Make an Image --------------------------------------------------------------- --------------------------------------------------------------- include 'mosaicwizard.g'; im:=imager(filename="q3c84.ms" ); - timing: took 33 seconds just to construct image tool ok:=im.advise(takeadvice=F, amplitudeloss=0.05, fieldofview=[value=1.0, unit="arcmin" ], pixels=3200, cell=[value=0.0195129378, unit="arcsec" ], facets=1, phasecenter=[type="direction" , refer="J2000" , m1=[value=0.724515775, unit="rad" ], m0=[unit="rad" , value=0.871803604]]); - aips++ recommended J2000: 03:19:48.160 +041.30.42.106 - with aips I used cell=0.01 and imsize=1024 - timing: Finished imager::advise 9.86 real 7.39 user 1.36 system # define what data can be imaged: ok:=im.setdata(mode="none" , nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], spwid=[1, 2] , fieldid=1, msselect=''); - timing: Finished imager::setdata 5.52 real 3.71 user 1.52 system # define image properties: ok:=im.setimage(nx=1024, ny=1024, cellx=[value=0.01, unit="arcsec" ], celly=[value=0.01, unit="arcsec" ], stokes="I" , doshift=F, phasecenter=[m0=[value=49.9506667, unit="deg" ], m1=[value=41.5116961, unit="deg" ], type="direction" , refer="J2000" ], shiftx=[value=0.0, unit="arcsec" ], shifty=[value=0.0, unit="arcsec" ], mode="mfs" , nchan=1, start=1, step=1, mstart=[value=0.0, unit="km/s" ], mstep=[value=0.0, unit="km/s" ], spwid=[1, 2] , fieldid=1, facets=1); -timing: no statistics given, took less than a minute. # with aips, I used natural weighting. Do so here. ok:=myimager.weight(type="natural" , rmode="norm" , noise=[value=0.0, unit="Jy" ], robust=0.0, fieldofview=[value=0.0, unit="rad" ], npixels=0); - timing: Finished imager::weight 14.06 real 10.97 user 1.53 system make a DIRTY image... ok:=im.makeimage(type="corrected" , image="q3c84.map2" , compleximage='' ); - Calculating image (without full skyequation) - Making dirty image from corrected data - Image is : q3c84.map - Retaining complex image: q3c84.cmplxmap - Frequency = 43.3149 GHz, synthesized continuum bandwidth = 0.1105 GHz - Image polarization = Stokes I - Frequency = 43.3149 GHz, synthesized continuum bandwidth = 0.1105 GHz - Image polarization = Stokes I - Fourier transforms will use image centers as tangent points - timing: Finished imager::makeimage 47.12 real 25.93 user 5.06 system view DIRTY image statistics: max = 9.023 Jy at position 513 513 1 03 19 48.160 +41 30 42.106 mean = -1.551e-4 rms = 2.061e-1 make a regionmask for clean from the dirty image: - From the wrench select "From image": you'll be asked for an image to display. Select that as usual. - On the resulting viewer, make a region and doubleclick. - It will then show up in the inputs for regionmask. - Set a name for the mask and press go. When it finishes, you can use the wrench View to verify that the mask is as you expect. ok:=myimager.regionmask(mask="q3c84.mask2" , region="myregion" ... - (scripter broke so can't give exact syntax here -- still broke as of 02mar01). - Run imager.clean with mask set to q3c84.mask. use same weight function 'na', same flux cutoff (0.00011), and same niter = 100000 ok:=im.clean(algorithm="clark" , niter=100000, gain=0.1, threshold=[value=0.00011, unit="Jy" ], displayprogress=T, model="q3c84.clean2" , fixed=F, complist='', mask="q3c84.mask2" , image="q3c84.cm2" , residual="q3c84.clean.resid2" ); - Iteration: 94783, Maximum residual=0.000109911 - Max Resid = -0.000108685 - Iteration 94784 flux [Jy] = 9.19961 - Clean used 94783 iterations to get to a max residual of 0.000109911 - Successfully deconvolved image - timing: Finished imager::clean 843.46 real 130.72 user 25.01 system so aips++ took 14 minutes to clean this field (the field was cleaned much deeper than it had to be). - restore the clean image. ok:=im.restore(model="q3c84.clean2" , complist='', image="q3c84.cm" , residual="q3c84.resid" ); - timing: Finished imager::restore 119.58 real 61.23 user 13.84 system ======================================================================= ======================================================================= view CLEAN aips++ image: max = 8.982 Jy at position 512 513 1 03 19 48.178 +41 30 41.986 mean = -1.464e04, min = -2.278e-2, max = 3.664e-2 rms = 6.407e-3 compare with DIRTY aips++ image: statistics: max = 9.022 Jy at position 513 513 1 03 19 48.160 +41 30 42.106 mean = 3.01e-3 rms = 2.25e-1 compare with aips: max = 8.781 Jy at position 512 513 1 mean = 9.20e-3 rms = 2.33e-2 rms seems to be better than image processed with aips, but there are obvious sidelobes and stripes that aips image doesn't have. bring up a viewer tool - contours not available unless the tool is brought up explicitely. - stats not available unless image is brought up via wrench or double click... make a contour plot and raster image. The rms in the image does seem to be better than aips but the stats do not seem to be correct: - contour levels plotted at 0.023 x -3,3,4,5,... (like aips) gives only central point, nothing in rest of map - contour levels plotted at 0.0079 x -3,3,4,5,... (rms measured in aips++ image) gives far too many contours in surrounding map and shows low-level, N-S, stripes. Either I have some bad antennas that I'm not able to find, or the stats are not working properly. Also, there is no weighting in this aips++ dataset and this could also be the reason why I see low level stripes.