! DJR 09.23.98
! Difmap script automap (version for NRAO)
!
! Usage: @automap.com <datadir> <source> <weight>
!
!        <datadir>  ex: /data/beavis_1/myers/CLASS/drusin/automap
!        <source>   0000000
!        <weight>   ex. 5527 or 11053
!
!-----------------------------------------------------------------------
! Hands off mapping procedure for reasonably well calibrated and edited
! data. This version has an extra outer loop to search for sources at
! large distances from the phase center in order to implement a crude
! form of multi-field mapping.
!
! Input:
!  uvfile  literal-string    The name of the UV FITS file.
!  large_size         int    A map dimension used to cover a wide search area.
!  large_cell       float    The map cell-size to use with 'large_size'.
!  field_size         int    The map dimension for a single field.
!  field_cell       float    The cell size to use with 'field_size'.
!  niter              int    The number of CLEAN-iterations per cycle.
!  gain             float    The CLEAN loop gain.
!  cutoff           float    The residual cut-off flux for.
!  soltime          float    Solution time for phase selfcal (minutes).
!  classwt          float    Weight scaling factor (given BW and band)
! Output Files:
!  source.win      Window file for rwin
!  source.mod      Clean model file
!  source.gmod     Modelfit model file
!  source.#.ps     Postscript maps each field
!  source.#.gif    Postscript maps each field (0=no sources found)
!
! Note: for gif file, be sure to set PGPLOT_GIF_WIDTH and _HEIGHT 
!       environment variables


string datadir;       datadir = "%%1"
string source;       source = "%%2"
float classwt;       classwt = %3             

string datadir;      datadir = datadir // "/"
string outdir;       outdir = "./" 


! float classwt;      classwt = 5527
! UV weight scale factor for I 50Mhz
! float classwt;      classwt = 11053

! Set up run parameters
string uvfile
uvfile =  datadir // "/" // source // ".UVF"

! Large field +/- 64"
integer large_size; large_size = 1024
float large_cell;   large_cell = 0.25
! Small field +/- 6.4"
integer field_size; field_size = 512
float field_cell;   field_cell = 0.05
! Default clean 25 inter at 5% gain
integer niter;      niter = 25
float gain;         gain = 0.05
! Clean window 75% beam
float win_size;     win_size = 0.75
! Absolute cutoff 10mJy
float cut_abs;      cut_abs = 0.010
! S/N and Dyn range cutoffs for big field
float cut_rms1;     cut_rms1 = 8.0
float cut_dyn1;     cut_dyn1 = 15.0
! S/N and Dyn range cutoffs for small field
float cut_rms2;     cut_rms2 = 6.0
float cut_dyn2;     cut_dyn2 = 20.0
! cutoff for cleaning to fraction of peak
float cut_cal;      cut_cal = 5.0
! Plot contours etc.
float plot_sig;     plot_sig = 3.0

logical dofields

! Modelfitting overlap and convergence parameters
float mod_dist_min; mod_dist_min = 0.2
float mod_chi_frac; mod_chi_frac = 0.001
float mod_chi_abs;  mod_chi_abs = 0.75

integer mod_iter_inc;   mod_iter_inc = 5
integer mod_iter_max;   mod_iter_max = 100
integer mod_iter

float chisq_prev
float chisq_next

float xfield
float yfield
float pkflux0
float pkflux1
float pkflux2
float cutoff1
float cutoff2
float limit1
float limit2
string psfile
string giffile
string mapfile
string modfile
string winfile

! Define some macros for maps

!
! Unshifted Big grid with super-uniform weighting 
!
#+search_map \
  unshift;\
  mapsize large_size, large_cell;\
  uvw 4,-2

!
! Small grid at current shift with uniform weighting
!
#+field_map \
  mapsize field_size, field_cell;\
  uvw 2,-2

! Set up variables to record the positions of the field centers
float xbig
float ybig
float xcent(25)
float ycent(25)
integer count; count = 0
integer field; field = 0

!
! Set up variables for final modelfitting
!
float peak_x(25)
float peak_y(25)
float peak_f(25)
integer n_peak; n_peak = 0
float p_x
float p_y
float p_f
logical found
float tol; tol = mod_dist_min^2
float d_x
float d_y

! Set units to arcsec (and Klambda)
mapunits arcsec

! Read the UV data.

print ""
print "*************************************************"
print "*************************************************"
print "!!! Starting Source ",source
print "!!! From UVfile ",uvfile
print "*************************************************"
print "*************************************************"
print ""

observe `uvfile`
select I

! Scale the weights
wtscal classwt

!
! This first loop locates the brightest peak field and uses it for
! an initial self-cal before the real searching
!
search_map

