$$$$   Procedures to reduce MERLIN data
$
$ This code is a mess.  I ahve butchered Phil's nice conditional
$ procedures with special cases and lost track of unused variables.
$ Sorry!  Will tidy and comment one day.... but it works!
$ Anita Aug 2003 updated Aug 2005
$ See MERLIN.HLP for user notes and limitations.
$------------------------------------------------------------------
vers 'runfil'
inn vers
tput addif
restore 0
clrmsg
dowait 1
$$$ * 


proc declare
$ No.Ants
  scalar uvdisc, seq, antref, uv2disc, aflux, ntel, gbnd, doflux
$           *,(UV)INSE,    *, (notused),  *,    *, (?see isline)
  scalar dobp, dophsc, doedit, isn, jsn, acell, dotgsc, scb, sce
$        *    , *      , *    ,SNVER,SNVER,CELLS,*(tarSC),*B/ECH
  scalar ntabs, k1, k2, disc, ii, jj, k, m, defuv, srchcl, gflux
$                      ,                ,UVRAN(1),largeCELLS
  scalar tseq, texist, kk, nn, line, fbchan, fechan,douvr, seetv
$ CUTO(undef?),                      ,(?isline             ),*     
  scalar doplot, dowobbl, dosearch, dorewt, mdisc, mseq, pkval
$         *    ,seeRUNIM,(*,SHIFT) ,*WTMO
  scalar dosmod, notuch, doplzn, ic, doload, trgvel, calok
$ scalar sunspot,
  array  wexist(4),mblc(3),mtrc(3),newwt(8),trgfreq(2)
  string*1 test
  string*2 tabtyp, tabt2, table, extype, ttype, oband
  string*3 mode
  string*4 tintp, stype, smmode, tsampty
  string*6 tclass, tclass1, tmpl, tver, lver, orver, class
  string*6 dir, multtb, uvtb, oclass
  string*12 sauce, tname, tname1, wname(4), target, acal, phcal
  string*12 pacal, mname
  string*60 plname
  string*40 pounds, text, excl
$
$  variables reserved for use in small procedures
  scalar v1, v2, itab, jtab
$
  scalar boxsi, cellf, decsi, loop, mbchan, mcells, mechan, mdotv;
  scalar mflux, mflx0, mgainu, mnboxes, mniter, mnit0, mrefant,
  scalar msolint, msol0, frate, rmsz, mbmaj, mbmin, mbpa, mdocal;
  scalar pniter, mdoband, dochans
  scalar mwtuv, nscal, rasi, rms, x1, x2, xy, y1, y2;
  array mbaddisk(10), mclbox(4,50), mimsize(2), timsize(2);
  array mzerosp(5), rmsb(4), muvt(2), mshift(2), muvrange(2);
  string*6 mapclas, mstokes; string*2 muvwtfn;
  string*50 outlin
ret; fin

proc checkin
  pounds='########################################'
  excl  ='!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  multtb='multtb'; uvtb='uvtb'; antref=refant; refant=0;
  line=-1; mseq=1; mdisc=indisk; defuv=0; mdoband=dobp
  scb=bchan; bchan=1; sce=echan; echan=0; timsize=imsize;
  mdotv=dotv;dotv=0
  if (imsize(1)=0) then; timsize(1)=512;end
  if (imsize(2)=0) then; timsize(2)=512;end
  imsize=512; msol0 = 5; mnit0 = 5000; muvrange = 0; mwtuv = 0.05;
  mshift=0; mbchan = 1; mechan = 0; mrefant = antref; mstokes = 'I'; 
  muvwtfn = 'NA'; mgainu = 0; mbaddisk = 0; mdocal = -1;
  dowobbl = 1; gbnd=1; dosmod = -1; trgvel=sysvel; sysvel=0;
  trgfreq=restfreq; restfreq=0;mbmaj=bmaj;mbmin=bmin;mbpa=bpa
  subarray 0;bmaj 0;bmin 0;bpa 0
$ *** till line fixed
  if (DOCHANS=1) then; DOPLOT=0;end
$
 mcells=acell; mapclas='icl001';
 if (stokes='LL') then mapclas='LCL001'; end;
 if (stokes='RR') then mapclas='RCL001'; end;
 if (stokes='Q')  then mapclas='QCL001'; end;
 if (stokes='U')  then mapclas='UCL001'; end;
$
$  Declare versions
$
  orver = version
  tver = 'TST';
  lver = 'APES';
ret; fin
$------------------------------------------------------------------

$------------------------------------------------------------------
$  Frequently used procedures
proc prpound; typ pounds; ret; fin
proc prexcl; typ excl; ret; fin

proc mess(text)
  prexcl; typ text, ; prexcl; 
ret; fin

$SETUVD: Set Name adverbs for input UV-Data
proc setuvd
 inname=mname;indisk=uvdisc;inclass=class;inseq=seq;intype='uv'
ret;fin
$SETMAP: Set Name adverbs for input UV-Data
proc setuvd
 inname=mname;indisk=uvdisc;inclass=class;inseq=seq;intype='ma'
ret;fin

$Get band from uvdata
proc GETBAND
SETUVD
task 'GETH';keystr '';keyw 'CRVAL3';keyval 0;GETH
if ((keyval(1)> 0.3e9)  & (keyval(1) <= 0.7e9)) then; oband='P';end
if ((keyval(1)> 0.7e9)  & (keyval(1) <= 3e9)) then; oband='L';end
$if ((keyval(1)> 3e9) & (keyval(1) <= 9e9)) then; oband='C';end
if ((keyval(1)> 3e9) & (keyval(1) <= 6.001e9)) then; oband='C5';end
if ((keyval(1)> 6e9) & (keyval(1) <= 9e9)) then; oband='C6';end
if ((keyval(1)> 19e9) &  (keyval(1) <= 26e9)) then; oband='K';end
fin

$To find total No antennas in AN
proc GETANTNO
task 'GETTHEAD';inext 'an';inver 1
keyword 'NUM ROW';GETT;NTEL = KEYVAL(1)
fin
$
$To set default refant to Mk
proc DEFREFANT
if antref = 0 then
 SETUVD
 GETANTNO
 task 'TABGET'
  for i=1:NTEL;inext 'an';inver 1;pixxy i 1 1;tabget
   if keystr = 'MK2' then; antref= i;mrefant=antref
  end
 end
end
if antref = 0 then
  for i=1:NTEL;inext 'an';inver 1
  pixxy i 1 1;tabget
   if (keystr='DARNHALL' ! keystr='TABLEY' ! keystr='KNOCKIN') then; 
   antref= i;mrefant=antref
  end
 end
end
fin
$$
$To set antwts
proc SETWTS
task 'TABGET'
mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc
SETUVD
GETANTNO
inext 'an';inver 1
if (oband='P') then
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DEFFORD' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'CAMBRIDG' 
then newwt(i)=1.64;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'KNOCKIN' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'WARDLE' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DARNHALL' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' 
then newwt(i)=9.25;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'TABLEY' 
then newwt(i)=1.0;end;end
end
if (oband='L') then
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DEFFORD' 
then newwt(i)=0.61;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'CAMBRIDG' 
then newwt(i)=1.74;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'KNOCKIN' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'WARDLE' 
then newwt(i)=0.73;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DARNHALL' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' 
then newwt(i)=30.25;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'TABLEY' 
then newwt(i)=0.77;end;end
end
if (oband='C5') then
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DEFFORD' 
then newwt(i)=0.1;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'CAMBRIDG' 
then newwt(i)=4.8;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'KNOCKIN' 
then newwt(i)=1.7;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'WARDLE' 
then newwt(i)=0.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DARNHALL' 
then newwt(i)=1.7;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' 
then newwt(i)=9.25;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'TABLEY' 
then newwt(i)=1.5;end;end
end
if (oband='C6') then
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DEFFORD' 
then newwt(i)=0.1;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'CAMBRIDG' 
then newwt(i)=4.8;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'KNOCKIN' 
then newwt(i)=1.7;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'WARDLE' 
then newwt(i)=0.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DARNHALL' 
then newwt(i)=1.7;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' 
then newwt(i)=9.25;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'TABLEY' 
then newwt(i)=1.5;end;end
end
if (oband='K') then
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DEFFORD' 
then newwt(i)=0.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'CAMBRIDG' 
then newwt(i)=1.64;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'KNOCKIN' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'WARDLE' 
then newwt(i)=0.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'DARNHALL' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' 
then newwt(i)=1.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' 
then newwt(i)=0.0;end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'TABLEY' 
then newwt(i)=1.0;end;end
end
fin
$
proc giveup
doconfrm 0
read doconfrm
if (doconfrm > 0) then;source(1) 'stop';inn source(1);imh;
inn ''
end
fin
$
proc tbfind (tabt2, ntabs)
   ntabs = 0
   for v1 = 1 to 32
      keyword = 'extype'!!char(v1); geth; extype = substr(keystrng,1,2)
      if (extype = tabt2) then
         keyword = 'extver'!!char(v1); geth; ntabs = keyval(1)
$         nplots=ntabs
         v1 = 32
         end
      keyword=' '; keyval=0; keystr=' '; 
      end 
ret; fin
$
proc rprint
   if (doplot = 1) then;docrt -1;outpr '';end
   if (doplot = 2) then;docrt 1;end
   if (doplot = 3) then;docrt -1;
   outpr  DIR!!':'!!plname!!'.txt';end
fin

