! final_clean -- script for making deeply cleaned images
!                by Dan Homan,  dhoman@nrao.edu
!
! This started out as the automatic calibration and mapping
! script by Taylor and Shepard, but I have turned it into
! a multi-resolution cleaning script for producing "final"
! deeply cleaned images in difmap.  No self-calibration
! is performed by this procedure... its' only purpose is
! to perform a deep multi-resolution clean on calibrated data.
!
! The motivation for this procedure comes from tests
! run by George Moellenbrock which showed that 'super-resolving'
! the highest brightness temperature features early on in
! the cleaning process often lead to a much more accurate
! clean of the lower surface brightness features in later
! clean iterations at lower resolution.  'super-resolving'
! in the early rounds seemed to avoid unrecoverable errors
! in cleaning the most compact VLBI structures.
! (Moellenbrock, Ph.D. Thesis (1999))
!
! During the self-calibration process, I usually implement
! George's ideas with carefully controlled 'super-resolution'
! or by modelfitting of the VLBI cores.  For making final, deep
! cleaned images I created this automated (reproducable)
! procedure.
!
! To run it: load the data, set your mapsize, then just type
!
! 0> @final_clean i
!
! Where "i" can be replaced by another Stokes parameter.
! Any exisiting clean components are automatically deleted, 
! uvtapers are removed, and clean windows are deleted.  
! Of course, you can change this by editing the script.
!
! The lowest resolution used is simply 'natural
! weighting' of "uvweight 0,-2".  Convolving with
! the same beam as this lowest resolution is probably
! the safest thing to do.  As 'natural weighting' is
! is the last weighting used by the procedure, the
! model will be automatically convolved with 
! this beam unless a specific change of beam or
! resolution is made first.
!
! The procedure just cleans, so you must write the
! FITS image out of difmap yourself...
!
! 0> wmap filename
! 
! Of course the model components (clean components) can
! be saved separately if desired...
!
! 0> wmod filename.mod
!


integer clean_niter; 
float clean_gain; clean_gain = 0.03
float dynam;
float flux_peak;

! Define the inner loop as a macro.

float flux_cutoff
float dyn_range
float last_dyn_range

#+map_residual \
flux_peak = peak(flux);\
flux_cutoff = imstat(rms) * dynam;\
while(abs(flux_peak)>flux_cutoff);\
 clean clean_niter,clean_gain;\
 flux_cutoff = imstat(rms) * dynam;\
 flux_peak = peak(flux);\
end while

! The following macro is similar to above, but stops
! when there is no improvement in the dynamic range 
! of the residual map

#+deep_map_residual \
last_dyn_range = abs(peak(flux))/imstat(rms);\
clean clean_niter, clean_gain;\
dyn_range = abs(peak(flux))/imstat(rms);\
while(last_dyn_range > dyn_range);\
 clean clean_niter,clean_gain;\
 last_dyn_range = dyn_range;\
 dyn_range = abs(peak(flux))/imstat(rms);\
end while


! select the stokes to clean
select %1

! clear previous model
clrmod true,true,true

! delete any windows
delwin

! remove any tapering
uvtaper 0

print "*********** FIRST TRY SUPER-UNIFORM WEIGHTING **********"
print "**** -- only if dynamic range is higher than 10 -- *****"

dynam = 10
clean_niter = 10
uvw 20,-1
map_residual
uvw 10,-1
map_residual
clean_niter = 50

print "*********** REGULAR UNIFORM WEIGHTING NEXT ***************"
uvw 2,-1
dynam = 6
map_residual
print "********** DEEP CLEANING AT NATURAL WEIGHTING **************"
uvw 0,-2
! set dynamic range and iteration cycle so the clean will go deep
clean_niter = 500  
deep_map_residual

print "********** FINAL CLEAN IS FINISHED **************"