!
! Difmap script automap
!
! Usage: @automap.com <classver> <source> <dointer>
!
!        <classver> CLASS1 CLASS2a CLASS2b
!        <source>   eg. 16156554
!        <dointer>  true/false
!
!-----------------------------------------------------------------------
! 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).
!  dointer        logical    Display intermediate maps
!  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
!
! Version:
!
! V.16Jun94 STM from CDF
! V.20Jun94 STM, second self-cal loop inside field loop
! V.05Jul94 STM, units arcsec
! V.07Jul94 STM, wtscal 4901.5 for VLA 25MHz 10s 1-mult vis
! V.08Jul94 STM, maprms-->imstat; uvstat for modelfit iterations
! V.14Jul94 STM, version for Jodrell Bank installation
! V.15Jul94 STM,CDF peakmax, 512 small maps, break out of search correctly
! V.16Jul94 STM,CDF add absolute cutoff on search of 10mJy, no FITS output
! V.09Aug94 CDF, modify to run with Caltech directory structure
!
! V.30May97 STM, modify to run with UPenn cosmo directories 
!                output files to current directory
! V.12Jun97 STM, modify to run with UPenn beavis directories 
! V.15Jun97 STM, output gif file
!
!-----------------------------------------------------------------------

! Assign command-file arguments to variables.

!string day;          fprint 0,"day:  ";       day = read_string
!string source;       fprint 0,"source:  ";    source = read_string
!integer large_size; fprint 0,"large_size: "; large_size = read_int
!float large_cell;   fprint 0,"large_cell: "; large_cell = read_float
!integer field_size; fprint 0,"field_size: "; field_size = read_int
!float field_cell;   fprint 0,"field_cell: "; field_cell = read_float
!integer niter;      fprint 0,"niter: ";      niter = read_int
!float gain;         fprint 0,"gain: ";       gain  = read_float
!float cutoff1;       fprint 0,"cutoff: ";     cutoff1= read_float
!float soltime;      fprint 0,"soltime: ";    soltime=read_float

string classver;     classver = "%%1"
string source;       source = "%%2"
logical dointer;     dointer = %3             

! UPenn cosmo directories
string datadir;      datadir = "/data/myers/CLASS/Data/"
string outdir;       outdir = "./"

! UV weight scale factor
float classwt;      classwt = 4901.5

! Set up run parameters
string uvfile
uvfile =  datadir // classver // "/" // source // ".UVF"
integer large_size; large_size = 1024
float large_cell;   large_cell = 0.250
integer field_size; field_size = 512
float field_cell;   field_cell = 0.050
integer niter;      niter = 25
float gain;         gain = 0.05
float cut_abs;      cut_abs = 0.010
float cut_cal;      cut_cal = 5.0
float cut_rms1;     cut_rms1 = 8.0
float cut_dyn1;     cut_dyn1 = 15.0
float cut_rms2;     cut_rms2 = 6.0
float cut_dyn2;     cut_dyn2 = 20.0
float plot_sig;     plot_sig = 3.0
float win_size;     win_size = 0.75

logical dofields

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;\
  uvrange 0,51.6/large_cell

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

! 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
if (dointer) mapl

!
! 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)
  if (dointer) mapl

!
! 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
  if (dointer) mapl

  while(peak(flux,max)>=(pkflux1/cut_cal) & peak(flux,max)>=cutoff1)
    peakwin win_size
    clean niter,gain
    if (dointer) mapl
  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
  if (dointer) mapl

  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
!
  cutoff1 = cut_rms1 * imstat(rms)
  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)

  if (dointer) mapl

  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)
	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
        if (dointer) mapl
        cutoff2 = cut_rms2 * imstat(rms)
        print "Residual Field Peak =",peak(flux,max)," RMS =",imstat(rms)
        print "New cutoffs Dyn-Range =",limit2," RMS = ",cutoff2
      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

      if (dointer) mapl

    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
!
    psfile =  "\""// outdir // source // "." // strnum(i) // ".ps\"/vps"
    if (!dointer) dev `psfile`
    mappl cln
    if (!dointer) 
        giffile =  "\""// outdir // source // "." // strnum(i) // ".gif\"/gif"
	dev `giffile`
    	mappl cln
    	dev "/null"
    end if
! 
! Save map
!
!    mapfile =  outdir // source // "." // strnum(i) // ".fits"
!    wmap `mapfile`

  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
  if (!dointer) 
    giffile =  "\""// outdir // source // ".0.gif\"/gif"
    dev `giffile`
    mappl 
    docont = true
    dev "/null"
  end if

end if

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