proc rplot(itab, jtab)
   if (doplot > 0) then
      dparm 0; outfile ''; lpen 3; aspmm 0; functype ' '; copies 1;
      userid 0;inver 0 
      if (doplot=1) then
         for plver = itab to jtab; go lwpla; end
      end
      if (doplot=2) then
         for plver = itab to jtab; go tkpl; end
      end
      if (doplot=3) then
      for plver = itab to jtab; 
      outfil DIR!!':'!!plname!!'_'!!char(plver)!!'.ps'
      go lwpl;end;end
   end
ret; fin

proc doextd (tabtyp, itab, jtab)
   inext = tabtyp;
   for invers = itab to jtab; extd; end;
ret; fin

$Fix multi files for source names > 8 char
proc SNAMFIX
if (length(inn) > 8) then
task 'tabed'
outn inn;outcl incl;outd ind;outse inse;inver 1;outver 0
bcount 1;ecount 1;opty 'repl';aparm 2 0;keystr inn
inext 'su';keywo '';keyval 0;timer 0;go tabed
end
fin

proc IFMULT
docal -1
if (substr(incl,1,3)='mul') then;docal 1;gainuse 0;end
fin

proc exist (disc, tname, tclass, tseq, texist)
$    disc,tname,tclass,tseq = disc, mname, class, seq. no.
   userid=0; inname=tname; inclass=tclass; inseq=tseq;
   intype=' '; indisk=disc; chkname; 
   if (tname=' ') then; error=1; end
   if (error=0) then; texist=true; else; texist=false; end;
   ret; fin

proc dozap (disc, tname, tclass, tseq, ttype)
$ delete a file if it exists
$ disc,tname,tclass,tseq = disc, mname, class, seq. no.
$ type = 'UV' or 'MA'
   exist (disc, tname, tclass, tseq, texist);
   if (texist) then
      userid=0; inname=tname; inclass=tclass; inseq=tseq;
      indisk=disc; intype=ttype; 
      type 'destroying: ',innam,inclass,inseq
      zap; recat;
      end;
   ret; fin

proc allzap (disc, tname, tclass, k1, k2, ttype)
$ delete file(s) from seq num k1 -> k2 if they exist
$ disc,tname,tclass = disc, mname, class, 
$ type = 'UV' or 'MA'
   for i = k1 to k2
      dozap (disc, tname, tclass, i, ttype)
   end
ret; fin

$ proc isline (line, fbchan, fechan, gbnd)
$ $ is this a spectral line file? If so, what is recommended
$ $ start channel for channel selection to avoid bandpass
$    v2 = 0; fbchan=1; gbnd=-1
$    for v1 = 1 to 3
$       keyword='CTYPE'!!char(v1); geth
$      smmode = substr(keystrng,1,4)
$      if (smmode = 'FREQ') then
$         v2=v1
$         keyword='NAXIS'!!char(v2); geth; v2=keyval(1)+keyval(2)
$         v1 = 3
$      end
$   end
$   if (v2>1) then; line=1; end
$   if (line=1) then 
$      fbchan=floor(v2/16)+1; fechan=v2; gbnd=1; 
$      if (mdoband >= 0) gbnd = -1
$   end
$   ret; fin
$
proc isline (line, fbchan, fechan, gbnd) 
$ Is this multi-channel data?   
   v2 = 0; fbchan=1; gbnd=-1
   for v1 = 1 to 3
      keyword='CTYPE'!!char(v1); geth
      smmode = substr(keystrng,1,4)
      if (smmode = 'FREQ') then
         v2=v1
         keyword='NAXIS'!!char(v2); geth; v2=keyval(1)+keyval(2)
         v1 = 3
      end
   if (v2>1) then; line=1; end
   if (line=1) then 
      gbnd=1; 
      mbchan=ceil(v2/32)+1;mechan=floor(v2*31/32)
   end
$ Is DOBAND expressly forbidden?
  if (MDOBAND < 0) then; gbnd = -1;
  end
$ If no spectral line specified:
   if (scb=1) then;scb=mbchan;end
   if (sce=0) then;sce=mechan;end
   end
   ret; fin

$ Look for 3C286 and set flux density
proc SET3C
$ *** check for K and P
innam=TARGET; incl='MULTTB'; ind=mdisc; ins=mseq; 
task 'GETTHEAD';inext 'SU';inver 1
keyword 'NUM ROW';GETT;ii = keyval(1)
i=1;m=0
inext 'SU';pixxy 0 2 1;keystr ''
for i=1:ii
pixxy(1)=i;go TABGET
if KEYSTR = '3C286' then
m=i;i=ii
source '3C286''';zerosp 0; optype 'CALC';aparm 0;go SETJY
$ Set flux of 3C286 to be 96% of VLA value at C-band 
if (oband='C5') then
inext 'SU';keystr '';pixxy m 5 1;go TABGET
tget SETJY
opty '';zerosp KEYVAL(1), 0; zerosp(1)=0.96*zerosp(1)
go SETJY
end
$
if (oband='C6') then
inext 'SU';keystr '';pixxy m 5 1;go TABGET
tget SETJY
opty '';zerosp KEYVAL(1), 0; zerosp(1)=0.92*zerosp(1)
go SETJY

end
end
end
fin


proc caldef
$ set default parameters of calibration parameters
   timer 0; selband -1; selfreq -1; freqid -1; suba 0; qual -1;
   calcode ''; docal -1; gainuse 0; dopol -1; blver -1; flagver 0; 
   doband -1; bpver 0; smooth 0; stokes ''; bchan 1; echan 0; bif 1;
   eif 0; antennas 0; polplot ''; 
   ret; fin

proc calprmdef
$ set default parameters for calib
   dofit 0; uvr 0; wtuv 0; clr2n; inver 0; ncomp 0; nmaps 0;
   cmethod ''; cmodel ''; smodel 0; clron; baddisk 0; soltype ''; 
   solmode 'p'; solcon 0; cparm 0; cparm(2)=1; snver 0; antwt 0; 
   gainerr 0; aparm 3 0; refant antref; solint 0.5; solmode 'p'; 
   aparm(7)=1; ret; fin

proc imagdef
$ set default parameters for IMAGR
   channel 0; clron; clr2n; nfield 1; fldsiz 0; 
   npoints 1;nchav 1;chinc 1
   if (DOCHANS=2) then;
      isline (line, fbchan, fechan, gbnd)
      npoints=(mechan-mbchan+1); nchav=(mechan-mbchan+1);
      chinc=(mecha-mbchan+1);
   end
   rashift 0; decshift 0; uvtaper 0; uvr 0; guard 0; rotate 0; zerosp 0
   uvwtfn ''; uvsize 0; robust 0; uvbox 0; uvbxfn 1; xtype 5; ytype 5; 
   xparm 0; yparm 0; niter 0; bcomp 0; nbox 0; clbox 0; boxfile ''; 
   gain 0.1; flux 0; minpatch 255; bmaj 0; bmin 0; bpa 0; phat 0; 
   factor 0; cmethod ''; cparm 0; maxpixel 0; dotv -1; outver 0;
   baddisk 0; antuse 0; ret; finish

proc runspl (sauce,oclass,i,j,k,jj,mbchan,mechan)
$ run SPLIT
$    sauce = SOURCE to split; oclass = o/p class mname
$    i = docalib; j = mdoband; k = aparm(1); jj=dopol
$   version = tver; 
   task 'SPLIT'; setuvd; caldef; 
   outclas=oclas; outd=ind; outs=ins; docal=i; gainuse 0; doband=j; 
   bpver 0; aparm 0; aparm(1)=k; bchan=mbchan; echan=mechan;
   douvc -1; ichansel mbchan, mechan, 1; source ''; aparm(6) 1
   source(1)=sauce; dopol=jj; 
   go
   ret; fin

proc runcl (i,j,k,tintp,tsampty,stype)
$    set input name parameters before calling runcl
$    i = SN to use; j = CL to update, new CL = k
$    tintp = interpolation method
$    tsampty = smoothing function
$    stype  = smotype
$                                                
$   version = tver; 
   task 'clcal'; setuvd
   opcode 'CALI';source=' '; soucode=' '; calsour=' '; qual -1;inver 0
   calcode ' '; timer 0; suba 0; anten 0; selband -1;doblank=-1 
   selfreq -1; freqid -1; opcode='cali'; interpol=tintp; 
   intparm 0; bparm 0;samptype=tsampty; smotype ' '; refant=antref; 
   snver=i; gainver=j; gainuse=k; baddisk 0;cutoff 0;icut 0 
   intparm=0.1667; bparm=intparm; smotype=stype;dobtween 0
   if ((smotype='AMPL') & (interpol<>'self')) then; 
   samptype = 'mwf';intparm(1)=0.5;bparm=intparm; end
   if (smotype=' ') then; intparm=0;bparm 0; end
   go; ret; fin

proc coptb (tname1,tclass1,j,jj,table, k1, k2)
$    set input name parameters before calling coptb
$    tname1, tclass1, j, jj = o/p file description
$    table = type of table; 
$    k1, k2 = inver, outver
$
$   version = tver; 
   task 'tacop'; setuvd
   outdisk=j; outnam=tname1; outcl=tclass1; outseq=jj;
   keyword ''; keyval 0; keystr ''; inext=table; inver=k1; 
   outver=k2; ncount=1; go; 
ret; fin
proc load (sauce)
$ sauce = source name
$   dozap (mdisc, sauce, uvtb, mseq, 'UV')
$   dozap (mdisc, sauce, multtb, mseq, 'UV')
   infile DIR!!':'!!sauce!!'.FITS'; outn=sauce; go; 
   inn sauce;  keyword='CDELT2';incl=uvtb; geth; 
   if (keyval(1)=2) then; keyval(1)=1; puth; end