!
! Set rms cutoff to CUT_RMS1 sigma of big map, and peak
!
cutoff1 = cut_rms1 * imstat(rms)
pkflux1 = peak(flux,max)
print "!!! Initial Peak =",pkflux1, "RMS =",imstat(rms), \
   "S/N =", pkflux1/imstat(rms)
print "!! Expected RMS =", imstat(noise), " Search cutoff =",cutoff1

!
! Check to see if there is a valid peak in the map
!
if (pkflux1 < cutoff1 & pkflux1 < cut_abs) 

  print "!!! Final Search Peak= ",peak(flux,max), " RMS =",imstat(rms)
  print "!!! No Valid peaks found above CUTOFF = ", cutoff1
!  print "!!! "
!  stop     ! this command didnt work
  dofields = false

else
  
  if ( pkflux1 <  cutoff1 ) print "!!! Peak above absolute cut = ", cut_abs

!
! Shift to the position of the brightest point in the wide map.
! Make smaller map then shift to that peak
!
  shift -peak(x,max), -peak(y,max)
  field_map
  shift -peak(x,max), -peak(y,max)

!
! Set rms cutoff to CUT_RMS2 sigma of small map, and peak
!
  cutoff1 = cut_rms2 * imstat(rms)
  pkflux1 = peak(flux,max)

  print "!! Initial Field Peak =",pkflux1, "RMS =",imstat(rms), \
     "S/N =", pkflux1/imstat(rms)
  print "!! Expected Field RMS =", imstat(noise), " Search cutoff =",cutoff1

!
! Start mapping the field
!
  peakwin win_size
  clean niter,gain
  clean niter,gain

  while(peak(flux,max)>=(pkflux1/cut_cal) & peak(flux,max)>=cutoff1)
    peakwin win_size
    clean niter,gain
  end while

  print "!! Residual Field Peak =",peak(flux,max), "RMS =",imstat(rms)

!
! 30-second initial phase-only selfcal
!
  selfcal false,false,0.5

  print "!! Residual Selfcal Peak =",peak(flux,max), "RMS =",imstat(rms)

  delwin
  clrmod true

  print "***********DONE WITH INITIAL PHASECAL************"
  
  dofields = true

end if

!
! THIS IS THE OUTER SEARCH LOOP
!

while(dofields)
  count = count + 1
!
! Allow a maximum of 25 iterations
!
  if(count > 25)
    break
  end if
!
! Set up large search map
!
  search_map

  pkflux1 = peak(flux,max)
! Store absolute map peak if first iteration
  if ( count == 1 ) 
    pkflux0 = pkflux1
    print "!!! SelfCal Peak =",pkflux1, "RMS =",imstat(rms), \
         "S/N =", pkflux1/imstat(rms)
  end if

!
! rms and dyn-range cutoffs for search maps
!
! CUT_RMS1 x map RMS cutoff
  cutoff1 = cut_rms1 * imstat(rms)
! Peak / CUT_DYN1 dyn range limit
  limit1  = pkflux0/cut_dyn1
  print " "
  print "!!",count,": Search Peak =", pkflux1, "RMS =",imstat(rms), \
      "S/N =", pkflux1/imstat(rms)
  print "!!",count,": Search Cutoffs Dyn-Range =",limit1,"RMS =", \
      cutoff1

  print " "

!
! If the peak alsolute flux in the wide map is lower than the cutoff(s), stop.
!
  if( pkflux1 < limit1 | ( pkflux1 < cutoff1 & pkflux1 < cut_abs ) )
    break 
  end if

!
! Shift to the position of the brightest point in the wide map.
! Set up small field map there
!
  xbig = peak(x,max)
  ybig = peak(y,max)
  shift -peak(x,max), -peak(y,max)
  field_map
  xfield = xbig + peak(x,max)
  yfield = ybig + peak(y,max)
  shift -peak(x,max), -peak(y,max)

!
! Record the peak flux of this field.  Note that pkflux0 is used to
!  determine a cutoff level.
!
  pkflux2 = peak(flux,max)

  limit2 = pkflux0/cut_dyn2
  cutoff2 = cut_rms2 * imstat(rms)
  print "!!",count,": Field Peak =", pkflux2, "RMS =",imstat(rms), \
     "S/N =", pkflux2/imstat(rms)
  print "!!",count,": Field cutoffs Dyn-Range =",limit2,"RMS =",cutoff2
  print "!!",count,": Expected RMS =", imstat(noise)

!
! Start mapping the field, using FIELD weighting.
!
  if( peak(flux,max) >= limit2 & peak(flux,max) >= cutoff2 )

! Increment counter
    field = field + 1

    print ""
    print "*************************************************"
    print "!!! Field number:",field," at ",xfield,yfield," Peak =", \
         peak(flux,max)
    print "*************************************************"
    print ""

!
! Inner FIELD loop 
!
    while( peak(flux,max) >= limit2 & peak(flux,max) >= cutoff2 )

