Date: Tue, 14 Jan 2003 16:27:59 -0700 (MST) From: Debra Shepherd To: Steven Myers Subject: new script to reduce gg tau Every thing works now. I got out a decent, consistent, 3mm continuum image from 07feb97 with this script. I think I'm ready to start drafting some sort of cookbook and then continue testing. We need to try the 10feb97 data because I understand that the phcorr will actually do something on it (data is worse than the 07feb data). We also need to reduce the 1mm data now that Kumar will have this working in tomorrow's daily. Deb ======================================================== Sorry for this being so verbose. Info needed for data reduction: ========================================================================= pbcal = 0528+134 phase and amplitude cals for 7feb97: 0528+134 0415+379 CRL618 ========================================================================= The first column is the flux at 3mm and the other is at 1mm filename="flux.txt" 0528 0415 CRL618 MWC349 07-feb-1997 2.63 2.05 5.81 4.63 1.55 2.20* 10-feb-1997 2.92 1.80 6.25 4.48 0.96 1.66* 20-feb-1997a 2.45 1.50* 5.43 3.56 0.96 1.66* 20-feb-1997b 2.41* 1.53* 5.38 3.96 1.60 2.50 25-mar-1997 2.30 1.54* 3.80* 2.23* 31-mar-1997 2.17 1.54 3.68 2.23 1.74 3.64 0.96* 1.66 18-sep-1997 1.98 1.25* 2.80* 2.07 1.95 2.99 16-oct-1997 1.92 1.09* 2.83 2.13 1.66 2.13 0.96 1.66* 18-oct-1997 1.98 1.17 2.87 2.03 0.96 1.66* ========================================================================= Copy GGTau FITS-IDI alma format: There are multiple files for each of 8 days of observation, spread over many months in 1997 & different iram configs. mkdir calfits cd calfits cp /home/maser2/alma-ti/cdrom-4/calfits/*.fits.gz . cp /home/maser2/alma-ti/cdrom-5/calfits/*.fits.gz . gunzip *.fits.gz Use the almatifiller global function to fill data: # script to fill g067 data (8 days) include 'logger.g' include 'os.g' include 'almati2ms.g'; # Fill data from CDROM-4 almatifiller(msfile='07feb97-g067.ms', fitsdir='./calfits', pattern='07-feb-1997*', append=F, compress=F, obsmode="CORR", chanzero="TIME_AVG"); Logger: Filling data to the output MS Processing primary HDU... Successfully closed empty server: almati2ms glish: ... filename= 07-feb-1997-g067-16.fits fullname= ./calfits/07-feb-1997-g067-16.fits F ***** BUG(): The filler should have put the IRAM phase-corrected data in the DATA column but it didn't. It put the non-phase-corrected data in the DATA column. Kumar knows about this and will fix it soon. In the mean time he has a script that will put the proper phase-corrected data in the DATA column. t:=table('07feb97-g067.ms', readonly=F); for (k in 1:t.nrows()){ dat:=t.getcell('ALMA_PHAS_CORR', k); t.putcell('DATA',k,dat); } t.done(); OK, things are as they should be: DATA = CORRECTED_DATA = ALMA_PHAS_CORR The ALMA_NO_PHAS_CORR row is not equal to the above and all checks on ALMA_PHAS_CORR_FLAG_ROW = F (plot (DATA_DESC_ID vs ALMA_PHAS_CORR_FLAG_ROW) ####################################################################### Work with 7feb97 data only for now using iramcalibrater GUI. #################################### Obtain a summary of the data: #################################### include 'ms.g'; myms:=ms(filename="07feb97-g067.ms" , readonly=T, lock=F); myms.summary(verbose=T); # summary saved in 07feb97.summary #################################### Start IRAM calibrator tool: #################################### include 'iramcalibrater.g'; ical:=iramcalibrater(msname="07feb97-g067.ms" ); # Execution will construct new tool ical # Starting server calibrater # Server started: /aips++/daily/linux_gnu/bin/calibrater # (AIPS++ version: 1.8 (build #431)) # Adding MODEL_DATA, CORRECTED_DATA and IMAGING_WEIGHT columns # Initializing MODEL_DATA (to unity) and CORRECTED_DATA (to DATA) #################################### Use CLIC quality metric to select on-line phase-corrected or uncorrected data: #################################### By default the DATA column is the phase-corrected data. However, if the on-line phase-corrected data is not good, then it is better to use the non-phae-corrected data. ical.phcor will look at the data, use the CLIC quality metric to determine whether it is good, if it is not then the non-phase corrected data will be selected. Phcor only works in daily right now: ok:=ical.phcor(trial=T); ok:=ical.phcor(trial=F); -----Begining method: phcor successful readonly open of default-locked table 07feb97-g067.ms/SOURCE: 14 columns, 4 rows successful readonly open of default-locked table 07feb97-g067.ms/FIELD: 9 columns, 4 rows successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows successful readonly open of default-locked table 07feb97-g067.ms/ANTENNA: 8 columns, 5 rows successful read/write open of default-locked table 07feb97-g067.ms: 29 columns, 80400 rows successful readonly open of default-locked table 07feb97-g067.ms/SOURCE: 14 columns, 4 rows successful readonly open of default-locked table 07feb97-g067.ms/FIELD: 9 columns, 4 rows successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows successful readonly open of default-locked table 07feb97-g067.ms/ANTENNA: 8 columns, 5 rows successful read/write open of default-locked table 07feb97-g067.ms: 29 columns, 80400 rows successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows -----Finished method: phcor Kumar says that the 07feb97 data are just fine. No data will be flagged. OK, there was no mention of flagging data here. proceed. #################################### Perform a bandpass calibration: #################################### ok:=ical.rf(fieldname="0528+134" , freqgrp="3mm-LSB" , visnorm=T, bpnorm=T, gibb=2, drop=0); -----Begining method: rf Initializing MODEL_DATA (to unity) and CORRECTED_DATA (to DATA) successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows The following calibration components will be applied: None. The following calibration components will be solved for: None. Selecting data Performing selection on MeasurementSet By selection 80400 rows are reduced to 660 The following calibration components will be solved for: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 preavg=0 phaseonly=F refant=1 append=F The following calibration components will be applied: None. The following calibration components will be solved for: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 preavg=0 phaseonly=F refant=1 append=F Initializing solvable bandpass terms (B-matrix) Solving for B Freq. group 3mm-LSB, spw= 2, nchan= 64 Freq. group 3mm-LSB, spw= 6, nchan= 256 -----Finished method: rf # two plots are produced and written to the current directory (no control over plotting parameters): 07feb97-g067.ms.RF_AMP.ps 07feb97-g067.ms.RF_PHASE.ps # a table with the bcal solutions is also written out: 07feb97-g067.ms.3mm-LSB.bcal The bcal table has 5 rows which specify the polynomial fits derived above. Not sure why there are only 5 rows (5 ants?) because the plots show solutions for 10 BASELINES, not 5 ANTS: # DATA = CORRECTED_DATA, solutions are not applied at this point. #################################### Phase calibration: #################################### ok:=ical.phase(fieldnames='0528+134 0415+379 CRL618', freqgrp="3mm-LSB" ); -----Beginning method: phase successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows The following calibration components will be applied: None. The following calibration components will be solved for: None. Selecting data Performing selection on MeasurementSet By selection 80400 rows are reduced to 1580 The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] The following calibration components will be solved for: GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 preavg=0 phaseonly=F refant=1 append=F The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] The following calibration components will be solved for: GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 preavg=0 phaseonly=F refant=1 append=F Applying BPOLY table from 07feb97-g067.ms.3mm-LSB.bcal Initializing solvable electronic gain terms (G-matrix) Solving for G Slot= 1, field= 0415+379, spw= 2, nchan= 64 Slot= 2, field= 0415+379, spw= 6, nchan= 256 Slot= 3, field= 0528+134, spw= 2, nchan= 64 Slot= 4, field= 0528+134, spw= 6, nchan= 256 Slot= 5, field= CRL618 , spw= 2, nchan= 64 Slot= 6, field= CRL618 , spw= 6, nchan= 256 The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 select=[ ] Selecting data Performing selection on MeasurementSet By selection 80400 rows are reduced to 6700 The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 select=[ ] Applying BPOLY table from 07feb97-g067.ms.3mm-LSB.bcal Applying GSPLINE table from 07feb97-g067.ms.3mm-LSB.phcal -----Finished method: phase # Phase solutions written out to a table in working directory: 07feb97-g067.ms.3mm-LSB.phcal/ # Plot of phase solutions and fit to data written out: 07feb97-g067.ms.PHAS.ps DATA does not equal CORRECTED_DATA column. So ical.phase modifies the corrected data. #################################### Establish absolute flux density scale: #################################### Calibrate the calibrators using CRL618 with the modeled flux density of 1.55. ok:=ical.flux(fieldnames='0528+134 0415+379 CRL618', freqgrp="3mm-LSB" , plot=F, gibb=0, drop=0); Glish interaction: - Number of antennas: 5 Give the field names and fluxes of fixed calibrators Field names = CRL618 Fluxes (e.g 1mJy) = 1.55 Logger messages: -----Beginning method: flux successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows successful readonly open of default-locked table 07feb97-g067.ms/ANTENNA: 8 columns, 5 rows Starting server ms Server started: /aips++/daily/linux_gnu/bin/ms (AIPS++ version: 1.8 (build #435)) Channel selection: #chan=64, start=1, width=1, incr=1 Polarization selection: [XX] Selection initialized ok Factor to multiply is 21.1016612 Flux found: 0528+134 = 2.48977771 Jy Flux found: 0415+379 = 5.5148381 Jy Flux found: CRL618 = 1.55 Jy Efficiency of antenna 1 is 21.0103899 +/- 3.02868542 Efficiency of antenna 2 is 20.8148687 +/- 3.02418419 Efficiency of antenna 3 is 21.0664088 +/- 3.00062946 Efficiency of antenna 4 is 20.3869189 +/- 2.95038559 Efficiency of antenna 5 is 20.7545524 +/- 2.95541108 # OK, flux of 0528 should be about 2.63, it is found to be 2.48 Jy 0415 should be about 5.81, it is found to be 5.82 Jy The efficiencies of the antennas at 3mm should be between 20 & 24 Jy/K and they are all around 22 Jy/K above. Solution looks good. Glish interaction: - Apply (y/n) y Selected Logger messages: Starting imager::setjy Fourier transforming: replacing MODEL_DATA column Performing interferometric gridding with convolution function SF Single facet Fourier transforms will use image center as tangent points Using single-field SkyEquation for single image 0528+134 spwid= 1 [I=2.49, Q=0, U=0, V=0] Jy, (user-specified) 0415+379 spwid= 1 [I=5.515, Q=0, U=0, V=0] Jy, (user-specified) CRL618 spwid= 1 [I=1.55, Q=0, U=0, V=0] Jy, (user-specified) 0528+134 spwid= 3 [I=2.49, Q=0, U=0, V=0] Jy, (user-specified) 0415+379 spwid= 3 [I=5.515, Q=0, U=0, V=0] Jy, (user-specified) CRL618 spwid= 3 [I=1.55, Q=0, U=0, V=0] Jy, (user-specified) 0528+134 spwid= 5 [I=2.49, Q=0, U=0, V=0] Jy, (user-specified) 0415+379 spwid= 5 [I=5.515, Q=0, U=0, V=0] Jy, (user-specified) CRL618 spwid= 5 [I=1.55, Q=0, U=0, V=0] Jy, (user-specified) 0528+134 spwid= 7 [I=2.49, Q=0, U=0, V=0] Jy, (user-specified) 0415+379 spwid= 7 [I=5.515, Q=0, U=0, V=0] Jy, (user-specified) CRL618 spwid= 7 [I=1.55, Q=0, U=0, V=0] Jy, (user-specified) All 3mm-LSB flux density values the MODEL_DATA column for the calibrators were set to the derived or specified fluxes. No table is written out and no ps file is generated (I chose plot=F here) #################################### Perform amplitude calibration: #################################### Solve for variations of amplitude with time. Use the flux and gain calibrators to derive the polynomial fits. Apply the corrections to all sources (including the target source). ok:=ical.amp(fieldnames='0528+134 0415+379 CRL618', freqgrp="3mm-LSB"); -----Beginning method: amp successful readonly open of default-locked table 07feb97-g067.ms/SPECTRAL_WINDOW: 14 columns, 24 rows The following calibration components will be applied: None. The following calibration components will be solved for: None. Selecting data Performing selection on MeasurementSet By selection 80400 rows are reduced to 1580 The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 select=[ ] Calibration table 07feb97-g067.ms.3mm-LSB.phcal exists, and append=F. Therefore, this table will be updated. The following calibration components will be solved for: GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 preavg=0 phaseonly=F refant=1 append=F The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 select=[ ] The following calibration components will be solved for: GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 preavg=0 phaseonly=F refant=1 append=F Applying BPOLY table from 07feb97-g067.ms.3mm-LSB.bcal Applying GSPLINE table from 07feb97-g067.ms.3mm-LSB.phcal Initializing solvable electronic gain terms (G-matrix) Solving for G Slot= 1, field= 0415+379, spw= 2, nchan= 64 Slot= 2, field= 0415+379, spw= 6, nchan= 256 Slot= 3, field= 0528+134, spw= 2, nchan= 64 Slot= 4, field= 0528+134, spw= 6, nchan= 256 Slot= 5, field= CRL618 , spw= 2, nchan= 64 Slot= 6, field= CRL618 , spw= 6, nchan= 256 Selecting data Performing selection on MeasurementSet By selection 80400 rows are reduced to 6700 The following calibration components will be applied: BPOLY table=07feb97-g067.ms.3mm-LSB.bcal t=0 select=[ ] GSPLINE table=07feb97-g067.ms.3mm-LSB.phcal t=0 select=[ ] Applying BPOLY table from 07feb97-g067.ms.3mm-LSB.bcal Applying GSPLINE table from 07feb97-g067.ms.3mm-LSB.phcal -----Finished method: amp Plot written to local directory: 07feb97-g067.ms.AMP.ps This is a plot of the residual amplitude solution. The absolute flux density differs slightly from the alma report generated by Athol and Robert - not sure why. The source data should be calibrated at this point. proceed. #################################### use 'uvt' function to write out calibrated source data: #################################### ok:=ical.uvt(fieldname="GG_TAU" , spwid=3, filename="ggtau64cal" , option="new" ); -----Beginning method: uvt Starting server ms Server started: /aips++/weekly/linux_gnu/bin/ms (AIPS++ version: 1.8 (build #433)) Selecting on field and spectral window ids By selection 80400 rows to be converted are reduced to 2560 Converting MeasurementSet tab24831_0 to FITS file 'scratch.fits' MS tab24831_0 is a subset of another MS Writing CORRECTED_DATA column Writing AIPS FQ table Found 1 spectral windows Writing AIPS AN table Found 5 antennas in array #1 Writing AIPS SU table Not setting velocity types writing 1 sources Successfully closed empty server: ms Starting server ms Server started: /aips++/weekly/linux_gnu/bin/ms (AIPS++ version: 1.8 (build #433)) Converting FITS file 'scratch.fits' to MeasurementSet 'ggtau64cal' Using tile shape [1, 7, 2340] for IRAM PdB with obstype=0 Reading and writing 2560 visibility groups Found binary table of type AIPS FQ following data Found binary table of type AIPS AN following data Found binary table of type AIPS SU following data Flushing MS to disk Successfully closed empty server: ms successful read/write open of default-locked table ggtau64cal: 23 columns, 2560 rows -----Finished method: uvt DONE ical and proceed to imaging to see what I've got... ========================================================================= Examination of the Table: Row uvw(J2000) = data flag = 1, 256 boolean flag_category = dummy weight = number sigma = number antenna1 = 0-3 antenna2 = 1-4 array_id = 0 data_desc_id = 0-23 - this is the spwid identifier exposure(s) = 30 or 60 feed1 = 0 feed2 = 0 field_id = 1-4 interval(s) = 60 observation_id = 0 processor_id = -1 scan_number = number like 853 state_id = -1 time(TAI) = number like 1997/02/07/16:38:42 time_centroid(TAI) = number like 1997/02/07/16:38:42 (same as above) data = [1,256]complex weight_spectrum = [1,256]float alma_phase_corr = [1,256]complex alma_no_phas_corr = [1,256]complex alma_phas_corr_flag_row = F (phcor did not modify flags in this data) model_data = [1,256]complex corrected_data = [1,256]complex imaging_weight = [1,256]float #################################### Examine data: #################################### include 'msplot.g'; msplt:=msplot(msfile="07feb97-g067.ms" , edit=T, flagfile=[i_am_unset="i_am_unset" ], nrows=[i_am_unset="i_am_unset" ], displaywidth=600, displayheight=350); msplt.done(); Amp = 0-10 Jy UV dist = 80 - 400 m * there are lots of bad data in here * Look at the calibrators in orginal dataset: include 'msplot.g'; msplt:=msplot(msfile="07feb97-g067.ms" , edit=T, flagfile=[i_am_unset="i_am_unset" ], nrows=[i_am_unset="i_am_unset" ], displaywidth=600, displayheight=350); msplt.done(); field 1 = 0415+379, 64 chans, spwid = 3 plot = msplot.0415.cal.ps flux = about 5.81 Jy but ranges from near 0 to 14 Jy field 2 = 0528+134, 64 chans, spwid = 3 plot = msplot.0528.cal.ps flux = about 2.63 Jy but ranges from 0 to 12 Jy field 4 = CRL618, 64 chans, spwid = 3 plot = msplot.crl618.cal.ps flux = about 1.55 Jy but ranges from 0 to 6.5 Jy The calibration was applied but there is bad data here. #################################### Flag central channels (as per alma script): #################################### Flag the 2 central Gibbs channels prior to imaging: include 'flagger.g'; flg:=flagger(msfile="ggtau64cal" ); MeasurementSet Name: ggtau64cal MS Version 2 Observer: Project: Observation: IRAM PdB(5 antennas) Data records: 2560 Total integration time = 23845 seconds Observed from 07-Feb-1997/16:38:42 to 07-Feb-1997/23:16:07 Fields: 1 ID Name Right Ascension Declination Epoch 1 GG_TAU 04:32:30.34 +17.31.40.52 J2000 Data descriptions: 1 (1 spectral windows and 1 polarization setups) ID Ref.Freq #Chans Resolution TotalBW Correlations 1 89189.8MHz 64 2500 kHz 160000 kHz XX Antennas: 5 ID= 1-5: 1=1 ok:=flg.setchan(chan=[32, 33] ); ok:=flg.query(query='TIME > 0', trial=F); # Flagging 2560 rows OK. Flagging done on each row. flg.done(); include 'msplot.g'; msplt:=msplot(msfile="ggtau64cal" , edit=T); center channels are flagged. Amp = 0-10 Jy UV dist = 80 - 400 m * I don't see any obvious change * all baseline appear to have high values at all times * end channels are causing the problem. flag channels 1,2 and 63,64 in msplot. print before and after data plots: ggtau.before.msplt.ps ggtau.after.msplt.ps msplt.done(); #################################### 3mm continuum imaging: #################################### include 'mosaicwizard.g'; imgr:=imager(filename="ggtau64cal" , compress=F, host='', forcenewserver=T); Logger: Opening MeasurementSet ggtau64cal Adding MODEL_DATA, CORRECTED_DATA and IMAGING_WEIGHT columns Initializing MODEL_DATA (to unity) and CORRECTED_DATA (to DATA) Initializing natural weights imgr.setdata(mode="channel" , nchan=60, start=3, step=1, spwid=1, fieldid=1, msselect=''); imgr.setimage(nx=256, ny=256, cellx=[value=0.2, unit="arcsec" ], celly=[value=0.2, unit="arcsec" ], stokes="I" , mode="mfs" , nchan=1, start=3, step=60, spwid=1, fieldid=1, facets=1); imgr.weight(type="briggs", rmode="norm", robust=0.0); imgr.clean(algorithm="clark" , niter=1000, gain=0.1, threshold=[value=0.0, unit="Jy" ], displayprogress=F, model="ggtau1.model" , fixed=F, complist='', mask="ggtau1.mask" , image="ggtau1.cm" , residual="ggtau1.resid" , interactive=T, npercycle=100); Interactive clean: dirty image is pretty bad - put a large poly around the center. src is at 17deg dec ==> N-S stripe difficult to get out. first 100: Initial maximum residual: 0.00868421 Iteration: 2, Maximum residual=0.00703421 Iteration: 8, Maximum residual=0.00543532 Iteration: 20, Maximum residual=0.00360974 Iteration: 49, Maximum residual=0.00248598 Iteration: 98, Maximum residual=0.0015489 Iteration: 100, Maximum residual=0.00142966 Clean used 100 iterations to get to a max residual of 0.00142966 second 100: Initial maximum residual: 0.00142965 Iteration: 3, Maximum residual=0.0014203 Iteration: 13, Maximum residual=0.00126088 Iteration: 44, Maximum residual=0.00102948 Iteration: 100, Maximum residual=0.000829319 Clean used 100 iterations to get to a max residual of 0.000829319 RMS = 5.72e-4 third 100: Initial maximum residual: 0.000829345 Iteration: 2, Maximum residual=0.000773633 Iteration: 10, Maximum residual=0.0007373 Iteration: 37, Maximum residual=0.000643592 Iteration: 100, Maximum residual=0.000554785 Clean used 100 iterations to get to a max residual of 0.000554785 RMS = 5.35e-4 fourth 100: nitial maximum residual: 0.000554796 Iteration: 2, Maximum residual=0.000516688 Iteration: 10, Maximum residual=0.000482603 Iteration: 38, Maximum residual=0.000444863 Iteration: 100, Maximum residual=0.000395171 RMS = 5.16e-4 Stop here to see what the image looks like. Beam used in restoration: 1.83417 by 0.962932 (arcsec) at pa 16.6164 (deg) Looks like I should have used a smaller mask and I went too deep for this big mask (there are small blips around real src structure). RMS = 5e-4 Try again: Since initial dirty is so bad, use clean image to make a good mask: imgr.regionmask(mask="ggtau1.mask2" , region=myregion, value=1.0); Initial maximum residual: 0.00868421 Iteration: 2, Maximum residual=0.00703421 Iteration: 7, Maximum residual=0.00535603 Iteration: 18, Maximum residual=0.00368502 Iteration: 39, Maximum residual=0.0024278 Iteration: 77, Maximum residual=0.00154037 Iteration: 131, Maximum residual=0.000877265 Iteration: 229, Maximum residual=0.000530984 Iteration: 231, Maximum residual=0.000480833 Clean used 231 iterations to get to a max residual of 0.000480833 Beam used in restoration: 1.83417 by 0.962932 (arcsec) at pa 16.6164 (deg) Much better. ggtau1.robust.ps RMS = 0.65 mJy/bm, max = 8.96 mJy/bm grey = 1.2 mJy/bm to max contours at -2,2,4,6,8,10 sigma See what the image looks like with uniform weighting: imgr.weight(type="uniform"); imgr.clean(algorithm="clark" , niter=1000, gain=0.1, threshold=[value=0.5, unit="mJy" ], displayprogress=F, model="ggtau1.model3" , fixed=F, complist='', mask="ggtau1.mask2" , image="ggtau1.cm3" , residual="ggtau1.resid3" , interactive=F, npercycle=100); Initial maximum residual: 0.00868421 Iteration: 2, Maximum residual=0.00703421 Iteration: 7, Maximum residual=0.00535603 Iteration: 18, Maximum residual=0.00368502 Iteration: 39, Maximum residual=0.0024278 Iteration: 77, Maximum residual=0.00154037 Iteration: 131, Maximum residual=0.000877265 Iteration: 229, Maximum residual=0.000530984 Iteration: 231, Maximum residual=0.000480833 Clean used 231 iterations to get to a max residual of 0.000480833 Clean did not reach threshold beam used in restoration: 1.83417 by 0.962932 (arcsec) at pa 16.6164 (deg) RMS = 0.65 mJy/bm Max = 8.96 mJy/bm Humm, there is absolutely no difference between these 2 images. Try natural weighting: imgr.weight(type="natural"); imgr.clean(algorithm="clark" , niter=1000, gain=0.1, threshold=[value=0.5, unit="mJy" ], displayprogress=F, model="ggtau1.model4" , fixed=F, complist='', mask="ggtau1.mask2" , image="ggtau1.cm4" , residual="ggtau1.resid4" , interactive=F, npercycle=100); nitial maximum residual: 0.0116784 Iteration: 1, Maximum residual=0.0105106 Iteration: 3, Maximum residual=0.00937658 Iteration: 8, Maximum residual=0.00767836 Iteration: 16, Maximum residual=0.00585018 Iteration: 28, Maximum residual=0.00432898 Iteration: 46, Maximum residual=0.00297798 Iteration: 68, Maximum residual=0.00195901 Iteration: 94, Maximum residual=0.00122027 Iteration: 122, Maximum residual=0.000734084 Iteration: 144, Maximum residual=0.000494288 Clean used 144 iterations to get to a max residual of 0.000494288 beam 1.98562 by 1.05086 (arcsec) at pa 12.9671 OK, this is slightly different. The image doesn't look exactly like the one in the alma/iram document but I suspect this is because I've only got one image here, not a slew of data from different configurations.