$ Fix source names > 8 chars
   if (doload = 1) then; tbfind ('cl', 1);
    if (exty <> 'cl') then; outcl='multi'; go multi; outcl=uvtb
      incl 'multi';SNAMFIX ;incl=uvtb
      dozap (mdisc, sauce, uvtb, mseq, 'UV'); else
    outn '';outse 0;outcl 'multi';renam;end
   end
ret; fin

proc dosnp (tabtyp, k1, tmpl, ttype)
$ run SNPLT
$ tabtyp = table type, k1 = version #
$ tmpl = optype
$ ttype = stokes
   task 'snplt'; setuvd;
   if (doplot > 0) then 
    tbfind ('PL', ii); nplots=ntabs
    inext=tabtyp; inver=k1; source ' '; qual -1; 
    timer 0; stokes=ttype; selband -1; selfreq -1; freqid 1;bif 1;eif 0
    antennas 0; pixr 0; nplot ntel; xinc 1; optyp=tmpl; 
    opcode ''; bcount 1; xaxis 0; doebar -1; cutoff 1e-6; 
    ltype 3; dotv -1; grchan 1; go; 
    tbfind ('PL', jj);nplots=ntabs
      if (jj > ii) then;
      rplot (ii+1, jj); doextd ('PL', ii+1, jj)
      end
   end
   ret; fin


$HEAD: Types next step.
proc head(outlin); prpound; typ '#',outlin; prpound; ret;fin

$GETSCAL: Varies MSOLINT solution interval for CALIB as a
$         function of the LOOP (KK) number
proc getscal(kk,msolint)
$ if (kk<9) then
$   msolint = msol0*(1-kk/10); 
$ else
$    msolint = msol0/10;
$ end
if (kk<=2) then; msolint = 5; end
if (kk=3)  then; msolint = 3; end
if (kk=4)  then; msolint = 2; end
if (kk>=5) then; msolint = 1; end
ret;fin

$GETUVR: Varies MUVRANGE for CALIB as a function of the LOOP (KK)
$number
proc getuvr(kk,muvrange)
 if ((dotgsc > 1) & (douvr>0)) then;
   muvrange(1)=defuv; muvrange(2)=0;
   if (kk=1)  then; muvrange(1)=muvrange(1); end
   if (kk=2)  then; muvrange(1)=muvrange(1)*0.8; end
   if (kk=3)  then; muvrange(1)=muvrange(1)*0.75; end
   if (kk=4)  then; muvrange(1)=muvrange(1)*0.666; end
   if (kk>=5) then; muvrange(1)=0;end
end
ret; fin

proc runscal (kk,msolint,muvrange,mode,mdocal,mgainu,isn)
$RUNSCAL: Sets up and runs CALIB
$         Only use Failed solutions for first 80% of iterations
$         (30% for PHCAL)
$         Set up name etc before running
$         
$    kk = current run in requested series
$    msolint = current solint
$    muvrange = current uv range
$    mode = self-cal mode, 'p' or 'a&p'
$    mdocal = docal; mgainu=gainuse
$    isn = version # of SN table generated
$
 task 'calib'; head('Self-cal uv data'); setuvd; caldef; calprmdef;
 bchan=mbchan; echan=mechan; uvrange=muvrange; wtuv=mwtuv; 
 solmode=mode; docal=mdocal; gainuse=mgainu; refant=mrefant; 
 solint=msolint; aparm 3 0 0 0 0 0 1 0 1
 in2na=mname; in2cl=mapclas; in2seq=0; in2di=mdisc;
 ncomp(1)=-mnit0; smodel=0; 
    if (mname=PACAL) then 
    clr2n;ncomp=0; smodel=1,0;end;
    if ((mname=PHCAL) & (dophsc<3) & (kk=1)) then 
    clr2n;ncomp=0; smodel=1,0;end;
    if ((mname=TARGET) & (kk=1) & (dophsc<1)) then 
    clr2n;ncomp=0;aparm 3 0 0 0 0 0 1 0 1; smodel=1,0;end;
    in2d=ind
type '1'
 if (dosmod>0) then; clr2n; ncomp=0; smodel=0; smodel(1)=pkval; end
 outna=mname; outcl='SCAL'; outse=0; outdi=mdisc; 
 aparm=3,0,0,0,0,0,1,0,1,0; 
if ((mname=PHCAL) & (DOPHSC=1)) then; aparm(9)=0;aparm(7)=3;end;
if ((mname=PHCAL) & (DOPHSC>1) & (kk>0.3*nscal)) then; 
 aparm(9)=0;aparm(7)=3;end;
$
if ((dotgsc>1) & (LOOP>2) & (kk>0.7*nscal)) then
 aparm(7)=3;aparm(9)=0; 
end; 
baddisk=mbaddisk; 
 test = substr(solmode,1,1); aparm(1)=3; soltype '';
 type '2'
 if (test = 'p') then
    aparm(1)=3; soltype '';
 end
 if (test = 'a') then
    aparm(1)=4; soltype 'L1';
 end
 type '3'
 tbfind('SN', isn); isn=isn+1; snver=isn
 if (docal>0) then
    if (gainuse=0) then; gainuse=isn-1; end;
 end      
 type '4'
 typ '# Solution Interval  ',solint
 typ '# Solution Mode      ',solmode
 typ '# O/P SN table       ',isn
 type'# UV range (klambda) ',uvrange; prpound; qh;
 go;
ret; fin

$GETIM: Varies MFLUX clean limit and MNITER iteration limit for IMAGR
$       as a function of the LOOP (KK) number
proc getim(kk,mflux,mniter)
 if (kk<3) then 
    mflux=-1; mniter=mnit0;
 else 
    if (kk=3) then
       if (rms < mflx0) then
          frate = 0.0
          rmsz = mflx0
       else
          frate = (mflx0 - rms) / (nscal - 4);
          rmsz = rms;
       end;
       type 'RATE, RMS0: ',frate, rmsz;
    end;
    mniter = mnit0;
    mflux = frate * (kk-4) + rmsz;
    type 'RMS FLUX LIMIT: ', mflux;
 end
ret;fin

$RUNIMG: Set up and run IMAGR
$       Wobbles CELLSIZE back and forth between LOOPs (KKs) if dowobbl
$       set to 1
$ Cannot use IMAGPRM etc.in case of old versions
proc runimg(kk,mflux,mniter,dowobbl,mbmaj,mbmin,mbpa)
 task 'imagr'
 head('Running Imagr'); setuvd; caldef; imagdef; IFMULT
 inse=0; in2di=mdisc; outdi=mdisc; 
 if (dowobbl=1) then
    if (mod(kk,3)=1) then 
       mcells=mcells*(cellf**2)
    else
       mcells=mcells/cellf
    end 
    if (kk=nscal) then
       mcells = acell
    end
 end
 if ((inn=TARGET) & (incl='MULTSC') & (dotgsc>1)) then
   getuvr(kk,muvrange)
 end
 if ((inn=PHCAL) & (MDOTV=2)) then; dotv 1;end
 if ((inn=TARGET)&((incl='MULTSC')!(incl='SPLIT'))&(MDOTV>0)) then; 
  dotv 1;end
 
 imagrprm 0;allok 0
 cellsize=mcells,mcells; imsize=mimsize; bchan=mbchan; echan=mechan; 
 stokes=mstokes; uvwtfn=muvwtfn; gain=0.1; flux=mflux; niter=mniter; 
 rashift(1)=mshift(1); decshift(1)=mshift(2); nbox=mnbox
 clbox=mclbox;
 bmaj=mbmaj; bmin=mbmin; bpa=mbpa;baddisk=mbaddisk; outse 0; qh; go;
 if (mniter>0) then; 
 setuvd; inclass=mapclas; inseq=0; inty='MA';inver 0;go ccmrg;end
ret;fin

$GETSIZE: return ra and dec size in pixels of current map
proc getsize
 keyw='NAXIS1';geth;rasi=keyv(1)
 keyw='NAXIS2';geth;decsi=keyv(1)
ret;fin
 
$SETNEAR(BOXSI,IC) : sets a box around PIXXY(1,2) location
$                 Limits box to within the image
$ PIXXY  - X and Y pixel number for box center
$ BOXSI  - radius of region to set blc and trc
$ IC     - channel number to set
proc setnear(boxsi,ic)
 blc=max(1,pixx(1)-boxsi),max(1,pixx(2)-boxsi);getsize
 trc=min(rasi,pixx(1)+boxsi),min(decsi,pixx(2)+boxsi)
 blc(3)=ic; trc(3)=ic
ret;fin

$GETBOX  - Returns the rms inside a box at X,Y
$ RMS    - OUTPUT RMS noise of box centered on X,Y
proc getbox(x,y,boxsi,ic,rms)
 pixxy=x,y,0; setnear(boxsi,ic); imstat; rms=pixstd
ret; fin

$GETRMS: Search four corners for second lowest RMS
$ More than one corner may have the source in it.
$ RMS   - OUTPUT average of middle 2 RMS values
$         RMS is put in map header.
$ X2, Y2- Scratch variables used nowhere else
proc getrms(ic,rms)
 head('Finding rms noise in image')
 setuvd; incl=mapclas; inty='MA'; ins=0; doinv=-1;
 getsize; boxsi=rasi/10; xy=1.5*boxsi
 getbox(xy,xy,boxsi,ic,rmsb(1)); x2=rasi-xy
 getbox(x2,xy,boxsi,ic,rmsb(2)); y2=decsi-xy
 getbox(x2,y2,boxsi,ic,rmsb(3))
 getbox(xy,y2,boxsi,ic,rmsb(4))
 rms=max(min(rmsb(1),rmsb(2)),min(rmsb(3),rmsb(4)))/2
 rms=min(max(rmsb(1),rmsb(2)),max(rmsb(3),rmsb(4)))/2+rms
 keyword='RMSNOISE';keyv=rms,0;keytyp='R';puthead
 prpound; typ '# RMS =',rms; prpound