!
! Inner-inner SELF-CAL loop
!
      while( peak(flux,max) >= limit2 & peak(flux,max) >= cutoff2 )
        peakwin win_size
! Store peaks for later modelfit
	p_x = xfield + peak(x,max)
	p_y = yfield + peak(y,max)
        p_f = peak(flux,max)

	print "***n_peak beginning= ", n_peak
	print "***p_ = ", p_x, p_y, p_f

	if( n_peak > 0)
	  found = false
	  do i = 1,n_peak
	    d_x = abs(p_x-peak_x(i))
	    if( d_x < mod_dist_min )
	      d_y = abs(p_y-peak_y(i))
	      if( d_y < mod_dist_min )
		if( (d_x^2 + d_y^2) < tol ) found = true
	      end if
            end if
	    if(found) break
	  end do
	  if(!found)
	    n_peak = n_peak + 1
	    peak_x(n_peak) = p_x
	    peak_y(n_peak) = p_y
	    peak_f(n_peak) = p_f
	  end if
	else
	  n_peak = 1
	  peak_x(n_peak) = p_x
	  peak_y(n_peak) = p_y
	  peak_f(n_peak) = p_f
	end if
! Clean
        clean niter,gain
        cutoff2 = cut_rms2 * imstat(rms)
        print "Residual Field Peak =",peak(flux,max)," RMS =",imstat(rms)
        print "New cutoffs Dyn-Range =",limit2," RMS = ",cutoff2

	print "***n_peak end= ", n_peak

      end while
!
! Phase Selfcal 10s
!
      selfcal false, false, 10./60.

      cutoff2 = cut_rms2 * imstat(rms)
      print "After Self-Cal Peak =",peak(flux,max)," RMS =",imstat(rms)
      print "New cutoffs Dyn-Range =",limit2," RMS =",cutoff2

    end while

! Remember some things
    xcent(field) = xfield
    ycent(field) = yfield

    print ""
    print "*************************************************"
    print "!! Finshed Field number:",field," Final RMS =",imstat(rms)
    print "*************************************************"
    print ""

  else
!
! Nothing found in this field - break search loop
!
    print ""
    print "!! Nothing found in field - breaking search loop"
    print ""
    break

  end if

end while


print ""
print "*************************************************"
print " Completed Field Searches, found ",field
print "*************************************************"
print ""

print "!!! Final Search Map Peak =",peak(flux,max), " RMS =",imstat(rms)

! Plot the contour plots and save field maps

if (field > 0)

  do i=1,field
    unshift 
    shift -xcent(i),-ycent(i)
    field_map
!
! Make lowest contour PLOT_SIG*sigma in absolute contours
!
    cmul = imstat(rms)
    logl plot_sig, pkflux0/cmul
!
! Make postscript and gif files
!

! Uncomment in order to make ps file

    psfile =  "\""// outdir // source // "." // strnum(i) // ".ps\"/vps"
    dev `psfile`
    mappl cln

    giffile =  "\""// outdir // source // "." // strnum(i) // ".gif\"/gif"
    dev `giffile`
    mappl cln
    dev "/null"

  end do

!
! Save model and windows
!
  modfile = outdir // source // ".mod"
  wmod `modfile`

  winfile = outdir // source // ".win"
  wwin `winfile`

  print ""
  print "*************************************************"
  print " Modelfit on ",n_peak," peaks"
  print "*************************************************"
  print ""

! Set up initial model at location of peaks
  clrmod true
  unshift
  do i = 1,n_peak
    addcmp peak_f(i),true,peak_x(i),peak_y(i),true,0.02,true,1,true,0,true
  end do
! Get initial chisq
  chisq_next = uvstat(chisq)
  chisq_prev = 2.0*chisq_next
  mod_iter = 0

  print "!!! Modelfit ",n_peak," components, initial chi-squared = ", \
     chisq_next

! Modelfit loop
  while ( chisq_next >= mod_chi_abs & (chisq_prev-chisq_next) >= \
       mod_chi_frac*chisq_prev & mod_iter < mod_iter_max )
    modelfit mod_iter_inc
    mod_iter = mod_iter + mod_iter_inc
    chisq_prev = chisq_next
    chisq_next = uvstat(chisq)
  end while

  print "!!! Modelfit ",mod_iter," iterations, final chi-squared = ", \
     chisq_next

  modfile = outdir // source // ".gmod"
  wmod `modfile`

else

  print "!!! NO SOURCES FOUND FOR SOURCE ",source
!
! Make gif file of large map for blank field
!
  unshift
  docont = false

  giffile =  "\""// outdir // source // ".0.gif\"/gif"
  dev `giffile`
  mappl 
  docont = true
  dev "/null"

end if

! Done
print ""
print "*************************************************"
print "!!! Completed Source ",source
print "*************************************************"
print "!!! "
print ""