visname='J2157-694.ms'

# target J2157-694
# time on source 90sec
# N antennas 25
# Band 3 ~ 90GHz
# thermal noise 0.12 mJy
# peak of the target ~1 Jy

os.system("rm -rf noselfcal_PKS*")
clean(vis = visname,
  imagename = 'noselfcal_PKS',
  mode = 'mfs',
  interactive = T,
  imsize = [1024, 1024],  
  cell = '0.2arcsec',
  weighting = 'natural',
  niter = 10000,
  threshold='0.5mJy')

# clean con una box generica che includa tutta l'emissione che compare dopo aver
# cleanato il picco
# Ho misurato un rms dell'ordine di 0.8 mJy.


# Clean iniziale con weighting uniform, per identificare per bene 
# le sorgenti puntiformi

os.system("rm -rf firstclean_PKS*")
clean(vis = visname,
  imagename = 'firstclean_PKS',
  mode = 'mfs',
  interactive = T,
  imsize = [1024, 1024],  
  cell = '0.2arcsec',
  weighting = 'uniform',
  niter = 10000,
  threshold='0.5mJy')


### Fermarsi dopo aver cleanato le prime tre sorgenti ptformi brillanti 
### che appaiono una alla volta!!!

### Cercare soluzioni in fase in un intervallo di soluzioni di 15 sec (meta' scan)

os.system("rm -rf p_halfscan.cal")
gaincal(vis=visname,
        caltable="p_halfscan.cal",
        solint='15s',
        calmode="p",
        refant="DV04",
        minblperant=4,
        gaintype="G")

plotcal(caltable="p_halfscan.cal",
        xaxis="time",
        yaxis="phase",
        subplot=331,
        iteration="antenna",
        timerange='',
        spw='',
        plotrange=[0,0,-50,50],
        markersize=5,
        fontsize=10.0)

# L'ordine di grandezza delle soluzioni e' di ~10 deg (20 deg per qualche antenna)
# Applicare le soluzioni 

applycal(vis=visname,
         gaintable=["p_halfscan.cal"],
         interp="linear")

# Creare un nuovo modello

os.system("rm -rf selfcal_phalfscan_PKS*")
clean(vis = visname,
  imagename = 'selfcal_phalfscan_PKS',
  mode = 'mfs',
  interactive = T,
  imsize = [1024, 1024],  
  cell = '0.2arcsec',
  weighting = 'natural',
  niter = 10000,
  threshold='0.5mJy')

### Fermarsi dopo aver cleanato le due sorgenti ptformi brillanti
### e aggiungendo qualche componente brillante di cui ci si fida abbastanza

# Ricalcolare soluzioni in fase riducendo il tempo di soluzione

os.system("rm -rf p_int.cal")
gaincal(vis=visname,
        caltable="p_int.cal",
        solint='int',
        calmode="p",
        refant="DV04",
        minblperant=4,
        gaintype="G")

plotcal(caltable="p_int.cal",
        xaxis="time",
        yaxis="phase",
        subplot=331,
        iteration="antenna",
        timerange='',
        spw='',
        plotrange=[0,0,-20,20],
        markersize=5,
        fontsize=10.0)

# di nuovo lo stesso ordine di grandezza nelle correzioni
# applicare le soluzioni

applycal(vis=visname,
         gaintable=["p_int.cal"],
         interp="linear")

# Rifare immagine

os.system("rm -rf selfcal_p_int_PKS*")
clean(vis = visname,
  imagename = 'selfcal_p_int_PKS',
  mode = 'mfs',
  interactive = T,
  imsize = [1024, 1024],  
  cell = '0.2arcsec',
  weighting = 'natural',
  robust = 0.5,
  niter = 10000,
  threshold='0.5mJy')


# Calcolare ultime correzioni in fase e ampiezza
# in questo caso e' importante applicare le soluzioni in fase migliori che si sono ottenute

os.system("rm -rf ap_int.cal")
gaincal(vis=visname,
        caltable="ap_int.cal",
        solint='int',
        calmode="ap",
        gaintable ='p_int.cal',
        refant="DV04",
        minblperant=4,
        gaintype="G")


plotcal(caltable="ap_int.cal",
        xaxis="time",
        yaxis="phase",
        subplot=331,
        iteration="antenna",
        timerange='',
        spw='',
        plotrange=[0,0,-180,180],
        markersize=5,
        fontsize=10.0)


# I plot delle soluzioni in fase questa volta mostrano le fasi molto vicine a zero.
# Perche' abbiamo applicato la correzione in fase calcolata nell'iterazione precedente

# E' necessario ora applicare entrambe le tabelle ap e p ai dati.

applycal(vis=visname,
         gaintable=["ap_int.cal","p_int.cal" ],
         interp="linear")

# Rifare l'immagine.
# Questa volta e' quella finale, possiamo cleanare piu' a lungo!

os.system("rm -rf selfcal_ap_int__PKS*")
clean(vis = visname,
  imagename = 'selfcal_ap_int__PKS',
  mode = 'mfs',
  interactive = T,
  imsize = [1024, 1024],  
  cell = '0.2arcsec',
  weighting = 'natural',
  niter = 10000,
  threshold='0.4mJy')

# Una rms dell'ordine di 0.2 mJy.