ret; fin

$GETSEQ:  GET sequence number of file
proc getseq;ins=0;keywor='imseq';geth;ins=keyv(1);ret;fin

$DELPRE:  Delete Previous map and self calibration
proc delpre(kk)
 if (kk>1) then
    setuvd; incl='SCAL'; chkname
    if (error <=0) then; getseq;
    if (ins>1) then ins=ins-1; clrst;zap; end;end
    setuvd; incl=mapclas; inty='MA'; getseq; 
    if (ins>1) then ins=ins-1; clrst;zap; end
 end
 setuvd; incl='?BM001'; inty='MA'; ins=0; clrst; 
 type 'zapping: ',inname,inclas,intype
 zap;recat
ret; fin

$POLIMG: produce images of I,Q & U through runimg
proc polimg (mbmaj, mbmin, mbpa)
  notuch=loop; loop=1; mstokes='I'
  runimg (loop,mflux,mniter,dowobbl,mbmaj,mbmin,mbpa)
  if (doplzn > 0) then
     pniter=0.1*mniter
     mstokes='Q';mapclas 'QCL001'
     runimg (loop,mflux,pniter,dowobbl,mbmaj,mbmin,mbpa)
     mstokes='U';mapclas 'UCL001'
     runimg (loop,mflux,pniter,dowobbl,mbmaj,mbmin,mbpa)
  end
  mstokes='I'; mapclas 'ICL001'
  dozap (mdisc, sauce, 'IBM001', 1, 'MA')
  dozap (mdisc, sauce, 'QBM001', 1, 'MA')
  dozap (mdisc, sauce, 'UBM001', 1, 'MA')
$$  dozap (mdisc, sauce, 'IMAGR', 1, 'UV')
$$  dozap (mdisc, sauce, 'IMAGR', 2, 'UV')
  loop=notuch
ret; fin

$POLCOM: run COMB on o/p of POLIMG
proc polcom (sauce)
   mname=sauce; class='QCL001'; seq=0; uvdisc=mdisc; setuvd
   intype 'MA'; task 'comb'; in2nam=innam; in2cl 'UCL001'; 
   blc 0; trc 0; outn=innam; aparm 0; doalign 1; bparm 0;
   clr3n;in3d ind;clr4n;in4d ind
   outcl 'POLI'; opcode 'POLI'; go
   opcode 'POLA'; outcl 'POLA'; go
   dozap (mdisc, PACAL, 'QCL001', 1, 'MA')
   dozap (mdisc, PACAL, 'UCL001', 1, 'MA')
ret; fin

proc pltmap (kk,mblc, mtrc, nn)
$ generate and plot the kntr map
$ kk = 0, use kntr & plot stokes I
$ kk = 1, use kntr & plot full polzn
$ mblc, mtrc - obvious
$ nn = 1, special override for parms for 3C286 polzn
    setuvd; intype 'MA'; tbfind ('PL',ii); nplots=ntabs
    notuch=loop; getrms(1,rms)
if (doplot > 0) then
    dotv -1
    blc=0; trc=0; blc=mblc; trc=mtrc; imstat; plev 0; ltype 3;
    cbplot 2;docont 1;dogrey -1;dovect -1;clev=3.0*rms;zinc=0;pixr=0
    if (rasi*decsi > 1024*1024) then; clev=5.0*rms;end
    levs -4,-2,-1,1,2,4,8,16,32,64,128,256,512,1024; 
    if (kk = 0) then; 
    clr2n;in2d ind;clr3n;in3d ind;clr4n;in4d ind;go kntr; end
    if (kk > 0) then
       blc=0; trc=0; dovect=1; rotate=0;clr2n;in2d ind;
       incl 'POLI'; imstat; incl 'ICL001'; dovect=1
       in3cl 'POLI'; in3n=innam; in3d=ind;
       in4cl 'POLA'; in4nam=innam; in4d=ind; 
       pcut=clev/2; icut=2*clev; xinc 2; yinc 2; 
       factor=(5/pixval)*acell/0.015; pixxy 0
       if (nn=1) then; factor=150; pcut=0.1; ; end
       blc=mblc; trc=mtrc; go kntr
    end
    loop=notuch; tbfind ('PL',jj); nplots=ntabs
    if (jj > ii) then; 
       rplot (ii+1, jj); doextd ('PL', ii+1, jj)
       end
    blc 0; trc 0
end
ret; fin

proc MKLOFLAG
task 'TABGET'
mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc
SETUVD
GETANTNO
inext 'an';inver 1
anten 0;base 0
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'MK2' then antenn i,0
end;end
for i=1:NTEL;pixxy i 1 1;tabget;if keystr = 'LOVELL' then basel i,0
end;end
if ((antenn(1) > 0) & (basel(1) > 0)) then
task 'UVFLG'; source ''; infile ''; timer 0; bchan 1; echan 0; 
bif 1; eif 0; aparm 0;stokes ''; opcode 'FLAG'; flagver 1; go
end
anten 0;base 0;opcod ''
fin
proc FNDAFLUX
$
 if ((ACAL = 'OQ208') ! (ACAL = '1404+286')) then; 
  if (oband='L') then; GFLUX = 1.0; end
  if (oband='C5') then; GFLUX = 2.4; end
  if (oband='C6') then; GFLUX = 2.4; end
 end
 if ((ACAL = '0552+398') ! (ACAL = 'DA193')) then; 
  if (oband='L') then; GFLUX = 2.2; end
  if (oband='C5') then; GFLUX = 5.9; end
  if (oband='C6') then; GFLUX = 5.9; end
 end
 if (ACAL = '2134+004') then; 
  if (oband='L') then; GFLUX = 4.5; end
  if (oband='C5') then; GFLUX = 8.7; end
  if (oband='C6') then; GFLUX = 8.9; end
 end
fin

$------------------------------------------------------
$ Begin 

proc merlin
$
$ Check inputs
$
   CHECKIN;
$
if (DIR='') then; DIR='PWD';end

$
   source ''
   if (target = ' ') then; 
         mess ('You must give the target name');
         source(1)='stop'
         end
   if ((dophsc >0) & (phcal = ' ')) then
         mess ('You must give phcal source if dophsc >0');
         source(1)='stop'
    end
         if (ACAL = ' ') then
         mess ('You must give a point source cal ACAL')
         source(1)='stop'
      end
   calok = 0;
   if ((ACAL='OQ208')!(ACAL='1404+286')!(ACAL='2134+004')) then
      calok = 1; end
   if ((ACAL ='0552+398')!(ACAL='DA193')) then
      calok = 1; end
  if (AFLUX = 0 ) then
         mess ('Should you set AFLUX for point cal ACAL?')
         if (calok > 0) then
          type 'Or, at L or C band defaults can be used for'
          type 'OQ208, 0552+398 or 2134+004, see HELP MERLIN'
          if (DOFLUX < 1) then
             mess ('Amplitude calibration will be poor')
             type 'You should set DOFLUX > 0 &/or set AFLUX'
          end
          type 'To stop now type 1; to continue press return'
          giveup
         else; 
           mess ('Amplitude calibration will be poor')
           type 'You must set AFLUX for ACAL=' ACAL 
           end
           type 'To stop now type 1; to continue press return'
         giveup
   end

    if (source(1)='stop') then;
         mess ('Source names must be exact MERLIN names')
    inn source(1);imh;end
       source(1) ''
if (PACAL = '') then; PACAL='3c286';
mess ('Setting default PACAL 3C286');
end
$
$---------------- # 1, load the data --------------------------------
   
$
$ read the FITS files using FITLD
$
   if (doload > 0) then
$
$ Stop if files exist (?change to question?)***
inn '';incl '';inty '';inse 0;chkn
if (error < 1) then;

prexcl
type 'If you want to load FITS files (doload >0), please empty'
type 'your AIPS directory. Pre-existing files may may be altered,'
type 'tasks may fail or you may get bad results.'
type 'To stop now type 1; to continue press return'
giveup
end
source (1) ''
$
      mess ('Loading and organizing the data')
$      version = tver; 
      task 'FITLD'; outdisk=mdisc; outcl=uvtb; 
      outseq=mseq; douvc -1; optype ' '; 
      sort 'TB'; rotate 0; incl=outcl; inseq=outs; ind=outd; 
      source ' '; aparm 0; aparm(1)=0.5;
         if (ACAL = ' ') then
         mess ('You must give point source cal ACAL')
      end
  if (AFLUX = 0 ) then
         mess ('Should you set AFLUX for point cal ACAL?')
         if (calok > 0) then
          type 'Or, at L or C band defaults can be used for'
          type 'OQ208, 0552+398 or 2134+004, see HELP MERLIN'
          if (DOFLUX < 1) then
             mess ('Amplitude calibration will be poor')
             type 'You should set DOFLUX > 0 &/or set AFLUX'
          end
          type 'To stop now type 1; to continue press return'
          giveup
         else; 
           mess ('Amplitude calibration will be poor')
           type 'You must set AFLUX for ACAL=' ACAL 
           end
           type 'To stop now type 1; to continue press return'
         giveup
   end
$
      if (doload >= 1) then
         mess ('If message server says FITLD1: ERROR...')
         mess ('check DIR and DOLOAD ')
         if (TARGET <> ' ') then
            load (TARGET)
            if (doload=2) then
               innam=TARGET; incl=uvtb; inseq=mseq; ind=mdisc;
               outd=indisk; outn=TARGET; outcl='MULTTB'; outs=mseq; 
               rename
            end
         mess ('If message server says FITLD1: ERROR...')
         mess ('check DIR and DOLOAD ')
         end
      end
$ load all files individually
      outcl 'uvtb'
      if (doload < 2) then
         if (PHCAL <> ' ') then
            load (PHCAL)
         end
         if (PACAL <> ' ') then
            load (PACAL)
         end
         if (ACAL <> ' ') then
            load (ACAL)
         end
$ check for existence 
$         dozap (mdisc, TARGET, 'MULTTB', mseq, 'UV')
         wname(1)=TARGET; wname(2)=PHCAL; wname(3)=ACAL; wname(4)=PACAL
         for ii=1 to 4;
            exist (mdisc, wname(ii), 'multi', mseq, wexist(ii));
         end; 
$
$ concatanate
         ind=mdisc; inseq=mseq; incl='multi'; in2d=ind; in2seq=inseq;
         in2cl='multi'; reweight=0; dopos=-2; doarray=1; outn='DBCON';
         outcl='TEMP'; outs=0; outd=mdisc; jj=0; nn=0; task 'DBCON'
         for ii=1 to 4;
            kk = 0;
            if (wexist(ii)) then
               if (jj=0) then
                  inname=wname(ii); kk = 1; nn = nn + 1;
               end;
               if (jj=1) then
                  in2name=wname(ii); kk = 0; nn = nn + 1; go
                  inname='DBCON'; incl='TEMP'; inseq=0;
               end;
               jj = jj + kk;
            end;
         end;
$
$ special case: only one source
         if (nn = 1) then
            userid=0; outname='DBCON'; outcl='TEMP';
            outdisk=mdisc; outseq 0; rename;
         end;
$
$ rename, tidy up.
$
         indisk=mdisc; innam='DBCON'; inclass='TEMP'; inseq=0; 
         outd=indisk; outn=TARGET; outcl='MULTTB'; outs=mseq; rename
$$$*** should this be multi?
         dozap (mdisc, TARGET, 'multi', mseq, 'UV')
         dozap (mdisc, PHCAL, 'multi', mseq, 'UV')
         dozap (mdisc, ACAL, 'multi', mseq, 'UV')
         dozap (mdisc, PACAL, 'multi', mseq, 'UV')
         allzap (mdisc, 'DBCON', 'TEMP', 1, 3, 'UV')
      end
   end
$
   if (doload <= 0) then
$ Stop if files exist (?change to question?)***
inn '';incl '';inty '';inse 0;chkn
if (error < 0) then;
prexcl
type 'If doload <=0, should be one (MULTTB) file only in aips;'
type 'remove others please. If you continue, these other files may'
type 'be altered, tasks may fail or you may get bad results.'
type 'To stop now type 1; to continue press return'
giveup
end
source (1) ''

      ind=mdisc; innam=TARGET; incl='MULTTB'; inseq=mseq; 
      doextd ('SN', -1,-1); doextd ('CL', 2, 20)
      doextd ('BP',-1,-1)
   end
mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc;setuvd
keyword 'sortord';keystr '';keyval 0;geth
 if (keystr <> 'tb') then;outn mname;outcl class;outse (seq+1)
 sort 'tb';rotate 0;go UVSRT;zap;inse (seq+1);outse seq;renam
end

$declare string*1 oband
$declare scalar gflux

FNDAFLUX
if (AFLUX = 0) then; AFLUX = GFLUX;end

DEFREFANT
GETBAND
$
  if (oband='P') then; acell=0.180; defuv=50; srchcl=0.3; end
  if (oband='L') then; acell=0.040; defuv=170; srchcl=0.05; end
  if (oband='C5') then; acell=0.015; defuv=1000; srchcl=0.018; end 
  if (oband='C6') then; acell=0.01; defuv=1000; srchcl=0.012; end 
$ Set smaller cells for 6-7 GHz
  if (oband='K') then; acell=0.003; defuv=5000; srchcl=0.004; end
$ index, set weights to 10, set pt src flux
$
   ind=mdisc; innam=TARGET; incl='MULTTB'; inseq=mseq; cparm 0 0 0.5;
   infile '';
go indxr; 
if (doload = 1) then
 task 'TAMRG';inext 'cl';inver 1;outver 1;aparm 1 0
 bparm 1 4 0;cparm 0;dparm 0;go TAMR
end 
  outn=inn; outcl=incl; outs=ins; outd=ind; aparm 0; antwt=10; 
go wtmod; 
SETWTS
   source=''; source(1)=ACAL; bif 0; eif 0; optyp '';  calcode ''; 
   sysvel 0; restf 0; veltyp ''; veldef ''; aparm 0; zerosp 0; 
   zerosp(1)=aflux; go setjy; 
   isline(line,fbchan,fechan,kk)
$ *** why kk above not gbnd?
   if (line=1 & restfreq (1) > 0 & restfreq(2) > 0) then
      aparm(1) = ceil(fechan/2); source(1)=TARGET; sysvel=trgvel
      restfreq=trgfreq; veltyp 'LSR'; veldef='RADIO'; zerosp 0;
      go setjy; aparm 0; sysvel 0; restf 0; veltyp ''; veldef ''; 
   end
   optyp 'scan'; docrt -1; source ' '; 
   plname = TARGET!!'_LISTR'
   if (doplot > 0) then; rprint; go listr; docrt 0;end
$
$ Flag Lo-Mk if present
MKLOFLAG
$
$ automatic flagging of points > 10 sigma from rms
$
ISLINE (line, fbchan, fechan, gbnd) 
$
   if (doedit > -1) then;
   task 'vplot'; caldef; clr2n; 
   source ACAL PHCAL PACAL; bchan=mbchan; echan=mechan;
   doarray 0; solint 3; bparm 0; aparm 0; aparm 1,10,0; optyp ''
$ in 31DEC03 on unix flagver must be 0
$   TBFIND('fg',0);flagver=ntabs;go;
$   source TARGET; TBFIND('fg',0);flagver=ntabs;bchan=scb;echan=sce;go
   flagver 0
   end
$
   if (doedit > 0) then
      mess ('Editing')
      type 'Important: For each source:'
      type '"FLAG TIMERANGE" will flag all sources,'
      type 'other FLAGS e.g. "BELOW" flag current source only.'
      type 'Set "STOKES FLAG" (1111 recommended)'
      type 'Usually, toggle "1-Ch" to flag all channels'
      innam=TARGET; incl='MULTTB'; ind=mdisc; ins=mseq
      task 'ibled'; caldef; clr2n; docat -1; dohist -1; clr3n;
      docal 1; gainuse 0; TBFIND('fg',0);flagver=ntabs
      bchan=mbchan;echan=mechan;dparm 0; source ''; nmaps 0
      if (PACAL <> ' ') then
         source(1)=PACAL; type 'Editing: ',PACAL; go
      end
      if (ACAL <> ' ') then
         source(1)=ACAL; type 'Editing: ',ACAL; go
      end
      if (PHCAL <> ' ') then
         source(1)=PHCAL; type 'Editing: ',PHCAL; go
      end
      if (TARGET <> ' ') then
         bchan=scb;echan=sce
         source(1)=TARGET; type 'Editing: ',TARGET; go
      end
   end
$
$--------------------------------------------------------------------
$
$
$   dozap (mdisc, ACAL, 'SPLIT', 1, 'UV')
$
$
$ self-calibrate ACAL and 3C286 if present
       mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc
       task 'calib'; calsour(1)=ACAL; calsour(2) '3C286'
       caldef; calprmdef
       bchan=mbchan;echan=mechan;cparm=0
       solint 0.5; flagver=0; snver 1; aparm(7)=5; go
$
$ CL1 + SN1 -> CL2 ACAL phase
       mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc
       plname = acal!!'_'!!pacal!!'_p_sn'
       source(1) ACAL
       dosnp ('SN', 1, 'PHAS', ' ')
       runcl (1, 1, 2, 'self', '', ' ')
       tget calib;docal 1;gainuse 2;aparm(1) 4;solint 2
       solmo 'a&p';snver 2;go
$ CL2 + SN2 -> CL3 ACAL a&p
       mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc
       source(1) ACAL
       plname = acal!!'_'!!pacal!!'_a_sn'
       dosnp ('SN', 2, 'AMP', ' ')
       runcl (2, 2, 3, 'self', '', 'ampl')

$ Use point source calibrator to determine BP corrections
   if (mdoband > -1) then
     innam=TARGET; incl='MULTTB'; ind=mdisc; ins=mseq; 
     isline(line,ii,fechan,gbnd)
$ ** why ii not fbchan above?
     if (gbnd=1) then
$
       mess ('Bandpass calibration using point source')
       task 'bpass'; calsour ''; calsour(1)=ACAL; caldef; clr2n; 
       clr3n;smodel 0; docal 1; gainuse 3; solint -1; bpver 1; 
       flagver 1; refant=antref; bpassprm 0; bpassprm(5)=1; 
       bpassprm(10)=3;bpassprm(11)=1; chansel 0; antwt 0;
       ichansel mbchan, mechan,1;go
$
$ examine BP solution using POSSM
  if (doplot > 0) then; 
       task 'possm'; ind=mdisc; innam=TARGET; incl='MULTTB'; 
       ins=mseq; tbfind ('PL', ii); nplots=ntabs
       source ''; source(1) ACAL; timerang=0; stokes 'i'; bif 1; 
       eif 0;bchan mbchan; echan mechan; docalib 1; flagver 1; doband=1
       bpver 0; smooth 0; bparm 0; aparm 0; aparm 1,1,1e-10,10,-30,30,0
       aparm(4)=2 * aflux; antennas=antref,0; baseline 0; solint 0
       stokes 'i';dotv -1; nplots (ntel-1); 
go POSS
       tbfind ('PL',jj); nplots=ntabs
       if (jj > ii) then; 
$       plname = acal!!'_chan'
        plname = acal!!'_bandpass'
          rplot (ii+1, jj); doextd ('PL', ii+1, jj)
          end
$       doextd ('SN',1,1); doextd('CL',2,2)         
     end
   end
end
$
$
$ Refine flux scale at C and L band
  if ((DOFLUX = 1) & ((oband='L') ! (oband='C5') ! (oband='C6'))) then
    SET3C
     if (m>0) then
     tget CALIB;keyword 'CRVAL3';keystrng '';GETH
      uvrang(1) 0;uvrang(2)=(KEYVAL(1)+KEYVAL(2))*0.00000007
      solint 10;clro;outd ind; bchan mbchan; echan mechan; 
      flagver 0;clr2n;in2d ind;anten 0;cparm 0;snver 3;docal 1;
      gainuse 2; aparm 2 0 0 0 0 0 3;solty '';solmo 'a&p';
      refant mrefant;doband=gbnd;calso ACAL '3C286';go CALIB
     task 'GETJY';calso '3C286''';source ACAL'';snver 3
     go GETJY
$
$ Apply flux scale to all sources
     tget CALIB;uvrang 0;aparm(1) 4;snver 4;go CALIB
     task 'SNCOR';source '';
     opcod 'ZPHS';snver 4;anten 0;stokes '';go SNCOR
     inext 'cl';inver 3;extd;inver 2;extd
     tget CLCAL;OPCODE='CALP';snver 4;gainver 1;gainuse 2;source ''
     interpol '';
     calso '';timer 0;anten 0; SAMPTYPE 'mwf';bparm 1000 0;go CLCAL
     end
      else
      inext 'cl';inver 3;extd;inver 2;extd
      inver 1;outver 2;ncount 1;keywor '';keyval 0;keystr ''
      go tacop
    end
$
$ After the next tidy-up CL2 either contains amplitude refinement
$ or is same as CL1 
$ --------------------------------------------------------------------
if (dophsc > 0) then
   mess ('Phase reference calibration')
$
$ Calib on phase calibrator -> SN1
$
$ delete any phase-cal image & temporary files
$
   allzap (mdisc, PHCAL, 'ICL001', 1, 3, 'MA')
   allzap (mdisc, PHCAL, 'SCAL', 1, 3, 'UV')
   dozap (mdisc, PHCAL, 'TCAL', 1, 'UV')
$
$ delete pre-existing SN tables & any CL table > 2
$ ** changed to allow for amp cal
$
   mname=TARGET; class='MULTTB'; uvd=mdisc; seq=mseq; setuvd
   task 'TASAV';outn ACAL; outcl '';outse 0;go TASAV
   tbfind ('SN', ii); if (ii > 0) then; doextd ('SN', 1, ii); end
   tbfind ('CL', ii); if (ii >=3) then; doextd ('CL', 3, ii); end
$
$ if polzn and L band perform ionospheric correction
$
$  if ((doplzn > 0) & (oband = 'L')) then
$     task 'farad'; infile ''; source ''; timer 0; antennas 0
$     subarray 1; gainver 1; opcode 'mod'; bparm 0; bparm(2)=sunspot;
$     go
$  end
$
$ split and average to make this step faster
$
   mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd
   isline (line,fbchan,fechan,gbnd);
   runspl (PHCAL,'TCAL',1,gbnd,2,-1,mbchan,mechan)
   mname=PHCAL; class='TCAL'; uvdisc=mdisc; seq=mseq; setuvd
$
   msolint = 1; loop = 1; muvrange = 0; mode='P'; mdocal=-1; mgainu=0
  mflux = 0; mniter = 1000; dowobbl = -1; mcells=acell;
  mimsiz = 512; mclbox = 251, 252, 260, 261
    if (dophsc > 1) then; mimsize = timsize; mclbox = 0; end
   mnbox = 1; mstokes 'I'; mapclas 'ICL001'
$
   nscal = 1
   if (dophsc = 1) then; nscal = 2; end
   if (dophsc > 1) then; nscal = 3; end
   prpound
   for loop = 1 to nscal by 1
      prpound; print '# LOOP',loop,'OF',nscal,'#'; prpound;
      mname = PHCAL; class = 'TCAL'; seq = mseq; uvdisc=mdisc;
      calsour(1)=PHCAL
      if ((dophsc > 0) & (loop=nscal)) then
         mode='A&P'; mdocal=-1; mgainu=0; msolint=5
         class 'SCAL'; seq = 2; 
         if (dophsc = 1) then; seq = 1; dosmod = 1; end
      end
$
      if (dophsc = 3 & loop = 1) then; 
      if (DIR='') then; DIR='PWD';end
      task 'FITLD';inf DIR!!':'!!PHCAL!!'.ICL001.FITS'
      outd uvdisc; opty '';source '';bch 0;ech 0;opcod ''
      outn PHCAL;outcl 'icl001';outse 1;go fitl
      inn PHCAL;incl 'icl001';inse 0;inver 0;outver 0;go CCMRG
      dosmod = -1;end
$ 
     runscal(loop,msolint,muvrange,mode,mdocal,mgainu,isn)
     dosmod = -1
$
       plname = phcal!!incl!!char(inse)!!'_'!!char(isn)!!'_p_sn'
      dosnp ('SN', isn, 'PHAS', ' ')
      if ((dophsc > 0) & (loop=nscal)) then
       plname = phcal!!incl!!char(inse)!!'_'!!char(isn)!!'_a_sn'
         dosnp ('SN', isn, 'AMP', ' ')
         end
$
      if (dophsc = 3 & loop = 1) then; 
      inn PHCAL;incl 'icl001';inse 1
      outn inn;outcl 'iclold';outse inse;rena;end
$
      if (loop=2) clrtemp
$
$ copy -> SN_loop on TARGET then CL_loop + SN_loop -> CL_loop+1 
$
      coptb (TARGET,'MULTTB',mdisc,1,'SN',isn,loop)
      jsn = loop + 2; kk=2; smmode='PHAS'
      if ((dophsc > 0) & (loop=nscal)) then 
$
        kk=loop+1; smmode='AMPL'; 
         if ((dophsc = 1) & (loop = nscal)) then; kk=4; jsn=5; end
      end
      mname=TARGET; class='MULTTB'; seq=1; uvdisc=mdisc
      if  ((dophsc > 0) &(loop = nscal)) then
        setuvd;task 'TASAV';outn PHCAL; outcl '';outse 0;go TASAV
      end  
      runcl (loop,kk,jsn,'2PT','',smmode)
      if ((dophsc = 1) & (loop = 1)) then
         jsn = 4; runcl (loop,kk,jsn,'2PT','',smmode)
      end
$ 
      mname=PHCAL; class='SCAL'; seq=loop; uvdisc=mdisc
      mimsiz=512; mclbox = 251, 252, 260, 261 
      if (dophsc > 1) then; mimsize = timsize; mclbox = 0; end
      runimg(loop,mflux,mniter,dowobbl,0,0,0)
$
      delpre(loop);
$
      mname=PHCAL; class='ICL001'; seq=loop; uvdisc=mdisc;setuvd
      inty 'ma';getsize;
      mblc = rasi*0.05,decsi*0.05, 1; mtrc = rasi*0.95,decsi*0.95, 1;
      plname = PHCAL!!class!!char(seq)
      pltmap(0,mblc,mtrc,0)
      blc=mblc; trc=mtrc;imst;blc 0;trc 0
      if ((dophsc = 1) & (loop = 1)) then; pkval=pixval; end
$
      print 'End of loop ', loop
   end
   dozap (mdisc, PHCAL, 'TCAL', 1, 'UV')
   allzap (mdisc, PHCAL, 'SCAL',1, 3, 'UV')
end
$
$
$--------------------------------------------------------------------
   if (doplzn > 0) then
   mess ('Polarization calibration')
if (dophsc <= 0) then
   mess ('Must have DOPHSC > 0 if doplzn > 0')
   source(1) 'stop';inn source(1);imh;inn ''
end

      dozap (mdisc, PHCAL, 'TCAL', 1, 'UV')
      dozap (mdisc, PACAL, 'TCAL', 1, 'UV')
$ delete pre-existing SN tables
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd
$      tbfind ('SN', ii); if (ii > 0); then; doextd ('SN', 1, ii); end
      doextd('SN', 1, 3); doextd ('CL', 6, 6)
$
$  split out the phase cal & create a multi source file
$
      isline (line,fbchan,fechan,gbnd);
      runspl (PHCAL,'TCAL',1,gbnd,2,-1,mbchan,mechan)
      mname=PHCAL; class='TCAL'; uvdisc=mdisc; seq=mseq; setuvd
      outn inn; outcl 'MULTI'; outd=mdisc; outs=mseq; aparm 1,0; 
      go multi; 
      inn outn;incl outcl;cparm 0 0 0.5;infile '';go indxr;SNAMFIX
$  
$  get map peak from phase-cal map & run PCAL
$
      mname=PHCAL; class='ICL001'; seq=0; uvdisc=mdisc; setuvd
      inty 'ma';getsize;
       imh
      mblc = rasi*0.25,decsi*0.25, 1; mtrc = rasi*0.75,decsi*0.75, 1
      blc = mblc;trc= mtrc;imstat; class='MULTI'; seq=mseq
      setuvd; source ''; source(1)=PHCAL; optyp ''; zerosp(1)=pixval;
      go setjy; task 'pcal'; docalib -1; clr2n; ncomp 0; pmodel 0; 
      cparm=0; solint=2; soltyp='appr'; inver 0;bparm 0;doband -1
      pmodel(1) zerosp(1);refant=antref; antennas=0; go pcal; 
      prtask 'pcal'; 
      if (doplot > 0) then; 
      plname = TARGET!!'_PCAL'
      rprint; prtmsg;end; prtask ''
$getn 1
$imst
$
$  copy polzn solutions to main file
$
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd
      doextd('AN', 1, 1)      
      mname=PHCAL; class='MULTI'; uvdisc=mdisc; seq=mseq; 
      coptb (TARGET,'MULTTB',mdisc,mseq,'AN',1,1)
      dozap (mdisc, PHCAL, 'TCAL', 1, 'UV')
      dozap (mdisc, PHCAL, 'MULTI', 1, 'UV')
$
$  calibrate 
$
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd
      isline (line,fbchan,fechan,gbnd); 
      runspl (PACAL,'TCAL',-1,gbnd,2,1,mbchan,mechan)
      mname=PACAL; class='TCAL'; uvdisc=mdisc; seq=mseq; 
      loop=1; msolint=1; muvr=0; mode='P'; mdocal=-1; mgainu=0;
      runscal(loop,msolint,muvrange,mode,mdocal,mgainu,isn)
$
$ copy -> SN_isn to TARGET then CL_2 + SN_isn -> CL_6
$
      plname = pacal!!'_p_sn'
      dosnp ('SN', isn, 'PHAS', ' ')
      coptb (TARGET,'MULTTB',mdisc,mseq,'SN',isn,1)
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
      runcl (1, 2, 6, '2PT', '', '')
      dozap (mdisc, PACAL, 'TCAL', 1, 'UV')
      clrtemp
$
$ split & map PACAL in I,Q & U
$
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
      mcells=srchcl; mflux=0; mniter=1000; mimsize=512; 
      mclbox 128, 129, 382, 383;dowobbl=-1
      runspl(PACAL,'TCAL',1,gbnd,2,1,mbchan,mechan)
      mname=PACAL; class='TCAL'; uvdisc=mdisc; seq=mseq; 
      doband 0;polimg(0,0,0)
      dozap (mdisc, PACAL, 'TCAL', 1, 'UV')
$      delpre(1)
$
$ form POLI, POLA & extract information
$
      polcom (PACAL)
      incl 'POLI'; blc 252 253 1; trc 260 261 1; imstat
      incl 'POLA'; imval
      dozap (mdisc, PACAL, 'ICL001', 1, 'MA')
      dozap (mdisc, PACAL, 'POLI', 1, 'MA')
      dozap (mdisc, PACAL, 'POLA', 1, 'MA')
$
$ correct CL 5 on TARGET.MULTTB for R-L rotation, also do CL 4
$ so can correct other sources
$
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd
      task 'clcor'; caldef; opcode 'polr'; stokes 'L'; gainver 6;
      gainuse gainver; clcorp(1)=2*(33-pixval); 
      go clcor; opcode 'PHAS'; gainver 5;gainuse gainver;
      go clcor; stokes ''
$
$ calibrate and map again
$
      mcells=srchcl;
      runspl (PACAL,'TCAL',1,gbnd,2,1,mbchan,mechan)
      mname=PACAL; class='TCAL'; uvdisc=mdisc; seq=mseq; 
      mniter 1000;doband 0;polimg (0,0,0)
$      delpre(1)
      dozap (mdisc, PACAL, 'TCAL', 1, 'UV')
      dozap (mdisc, PACAL, 'SCAL', 1, 'UV')
      polcom(PACAL)
      class='POLI'; mblc=64,64,1; mtrc=448, 448,1; 
$      pltmap (0, mblc, mtrc, 0)
      plname = PACAL!!class!!'_Pol'
      pltmap (1, mblc, mtrc, 1)
$
$ map phase cal in polzn
$
      mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; setuvd 
      doextd ('CL',6,6); mcells=acell;
      runspl(PHCAL,'TCAL',1,gbnd,2,1,mbchan,mechan)
      mname=PHCAL; class='TCAL'; uvdisc=mdisc; seq=mseq; 
      mniter 1000;doband 0;polimg(0,0,0)
$      delpre(1)
      dozap (mdisc, PHCAL, 'TCAL', 1, 'UV')
      polcom(PHCAL)
      inn PHCAL; incl 'POLA';inse 0;getsize;
      mblc = rasi*0.05,decsi*0.05, 1; mtrc = rasi*0.95,decsi*0.95, 1
$      pltmap (0, mblc, mtrc, 0)
      plname = PHCAL!!class!!char(seq)!!'_pol'
      pltmap (1, mblc, mtrc, 0)
$      class 'POLI'; seq=0; pltmap (0, mblc, mtrc, 0)
      dozap (mdisc, PHCAL, 'QCL001', 1, 'MA')
      dozap (mdisc, PHCAL, 'UCL001', 1, 'MA')
   end
$      
$--------------------------------------------------------------------
   mess ('Initial imaging of the target source')
      dozap (mdisc, TARGET, '1CAL', 1, 'UV')
      dozap (mdisc, TARGET, 'IIM001', 1, 'MA')
      dozap (mdisc, TARGET, 'ICL001', 1, 'MA')
      dozap (mdisc, TARGET, 'QCL001', 0, 'MA')
      dozap (mdisc, TARGET, 'UCL001', 0, 'MA')
      dozap (mdisc, TARGET, 'IBM001',1, 'MA')
$$      dozap (mdisc, TARGET, 'IMAGR',0, 'UV')
      dozap (mdisc, TARGET, '2CAL', 1, 'UV')
      dozap (mdisc, TARGET, 'SCAL', 0, 'UV')
$
$  average and image - create multi file
$
$ Skip first imaging if no phase-ref
$
      mname=TARGET; class='MULTTB'; seq=mseq; uvdisc=mdisc; setuvd
      doextd('SN',1,5); isline (line,mbchan,fechan,gbnd); 
      v1=-1; if (doplzn>0) then; v1=1; end
      runspl (TARGET,'',1,gbnd,2,v1,scb,sce)  
      task 'MULTI' 
      inn=TARGET; inclass='SPLIT';source '';inse 0;aparm 0.5 0
      ind uvdis;outn inn;outcl 'multsc';outse 0;outd ind;go MULTI
      incl 'MULTSC';cparm 0 0 0.5;infil '';go INDXR;SNAMFIX
      dozap (uvdisc, TARGET, 'SPLIT', 0, 'UV')
      mname=TARGET; class='MULTSC';seq=0; uvdisc=mdisc; setuvd
$ 
      loop=1; mflux=0; mniter=5000; mimsize=timsize; dowobbl=-1;
      mbchan=1; mechan=1; mcells=acell; mnbox=0; mclbox=0
      nscal=loop; kk=loop;
      doband -1
      if ((dosearch>0) & (dophsc > 0)) then 
         mess ('Searching for peak in image')
         mimsize=2048; mcells=srchcl; mniter=0
         source(1)=TARGET;
         runimg (loop,mflux,mniter,dowobbl,0,0,0)
         class='IIM001'; seq=1; setuvd; intype 'MA';
         blc mimsize(1)*0.05,mimsize(1)*0.05,1; 
         trc mimsize(1)*0.95,mimsize(1)*0.95,1; imstat
         mshift(1)=(pixxy(1)-mimsize(1)/2)*-srchcl
         mshift(2)=(pixxy(2)-mimsize(1)/2)*srchcl
         dozap (mdisc,TARGET,'IIM001',0,'MA')
         dozap (mdisc,TARGET,'IBM001',0,'MA')
         mimsize=timsize; mcells=acell; mniter=5000
         mname=TARGET; class='MULTSC'; seq=mseq; uvdisc=mdisc; setuvd
      end
$
$$$$$$$$$$$$$$$$$$$$
      if (dotgsc >= 1) then
         loop=1; mflux=0; mniter=1000; dowobbl=-1;mimsiz=timsiz 
         mname=target;seq 0;docal 1;gainuse 0;class 'MULTSC'
         source(1)=TARGET;doband -1
         if (dophsc>0) then
            runimg (loop,mflux,mniter,dowobbl,0,0,0)
$print 'uvran' muvrang uvrang dotgsc defuv douvr 
         delpre (loop);end
$
         mess ('Self-cal of the target source')
$
         mname=TARGET; class='MULTSC'; seq=mseq; uvd=mdisc; setuvd
         tbfind ('SN', ii); if (ii > 0) then; doextd ('SN', 1, ii); end
$*** changed from loop=2 kk=loop-1
         loop=1; msolint=5; mode='P'; mdocal=1; mgainu=0;kk=loop 
         getuvr(loop,muvrange); calsour(1)=TARGET
         runscal(loop,msolint,muvrange,mode,mdocal,mgainu,isn)
         plname = target!!char(loop)!!'_p_sn'
         dosnp ('SN', 1, 'PHAS', ' ')
         runcl(1,1,2,'2pt','','')
         coptb (TARGET,'MULTTB',mdisc,1,'SN',1,0)
         tbfind ('SN', ii); if (ii > 0) then; doextd ('SN', 1, ii); end
$
      end
$
      if (dotgsc > 1) then
$
         if (dotgsc = 2)  then; nscal = 5; end
         if (dotgsc >= 3) then; nscal = 10; end
         cellf=1.02; dowobbl=1;
$
$ repeat self-calibration & mapping loop a number of times
$
$$$$$$$$$$$$ *** check loop, kk nos etc.
 for loop = 2 to nscal; 
    mname=TARGET; class='MULTSC'; seq=mseq; uvdisc=mdisc
    kk = loop; jj = nscal; 
    prpound; print '# LOOP',kk,'OF',jj,'#'; prpound;
$
$ Set MSOLINT & MUVRANGE as a function of LOOP
    getscal(loop,msolint); 
    getuvr(loop,muvrange);
$
$ Do Map/Clean Iteration
    mname=target;seq 0;docal 1;gainuse 0;class 'MULTSC'
    source(1)=TARGET;mimsize=timsize;doband -1
       if (loop <= nscal) then; 
         runimg(loop,mflux,mniter,dowobbl,0,0,0)
       end
$ Set MFLUX, MNITER as function of LOOP
       getim(loop,mflux,mniter);
$
$ Clean up after CALIB and IMAGR
       delpre(loop)
       mname=target;seq 0;docal 1;gainuse 0;class 'MULTSC'
       setuvd 
$
$    Do Self-Calibration
       mode='P';mgainu=2;mdocal=1;
       if ((dotgsc > 3) & (loop=nscal-2)) then
          mode='A&P';mgainu=3; msolint=5
       end
       if ((dotgsc > 3) & (loop=nscal-1)) then
         mode='A&P';mgainu=3; msolint=3
       end
       if ((dotgsc > 3) & (loop=nscal)) then
         mode='A&P';mgainu=3; msolint=2
       end
       tbfind ('SN', ii); if (ii > 0) then; doextd ('SN', 1, ii); end
       calsour(1)=TARGET
       runscal(loop,msolint,muvrange,mode,mdocal,mgainu,isn)
       plname = target!!char(loop)!!'_p_sn'
       dosnp ('SN', 1, 'PHAS', ' ')
       if (mode='p') then; runcl(1,2,3,'2pt','','');end
       if (mode='A&P') then
         plname = target!!char(loop)!!'_a_sn'
         dosnp ('SN', 1, 'AMP', ' ')
         runcl(1,3,4,'2pt','','ampl')
       end
$ Copy back appropriate SN's  
       if ((dotgsc>3) & (LOOP=nscal-3)) then      
       coptb (TARGET,'MULTTB',mdisc,1,'SN',1,0);end
       if ((dotgsc>1) & (LOOP=nscal)) then      
       coptb (TARGET,'MULTTB',mdisc,1,'SN',1,0);end
       tbfind ('SN', ii); if (ii > 0) then; doextd ('SN', 1, ii); end
$
     if ((nscal>5) & (loop=5)) clrtemp
         end
clrtemp
      recat
      end
$$$$$$$$$$$$$ After self-cal
 if (loop >= nscal) then
$ apply SN tables if any to MULTTB
 if (dotgsc>0) then;
   mname TARGET; class 'MULTTB';seq 1;setuvd
   if (dophsc > 0) then 
    runcl (1,5,6,'2pt','','');end
   if (dophsc <1) then 
    runcl (1,2,6,'2pt','','');end
   if (dotgsc>1) then
   runcl (2,6,7,'2pt','','');end
   if (dotgsc>3) then
   runcl (3,7,8,'2pt','','ampl');end
 end  

$    Re-weight for final image if requested
$    and do wide-field if requested
      mname=target;seq 0;docal 1;gainuse 0;class 'MULTTB'
      source(1)=TARGET; setuvd
isline (line, fbchan, fechan, gbnd)
k=2;jj=-1
if (DOCHANS >0) then;k=0;end
if (DOPLZN >0) then; jj=1;end
    runspl (TARGET,'split',1,gbnd,k,jj,mbchan,mechan)
    if (dorewt > 0) then
     inn TARGET;incl 'split';inse 1;ind uvdisc           
     outn=innam; outcl=incl; outs=inseq; outd=ind; aparm 0
     antwt=newwt; go wtmod
    end 
$ *** NEED TO FIX DOTV
      mname=target;seq 0;docal -1;gainuse 0;class 'SPLIT'
      source(1)=TARGET;dowobbl=-1; mcells=acell;
      mimsize=timsize;doband -1;mniter 5000
            if (doplzn <=0) then; 
                runimg(loop,mflux,mniter,dowobbl,mbmaj,mbmin,mbpa); 
            end
            if (doplzn > 0) then; polimg(mbmaj,mbmin,mbpa);end
$
$    Calculate image RMS
$
      mname=TARGET; class='ICL001'; seq=0; uvdisc=mdisc; setuvd
            inse 0
$getrms(1,rms);
      mblc=0; mtrc=0; 
      plname = TARGET!!class!!char(seq)      
      pltmap(0,mblc,mblc,0); pkval=pixval;
$      inver=1; outver=1; go ccmrg
      if (doplzn > 0) then
         polcom(TARGET); 
      mname=TARGET; class='ICL001'; seq=0; uvdisc=mdisc; setuvd
         plname = TARGET!!class!!char(seq)!!'_Pol'
         pltmap(1,mblc,mblc,0)
         dozap (mdisc, TARGET, 'QCL001', 0, 'MA')
         dozap (mdisc, TARGET, 'UCL001', 0, 'MA')
      end
end
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 $ plot data fit & uv-coverage
 $
    if (doplot > 0) then
       task 'vplot'; uvdi=mdisc;mname=TARGET;class='MULTSC';seq=0
       setuvd; tbfind ('PL', ii); nplots=ntabs
       caldef;IFMULT; bchan=scb; echan=sce
       in2d=ind;in2name=innam;in2cl='ICL001'; in2seq=0; ncomp=-mniter,0
       nmaps=1; xinc 1; aparm 0; bparm 0; bparm(2)=-2;dotv =-1
       nplots (ntel-1); doarray 0; optyp ''; anten mrefant,0
       if (pkval > 0.3) then; solint=0; end
       if ((pkval <= 0.3) & (pkval > 0.05)) then; solint = 1; end
       if (pkval < 0.05) then; solint = 3; end; go;anten 0
       task 'uvplt'; docal -1; dopol -1; doband -1; xinc 1; factor 0;
       bparm 0; antennas 0; baseline 0; bparm 6 7 0; bchan=scb; 
       echan=bchan; go; tbfind ('PL',jj); nplots=ntabs;dotv =-1
        if (jj > ii) then; 
        plname = target!!'_uv'
          rplot (ii+1, jj); doextd ('PL', ii+1, jj)
        end
     end
       line = 0
 $
 $--------------------------------------------------------------------
 $
    if (line=1) then
       mess ('Generate and plot line cube')
       mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
       setuvd; doextd ('SN', 1, 5)
 $      exist (mdisc, TARGET, '1CAL', 1, texist)
 $      if (texist=true) then
 $         mname=TARGET; class='1CAL'; seq=mseq; uvd=mdisc; setuvd
 $         tbfind ('SN', ii)
 $         coptb (TARGET,'MULTTB',mdisc,mseq,'SN',ii,1)
 $         mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
 $         doextd ('CL', 5, 6)
$ *** CL Nos much be incremented by 1
$          runcl (1, 4, 5, '2PT', '', 'BOTH')
$       end
 $
       clrtemp
$       exist (mdisc, TARGET, '2CAL', 1, texist)
$       if (texist=true) then
$          mname=TARGET; class='2CAL'; seq=mseq; uvd=mdisc; setuvd
$          coptb (TARGET,'MULTTB',mdisc,mseq,'SN',1,2)
$          mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
$          runcl (2, 5, 6, '2PT', '', 'BOTH')
$       end
 $
$** This may need mult and source(1)=TARGET
$ ** set the velocity etc?  POSSM?
$       dozap (mdisc, TARGET, 'FCAL', 0, 'UV')
       mname=TARGET; class='MULTTB'; uvdisc=mdisc; seq=mseq; 
       v1=-1; if (doplzn>0) then; v1=1; end
       gbnd = 1; if (mdoband=0) then; gbnd = -1; end
       runspl (TARGET,'FCAL',1,gbnd,0,v1,1,0)
       mname=TARGET; class='FCAL'; uvdisc=mdisc; seq=mseq; 
       loop=1; mflux=0; mniter=1000; mimsize=timsize; dowobbl=-1;
       mcells=acell; mnbox=0; mclbox=0
       mbchan 1;mechan 0
       doband -1
       runimg (loop,mflux,mniter,dowobbl,0,0,0)
 $
 $  plot the cube
 $
 $ ** customise this for No chans up to 32 plots of 5x6 for 1023
isline(line,fbchan,fechan,gbnd)
       mname=TARGET; class='ICL001'; uvdisc=mdisc; seq=0; setuvd
       intype 'MA'; inseq 0; getseq; tbfind('PL', ii); nplots=ntabs
       v1 = ceil(fechan/12);task 'kntr';dogrey -1;dovect -1;zinc 0;
       clr2n;in2d ind;clr3n;in3d ind;clr4n;in4d ind;dotv -1
       xyratio 1; pixr 0; ltype 3; doalign -1; plev 0; doblank 1;
       dowedge -1; docircle -1; inver 0; stfactor 0; cbplot 12;
       altsw; label 1; dotv -1
       for kk = 1 to v1
          bchan = ((kk-1)*12) + 1; echan = bchan + 11
          if (echan > fechan) then; echan = fechan; end
          j = (bchan+echan)/2; getrms(j,rms)
          blc = 32, 32, bchan; trc = 97, 97, echan
          clev = 3 * rms; 
          levs -4,-2,-1,1,2,4,8,16,32,64,128,256,512,1024; 
          go kntr
       end
       tbfind ('PL', jj); nplots=ntabs
       if (jj > ii) then; 
        plname = target!!'_cub'
          rplot (ii+1, jj); doextd ('PL', ii+1, jj)
       end
    end
 $   version = orver
 $
 $------------------------------------------------------------
 recat
 return; finish
 tget addi
 vers inn
 inn '';incl '';inty '';inse 0