<?xml version="1.0"?>
<feed xmlns="http://www.w3.org/2005/Atom" xml:lang="en">
	<id>http://ovsa.njit.edu//wiki/index.php?action=history&amp;feed=atom&amp;title=Tohban_EOVSA_Flare_Imaging_Tutorial_A-Z</id>
	<title>Tohban EOVSA Flare Imaging Tutorial A-Z - Revision history</title>
	<link rel="self" type="application/atom+xml" href="http://ovsa.njit.edu//wiki/index.php?action=history&amp;feed=atom&amp;title=Tohban_EOVSA_Flare_Imaging_Tutorial_A-Z"/>
	<link rel="alternate" type="text/html" href="http://ovsa.njit.edu//wiki/index.php?title=Tohban_EOVSA_Flare_Imaging_Tutorial_A-Z&amp;action=history"/>
	<updated>2026-04-14T09:54:28Z</updated>
	<subtitle>Revision history for this page on the wiki</subtitle>
	<generator>MediaWiki 1.38.1</generator>
	<entry>
		<id>http://ovsa.njit.edu//wiki/index.php?title=Tohban_EOVSA_Flare_Imaging_Tutorial_A-Z&amp;diff=6033&amp;oldid=prev</id>
		<title>Sshaik: Created page with &quot;=== Connection details to pipeline server === ==== Step 1: Downloading raw data (IDB) on Pipeline machine==== On CASA of the Pipeline machine, which has the complete EOVSA SQL...&quot;</title>
		<link rel="alternate" type="text/html" href="http://ovsa.njit.edu//wiki/index.php?title=Tohban_EOVSA_Flare_Imaging_Tutorial_A-Z&amp;diff=6033&amp;oldid=prev"/>
		<updated>2022-02-10T21:30:41Z</updated>

		<summary type="html">&lt;p&gt;Created page with &amp;quot;=== Connection details to pipeline server === ==== Step 1: Downloading raw data (IDB) on Pipeline machine==== On CASA of the Pipeline machine, which has the complete EOVSA SQL...&amp;quot;&lt;/p&gt;
&lt;p&gt;&lt;b&gt;New page&lt;/b&gt;&lt;/p&gt;&lt;div&gt;=== Connection details to pipeline server ===&lt;br /&gt;
==== Step 1: Downloading raw data (IDB) on Pipeline machine====&lt;br /&gt;
On CASA of the Pipeline machine, which has the complete EOVSA SQL database.&lt;br /&gt;
&lt;br /&gt;
If you are not using mobaxterm, directly SSH to Pipeline through your Linux or Mac computer and start tcsh to run CASA.&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
from astropy.time import Time&lt;br /&gt;
import os&lt;br /&gt;
trange = Time(['2017-08-21 20:15:00', '2017-08-21 20:25:00'])                               ###Change accordingly###&lt;br /&gt;
#### (Optional) change output path, default current directory &amp;quot;./&amp;quot; #####&lt;br /&gt;
outpath = './msdata/'                                                                       ###Change accordingly###&lt;br /&gt;
if not os.path.exists(outpath):&lt;br /&gt;
    os.makedirs(outpath)&lt;br /&gt;
######################################################&lt;br /&gt;
msfiles = importeovsa(idbfiles=trange, ncpu=1, visprefix=outpath)&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
NOTE: The above step currently gives an error with importeovsa. This is because the data location was changed and the code doesn't know where to look for the IDB files. Sijie is looking into the problem. We will update the page when it is fixed. For now, please use the method below for all time ranges. &lt;br /&gt;
&lt;br /&gt;
OR (for a time range longer than 10 minutes)&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
from suncasa.tasks import task_calibeovsa as calibeovsa&lt;br /&gt;
from suncasa.tasks import task_importeovsa as timporteovsa&lt;br /&gt;
from split_cli import split_cli as split&lt;br /&gt;
import dump_tsys as dt&lt;br /&gt;
from util import Time&lt;br /&gt;
import numpy as np&lt;br /&gt;
import os&lt;br /&gt;
from glob import glob&lt;br /&gt;
from eovsapy import util&lt;br /&gt;
&lt;br /&gt;
trange = Time(['2017-08-21 20:15:00', '2017-08-21 20:35:00'])                               ###Change accordingly###&lt;br /&gt;
idbdir = util.get_idbdir(trange[0])&lt;br /&gt;
&lt;br /&gt;
info = dt.rd_fdb(trange[0])&lt;br /&gt;
for k, v in info.iteritems():&lt;br /&gt;
    info[k] = info[k][~(info[k] == '')]&lt;br /&gt;
&lt;br /&gt;
sidx = np.where(&lt;br /&gt;
    np.logical_and(info['SOURCEID'] == 'Sun', info['PROJECTID'] == 'NormalObserving') &amp;amp; np.logical_and(&lt;br /&gt;
        info['ST_TS'].astype(np.float) &amp;gt;= trange[0].lv,&lt;br /&gt;
        info['ST_TS'].astype(np.float) &amp;lt;= trange[&lt;br /&gt;
            1].lv))&lt;br /&gt;
filelist = info['FILE'][sidx]&lt;br /&gt;
&lt;br /&gt;
outpath = './msdata/'                                                                       ###Change accordingly###&lt;br /&gt;
if not os.path.exists(outpath):&lt;br /&gt;
    os.makedirs(outpath)&lt;br /&gt;
inpath = idbdir + '{}/'.format(trange[0].datetime.strftime(&amp;quot;%Y%m%d&amp;quot;))&lt;br /&gt;
ncpu = 1&lt;br /&gt;
&lt;br /&gt;
msfiles = timporteovsa.importeovsa(idbfiles=[inpath + ll for ll in filelist], ncpu=ncpu, timebin=&amp;quot;0s&amp;quot;, width=1,&lt;br /&gt;
                                   visprefix=outpath,&lt;br /&gt;
                                   nocreatms=False, doconcat=False,&lt;br /&gt;
                                   modelms=&amp;quot;&amp;quot;, doscaling=False, keep_nsclms=False, udb_corr=True)&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
If you come across errors with calibeovsa, add following lines to your ~/.casa/init.py file.&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
import sys&lt;br /&gt;
sys.path.append('/common/python')&lt;br /&gt;
sys.path.append('/common/python/packages/pipeline_casa')&lt;br /&gt;
execfile('/common/python/suncasa/tasks/mytasks.py')&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
==== Step 2: Concatenate all the 10 mins data, if there are any====&lt;br /&gt;
Follow this step if there are more than one .ms files, if not run step 3 directly. If doimage=True, a quicklook image will be produced (by integrating over the entire time) as shown below.&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
# This is to set the path/name for the concatenated files&lt;br /&gt;
concatvis = os.path.basename(msfiles[0])[:11] + '_concat_cal.ms'&lt;br /&gt;
vis = calibeovsa.calibeovsa(msfiles, doconcat=True, concatvis=concatvis, caltype=['refpha','phacal'], doimage=True)&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
[[file:Figure1_imagingtutorial.png|frame|center|800px|'''Figure 1:''' Quick-look full-Sun image after the initial calibration.]]&lt;br /&gt;
&lt;br /&gt;
==== Step 3: Calibration ====&lt;br /&gt;
This will calibrate the input visibility, write out calibration tables under /data1/eovsa/caltable/, and apply the calibration. &lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
vis = calibeovsa.calibeovsa('IDB20170821202020.ms', caltype=['refpha','phacal'], doimage=True)            ###Change the vis filename accordingly###&lt;br /&gt;
# Append '_cal' to the ms filename and split the corrected column to the new caled ms&lt;br /&gt;
vis_str = str(' '.join(vis))&lt;br /&gt;
caled_vis=vis_str.replace('.ms','_cal.ms')&lt;br /&gt;
split(vis=' '.join(vis),outputvis=caled_vis,datacolumn='corrected',timerange='',correlation='XX')&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
One needs to transfer the created caled ms data files (xxx_concat_cal.ms or xxx_cal.ms) from pipeline to inti server, which has all the casatasks installed, in order to run the rest of the imaging steps.&lt;br /&gt;
&lt;br /&gt;
=== Connection details to Inti server ===&lt;br /&gt;
For Windows, on Mobaxterm,&lt;br /&gt;
&lt;br /&gt;
#With your NJIT VPN connected, connect to one of the afsconnect servers (for example, afsaccess3.njit.edu) using your UCID and password.&lt;br /&gt;
#ssh -X UCID@inti.hpcnet.campus.njit.edu &lt;br /&gt;
&lt;br /&gt;
To avoid typing the full inti address each time you attempt for ssh, you may wish to add the following lines with your username to C:\Program Files\Git\etc\ssh\ssh_config on Windows and /.ssh/config on Mac and Linux machines.&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
Host inti&lt;br /&gt;
Hostname inti.hpcnet.campus.njit.edu&lt;br /&gt;
User USERNAME&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
To have sufficient disk space with EOVSA data analysis over Inti, use your dedicated directory YOURDIRECTORY at the location given below. If you do not have a directory, Please take help from Sijie in creating one. &lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
cd /inti/data/users/YOURDIRECTORY     &lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
For Linux or Mac machine,&lt;br /&gt;
ssh -X UCID@inti.hpcnet.campus.njit.edu&lt;br /&gt;
&lt;br /&gt;
=== Transferring data files between servers ===&lt;br /&gt;
For directly transferring your calibrated .ms data between the Pipeline and Inti servers, follow the below given steps. &lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
1. Log into Inti using your username.&lt;br /&gt;
ssh -X USERNAME@inti.hpcnet.campus.njit.edu&lt;br /&gt;
&lt;br /&gt;
Then create a tunnel into Pipeline from Inti.&lt;br /&gt;
ssh -L 8888:pipeline.solar.pvt:22 guest@ovsa.njit.edu&lt;br /&gt;
&lt;br /&gt;
2. Log into Inti again from a new terminal.&lt;br /&gt;
&lt;br /&gt;
Change to your working directory and give this command to copy your data on Pipeline.&lt;br /&gt;
scp -v -C -r -P 8888 USERNAMEofPipeline@localhost:PATHofMSDATA/MSfilename ./&lt;br /&gt;
Eg: scp -v -C -r -P 8888 shaheda@localhost:/data1/shaheda/IDB20220118173922_cal.ms ./&lt;br /&gt;
where, MSfilename, PATHofMSDATA are your .ms data and its path on Pipeline machine.&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Alternatively, one can follow the below procedure to do the transfer.&lt;br /&gt;
&lt;br /&gt;
For Windows, on mobaxterm, drag and drop the ms file to your local machine or use scp command as given below.&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
scp -r -C userid@pipeline:/your/folder/msfile /localmachine/destination/folder/&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
On mobaxterm and on your local terminal, use the following command to finally copy the ms file to Inti.&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
scp -r -C /localmachine/destination/folder/msfile UCID@inti.hpcnet.campus.njit.edu:/inti/data/users/YOURDIRECTORY&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
For Mac and Linux, SCP can be used in the same way. Add the following lines to your SSH config file, to bypass ovsa.njit.edu from copying.&lt;br /&gt;
vi ~/.ssh/config&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
Host ovsa&lt;br /&gt;
        HostName ovsa.njit.edu&lt;br /&gt;
        User guest&lt;br /&gt;
Host pipeline&lt;br /&gt;
        Hostname pipeline.solar.pvt&lt;br /&gt;
        User userid&lt;br /&gt;
        ProxyCommand ssh -W %h:%p guest@ovsa.njit.edu&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
=== Software details on the servers ===&lt;br /&gt;
On Inti,&lt;br /&gt;
when logging in for the first time, please add the following lines to your accounts ~/.bashrc file.&lt;br /&gt;
&lt;br /&gt;
&amp;gt;&amp;gt;vi ~/.bashrc  &lt;br /&gt;
Insert the text given below and save it.&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
#### setting start ####&lt;br /&gt;
if [ $HOSTNAME == &amp;quot;baozi.hpcnet.campus.njit.edu&amp;quot; ]; then&lt;br /&gt;
    source /srg/.setenv_baozi&lt;br /&gt;
fi&lt;br /&gt;
if [ $HOSTNAME == &amp;quot;inti.hpcnet.campus.njit.edu&amp;quot; ]; then&lt;br /&gt;
    source /inti/.setenv_inti&lt;br /&gt;
fi&lt;br /&gt;
if [ $HOSTNAME == &amp;quot;guko.resource.campus.njit.edu&amp;quot; ]; then&lt;br /&gt;
    source /data/data/.setenv_guko&lt;br /&gt;
fi&lt;br /&gt;
#### setting end ####&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Both CASA 5 and 6 are available on Inti. &lt;br /&gt;
&lt;br /&gt;
Enter the bash environment on inti, and load the desired casa environment.&lt;br /&gt;
To load CASA 5, enter bash environment by giving &amp;gt;&amp;gt;bash &lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
&amp;gt;&amp;gt; loadcasa5&lt;br /&gt;
&amp;gt;&amp;gt; suncasa              #This should load the software making you ready for the analysis&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
Or to load CASA 6, enter bash environment by giving &amp;gt;&amp;gt;bash &lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
&amp;gt;&amp;gt; loadcasa6&lt;br /&gt;
&amp;gt;&amp;gt; ipython              #This should load the software making you ready for the analysis&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
Here, for example, to use clean, first start ipython as given above, then type in &amp;gt;&amp;gt;from casatasks import tclean&lt;br /&gt;
&lt;br /&gt;
====Step 4: Self-calibration on Inti server ====&lt;br /&gt;
https://github.com/binchensun/casa-eovsa/blob/master/slfcal_example.py&lt;br /&gt;
[[file:Figure2_imagingtutorial.png|thumb|right|500px|Figure 2: Full-Sun image after initial clean to find the flare location.]]&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
from suncasa.utils import helioimage2fits as hf&lt;br /&gt;
import os&lt;br /&gt;
import numpy as np&lt;br /&gt;
import pickle&lt;br /&gt;
from matplotlib import gridspec as gridspec&lt;br /&gt;
from sunpy import map as smap&lt;br /&gt;
from matplotlib import pyplot as plt&lt;br /&gt;
from split_cli import split_cli as split&lt;br /&gt;
import time&lt;br /&gt;
&lt;br /&gt;
# =========== task handlers =============&lt;br /&gt;
dofullsun = 1 # initial full-sun imaging                                         ###Change accordingly###&lt;br /&gt;
domasks=1 # get masks                                                            ###Change accordingly###&lt;br /&gt;
doslfcal=1 # main cycle of doing selfcalibration                                 ###Change accordingly###&lt;br /&gt;
doapply=1 # apply the results                                                    ###Change accordingly###&lt;br /&gt;
doclean_slfcaled=1 # perform clean for self-calibrated data                      ###Change accordingly###&lt;br /&gt;
&lt;br /&gt;
# ============ declaring the working directories ============&lt;br /&gt;
workdir = os.getcwd()+'/' #main working directory. Using current directory in this example&lt;br /&gt;
slfcaldir = workdir+'slfcal/' #place to put all selfcalibration products&lt;br /&gt;
imagedir = slfcaldir+'images/' #place to put all selfcalibration images&lt;br /&gt;
maskdir = slfcaldir+'masks/' #place to put clean masks&lt;br /&gt;
imagedir_slfcaled = slfcaldir+'images_slfcaled/' #place to put final self-calibrated images&lt;br /&gt;
caltbdir = slfcaldir+'caltbs/' # place to put calibration tables&lt;br /&gt;
# make these directories if they do not already exist&lt;br /&gt;
dirs = [workdir, slfcaldir, imagedir, maskdir, imagedir_slfcaled, caltbdir]&lt;br /&gt;
for d in dirs:&lt;br /&gt;
    if not os.path.exists(d):&lt;br /&gt;
        os.makedirs(d)&lt;br /&gt;
&lt;br /&gt;
# ============ Split a short time for self-calibration ===========&lt;br /&gt;
# input visibility&lt;br /&gt;
ms_in = workdir + 'IDB20170821202020_cal.ms'                                    ###Change the initial calibrated (through calibeovsa) vis accordingly###&lt;br /&gt;
# output, selfcaled, visibility&lt;br /&gt;
ms_slfcaled = workdir + os.path.basename(ms_in).replace('cal','slfcaled') &lt;br /&gt;
# intermediate small visibility for selfcalbration &lt;br /&gt;
# selected time range for generating self-calibration solutions&lt;br /&gt;
trange='2017/08/21/20:21:10~2017/08/21/20:21:30'                                ###Change accordingly###&lt;br /&gt;
slfcalms = slfcaldir+'slfcalms.XX.slfcal'&lt;br /&gt;
slfcaledms = slfcaldir+'slfcalms.XX.slfcaled'&lt;br /&gt;
if not os.path.exists(slfcalms):&lt;br /&gt;
    split(vis=ms_in,outputvis=slfcalms,datacolumn='data',timerange=trange,correlation='XX')&lt;br /&gt;
&lt;br /&gt;
# ============ Prior definitions for spectral windows, antennas, pixel numbers =========&lt;br /&gt;
spws=[str(s+1) for s in range(30)]&lt;br /&gt;
antennas='0~12' &lt;br /&gt;
npix=512&lt;br /&gt;
nround=3 #number of slfcal cycles&lt;br /&gt;
&lt;br /&gt;
# =========== Step 1, doing a full-Sun image to find out phasecenter and appropriate field of view =========&lt;br /&gt;
if dofullsun:&lt;br /&gt;
    #initial mfs clean to find out the image phase center&lt;br /&gt;
    im_init='fullsun_init'&lt;br /&gt;
    os.system('rm -rf '+im_init+'*')&lt;br /&gt;
    tclean(vis=slfcalms,&lt;br /&gt;
            antenna='0~12',&lt;br /&gt;
            imagename=im_init,&lt;br /&gt;
            spw='1~15',&lt;br /&gt;
            specmode='mfs',&lt;br /&gt;
            timerange=trange,&lt;br /&gt;
            imsize=[npix],&lt;br /&gt;
            cell=['5arcsec'],&lt;br /&gt;
            niter=1000,&lt;br /&gt;
            gain=0.05,&lt;br /&gt;
            stokes='I',&lt;br /&gt;
            restoringbeam=['30arcsec'],&lt;br /&gt;
            interactive=False,&lt;br /&gt;
            pbcor=True)&lt;br /&gt;
&lt;br /&gt;
    hf.imreg(vis=slfcalms,imagefile=im_init+'.image.pbcor',fitsfile=im_init+'.fits',&lt;br /&gt;
             timerange=trange,usephacenter=False,verbose=True)&lt;br /&gt;
    clnjunks = ['.flux', '.mask', '.model', '.psf', '.residual','.sumwt','.pb','.image']&lt;br /&gt;
    for clnjunk in clnjunks:&lt;br /&gt;
        if os.path.exists(im_init + clnjunk):&lt;br /&gt;
            os.system('rm -rf '+im_init + clnjunk)&lt;br /&gt;
&lt;br /&gt;
    from sunpy import map as smap&lt;br /&gt;
    from matplotlib import pyplot as plt&lt;br /&gt;
    fig = plt.figure(figsize=(6,6))&lt;br /&gt;
    ax = fig.add_subplot(111)&lt;br /&gt;
    eomap=smap.Map(im_init+'.fits')&lt;br /&gt;
    #eomap.data=eomap.data.reshape((npix,npix))&lt;br /&gt;
    eomap.plot_settings['cmap'] = plt.get_cmap('jet')&lt;br /&gt;
    eomap.plot(axes = ax)&lt;br /&gt;
    eomap.draw_limb()&lt;br /&gt;
    plt.show()&lt;br /&gt;
    viewer(im_init+'.image.pbcor')&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
# parameters specific to the event (found from step 1)&lt;br /&gt;
phasecenter='J2000 10h02m59 11d58m07'                                               ###Change accordingly###&lt;br /&gt;
xran=[280,480]                                                                      ###Change accordingly###&lt;br /&gt;
yran=[-50,150]                                                                      ###Change accordingly###&lt;br /&gt;
&lt;br /&gt;
# =========== Step 2 (optional), generate masks =========&lt;br /&gt;
# if skipped, will not use any masks&lt;br /&gt;
if domasks:&lt;br /&gt;
    clearcal(slfcalms)&lt;br /&gt;
    delmod(slfcalms)&lt;br /&gt;
    antennas='0~12' &lt;br /&gt;
    pol='XX'&lt;br /&gt;
    imgprefix=maskdir+'slf_t0'&lt;br /&gt;
&lt;br /&gt;
    # step 1: set up the clean masks&lt;br /&gt;
    img_init=imgprefix+'_init_ar_'&lt;br /&gt;
    os.system('rm -rf '+img_init+'*')&lt;br /&gt;
    #spwrans_mask=['1~5','6~12','13~20','21~30']&lt;br /&gt;
    spwrans_mask=['1~12']&lt;br /&gt;
    #convert to a list of spws&lt;br /&gt;
    spwrans_mask_list = [[str(i) for i in (np.arange(int(m.split('~')[0]),int(m.split('~')[1])))] for m in spwrans_mask]&lt;br /&gt;
    masks=[]&lt;br /&gt;
    imnames=[]&lt;br /&gt;
    for spwran in spwrans_mask:&lt;br /&gt;
        imname=img_init+spwran.replace('~','-')&lt;br /&gt;
        try:&lt;br /&gt;
            tclean(vis=slfcalms,&lt;br /&gt;
                    antenna='0~12',&lt;br /&gt;
                    imagename=imname,&lt;br /&gt;
                    spw=spwran,&lt;br /&gt;
                    specmode='mfs',&lt;br /&gt;
                    timerange=trange,&lt;br /&gt;
                    imsize=[npix],&lt;br /&gt;
                    cell=['2arcsec'],&lt;br /&gt;
                    niter=1000,&lt;br /&gt;
                    gain=0.05,&lt;br /&gt;
                    stokes='XX',&lt;br /&gt;
                    restoringbeam=['20arcsec'],&lt;br /&gt;
                    phasecenter=phasecenter,&lt;br /&gt;
                    weighting='briggs',&lt;br /&gt;
                    robust=1.0,&lt;br /&gt;
                    interactive=True,&lt;br /&gt;
		    datacolumn='data',&lt;br /&gt;
                    pbcor=True,&lt;br /&gt;
                    savemodel='modelcolumn')&lt;br /&gt;
            imnames.append(imname+'.image')&lt;br /&gt;
            masks.append(imname+'.mask')&lt;br /&gt;
            clnjunks = ['.flux', '.model', '.psf', '.residual']&lt;br /&gt;
            for clnjunk in clnjunks:&lt;br /&gt;
                if os.path.exists(imname + clnjunk):&lt;br /&gt;
                    os.system('rm -rf '+ imname + clnjunk)&lt;br /&gt;
        except:&lt;br /&gt;
            print('error in cleaning spw: '+spwran)&lt;br /&gt;
&lt;br /&gt;
    pickle.dump(masks,open(slfcaldir+'masks.p','wb'))&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
# =========== Step 3, main step of selfcalibration =========&lt;br /&gt;
if doslfcal:&lt;br /&gt;
    if os.path.exists(slfcaldir+'masks.p'):&lt;br /&gt;
        masks=pickle.load(open(slfcaldir+'masks.p','rb'))&lt;br /&gt;
    if not os.path.exists(slfcaldir+'masks.p'):&lt;br /&gt;
        print 'masks do not exist. Use default mask'&lt;br /&gt;
        masks=[]&lt;br /&gt;
    os.system('rm -rf '+imagedir+'*')&lt;br /&gt;
    os.system('rm -rf '+caltbdir+'*')&lt;br /&gt;
    #first step: make a mock caltable for the entire database&lt;br /&gt;
    print('Processing ' + trange)&lt;br /&gt;
    slftbs=[]&lt;br /&gt;
    calprefix=caltbdir+'slf'&lt;br /&gt;
    imgprefix=imagedir+'slf'&lt;br /&gt;
    tb.open(slfcalms+'/SPECTRAL_WINDOW')&lt;br /&gt;
    reffreqs=tb.getcol('REF_FREQUENCY')&lt;br /&gt;
    bdwds=tb.getcol('TOTAL_BANDWIDTH')&lt;br /&gt;
    cfreqs=reffreqs+bdwds/2.&lt;br /&gt;
    tb.close()&lt;br /&gt;
    # starting beam size at 3.4 GHz in arcsec&lt;br /&gt;
    sbeam=40.&lt;br /&gt;
    strtmp=[m.replace(':','') for m in trange.split('~')]&lt;br /&gt;
    timestr='t'+strtmp[0]+'-'+strtmp[1]&lt;br /&gt;
    refantenna='0'&lt;br /&gt;
    # number of iterations for each round&lt;br /&gt;
    niters=[100, 300, 500]&lt;br /&gt;
    # roburst value for weighting the baselines&lt;br /&gt;
    robusts=[1.0, 0.5, 0.0]&lt;br /&gt;
    # apply calibration tables? Set to true for most cases&lt;br /&gt;
    doapplycal=[1,1,1]&lt;br /&gt;
    # modes for calibration, 'p' for phase-only, 'a' for amplitude only, 'ap' for both&lt;br /&gt;
    calmodes=['p','p','a']&lt;br /&gt;
    # setting uvranges for model image (optional, not used here)&lt;br /&gt;
    uvranges=['','',''] &lt;br /&gt;
    for n in range(nround):&lt;br /&gt;
        slfcal_tb_g= calprefix+'.G'+str(n)&lt;br /&gt;
        fig = plt.figure(figsize=(8.4,7.))&lt;br /&gt;
        gs = gridspec.GridSpec(5, 6)&lt;br /&gt;
        for s,sp in enumerate(spws):&lt;br /&gt;
            print 'processing spw: '+sp&lt;br /&gt;
            cfreq=cfreqs[int(sp)]&lt;br /&gt;
            # setting restoring beam size (not very useful for selfcal anyway, but just to see the results)&lt;br /&gt;
            bm=max(sbeam*cfreqs[1]/cfreq, 6.)&lt;br /&gt;
            slfcal_img = imgprefix+'.spw'+sp.zfill(2)+'.slfcal'+str(n)&lt;br /&gt;
            # only the first round uses nearby spws for getting initial model&lt;br /&gt;
            if n == 0:&lt;br /&gt;
                spbg=max(int(sp)-2,1)&lt;br /&gt;
                sped=min(int(sp)+2,30)&lt;br /&gt;
                spwran=str(spbg)+'~'+str(sped)&lt;br /&gt;
                print('using spw {0:s} as model'.format(spwran))&lt;br /&gt;
                if 'spwrans_mask' in vars():&lt;br /&gt;
                    for m, spwran_mask in enumerate(spwrans_mask):&lt;br /&gt;
                        if sp in spwran_mask:&lt;br /&gt;
                            mask = masks[m]&lt;br /&gt;
                            print('using mask {0:s}'.format(mask))&lt;br /&gt;
                            findmask = True&lt;br /&gt;
                    if not findmask:&lt;br /&gt;
                        print('mask not found. Do use any masks')&lt;br /&gt;
            else:&lt;br /&gt;
                spwran = sp&lt;br /&gt;
                if 'spwrans_mask' in vars():&lt;br /&gt;
                    for m, spwran_mask in enumerate(spwrans_mask):&lt;br /&gt;
                        if sp in spwran_mask:&lt;br /&gt;
                            mask = masks[m]&lt;br /&gt;
                            print 'using mask {0:s}'.format(mask)&lt;br /&gt;
                            findmask = True&lt;br /&gt;
                    if not findmask:&lt;br /&gt;
                        print('mask not found. Do use any masks')&lt;br /&gt;
            try:&lt;br /&gt;
                tclean(vis=slfcalms,&lt;br /&gt;
                        antenna=antennas,&lt;br /&gt;
                        imagename=slfcal_img,&lt;br /&gt;
                        uvrange=uvranges[n],&lt;br /&gt;
                        spw=spwran,&lt;br /&gt;
                        specmode='mfs',&lt;br /&gt;
                        timerange=trange,&lt;br /&gt;
                        imsize=[npix],&lt;br /&gt;
                        cell=['2arcsec'],&lt;br /&gt;
                        niter=niters[n],&lt;br /&gt;
                        gain=0.05,&lt;br /&gt;
                        stokes='XX', #use pol XX image as the model&lt;br /&gt;
                        weighting='briggs',&lt;br /&gt;
                        robust=robusts[n],&lt;br /&gt;
                        phasecenter=phasecenter,&lt;br /&gt;
                        mask=mask,&lt;br /&gt;
                        restoringbeam=[str(bm)+'arcsec'],&lt;br /&gt;
                        pbcor=False,&lt;br /&gt;
                        interactive=False,&lt;br /&gt;
                        savemodel='modelcolumn')&lt;br /&gt;
                if os.path.exists(slfcal_img+'.image'):&lt;br /&gt;
                    fitsfile=slfcal_img+'.fits'&lt;br /&gt;
                    hf.imreg(vis=slfcalms,imagefile=slfcal_img+'.image',fitsfile=fitsfile,&lt;br /&gt;
                             timerange=trange,usephacenter=False,toTb=True,verbose=False,overwrite=True)&lt;br /&gt;
                clnjunks = ['.mask','.flux', '.model', '.psf', '.residual', '.image','.pb','.image.pbcor','.sumwt']&lt;br /&gt;
                for clnjunk in clnjunks:&lt;br /&gt;
                    if os.path.exists(slfcal_img + clnjunk):&lt;br /&gt;
                        os.system('rm -rf '+ slfcal_img + clnjunk)&lt;br /&gt;
                ax = fig.add_subplot(gs[s])&lt;br /&gt;
                eomap=smap.Map(fitsfile)&lt;br /&gt;
                eomap.plot_settings['cmap'] = plt.get_cmap('jet')&lt;br /&gt;
                eomap.plot(axes = ax)&lt;br /&gt;
                eomap.draw_limb()&lt;br /&gt;
                #eomap.draw_grid()&lt;br /&gt;
                ax.set_title(' ')&lt;br /&gt;
		ax.get_xaxis().set_visible(False)&lt;br /&gt;
		ax.get_yaxis().set_visible(False)&lt;br /&gt;
                ax.set_xlim(xran)&lt;br /&gt;
                ax.set_ylim(yran)&lt;br /&gt;
                os.system('rm -f '+ fitsfile)&lt;br /&gt;
&lt;br /&gt;
            except:&lt;br /&gt;
                print 'error in cleaning spw: '+sp&lt;br /&gt;
                print 'using nearby spws for initial model'&lt;br /&gt;
                sp_e=int(sp)+2&lt;br /&gt;
                sp_i=int(sp)-2&lt;br /&gt;
                if sp_i &amp;lt; 1:&lt;br /&gt;
                    sp_i = 1&lt;br /&gt;
                if sp_e &amp;gt; 30:&lt;br /&gt;
                    sp_e = 30&lt;br /&gt;
                sp_=str(sp_i)+'~'+str(sp_e)&lt;br /&gt;
                try:&lt;br /&gt;
                    tget(tclean)&lt;br /&gt;
                    spw=sp_&lt;br /&gt;
                    print('using spw {0:s} as model'.format(sp_))&lt;br /&gt;
                    tclean()&lt;br /&gt;
                except:&lt;br /&gt;
                    print 'still not successful. abort...'&lt;br /&gt;
                    break&lt;br /&gt;
&lt;br /&gt;
            gaincal(vis=slfcalms, refant=refantenna,antenna=antennas,caltable=slfcal_tb_g,spw=sp, uvrange='',\&lt;br /&gt;
                    gaintable=[],selectdata=True,timerange=trange,solint='inf',gaintype='G',calmode=calmodes[n],\&lt;br /&gt;
                    combine='',minblperant=4,minsnr=2,append=True)&lt;br /&gt;
            if not os.path.exists(slfcal_tb_g):&lt;br /&gt;
                print 'No solution found in spw: '+sp&lt;br /&gt;
        figname=imagedir+'slf_t0_n{:d}.png'.format(n)&lt;br /&gt;
	plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)&lt;br /&gt;
        plt.savefig(figname)&lt;br /&gt;
        time.sleep(10)&lt;br /&gt;
        plt.close()&lt;br /&gt;
&lt;br /&gt;
        if os.path.exists(slfcal_tb_g):&lt;br /&gt;
            slftbs.append(slfcal_tb_g)&lt;br /&gt;
            slftb=[slfcal_tb_g]&lt;br /&gt;
            os.chdir(slfcaldir)&lt;br /&gt;
            if calmodes[n] == 'p': &lt;br /&gt;
                plotcal(caltable=slfcal_tb_g,antenna='1~12',xaxis='freq',yaxis='phase',\&lt;br /&gt;
                        subplot=431,plotrange=[-1,-1,-180,180],iteration='antenna',figfile=slfcal_tb_g+'.png',showgui=False)&lt;br /&gt;
            if calmodes[n] == 'a':&lt;br /&gt;
                plotcal(caltable=slfcal_tb_g,antenna='1~12',xaxis='freq',yaxis='amp',\&lt;br /&gt;
                        subplot=431,plotrange=[-1,-1,0,2.],iteration='antenna',figfile=slfcal_tb_g+'.png',showgui=False)&lt;br /&gt;
            os.chdir(workdir)&lt;br /&gt;
&lt;br /&gt;
&lt;br /&gt;
        if doapplycal[n]:&lt;br /&gt;
            clearcal(slfcalms)&lt;br /&gt;
            delmod(slfcalms)&lt;br /&gt;
            applycal(vis=slfcalms,gaintable=slftb,spw=','.join(spws),selectdata=True,\&lt;br /&gt;
                     antenna=antennas,interp='nearest',flagbackup=False,applymode='calonly',calwt=False)&lt;br /&gt;
&lt;br /&gt;
        if n &amp;lt; nround-1: &lt;br /&gt;
            prompt=raw_input('Continuing to selfcal?')&lt;br /&gt;
            #prompt='y'&lt;br /&gt;
            if prompt.lower() == 'n':&lt;br /&gt;
                if os.path.exists(slfcaledms):&lt;br /&gt;
                    os.system('rm -rf '+slfcaledms)&lt;br /&gt;
                split(slfcalms,slfcaledms,datacolumn='corrected')&lt;br /&gt;
                print 'Final calibrated ms is {0:s}'.format(slfcaledms)&lt;br /&gt;
                break&lt;br /&gt;
            if prompt.lower() == 'y':&lt;br /&gt;
                slfcalms_=slfcalms+str(n)&lt;br /&gt;
                if os.path.exists(slfcalms_):&lt;br /&gt;
                    os.system('rm -rf '+slfcalms_)&lt;br /&gt;
                split(slfcalms,slfcalms_,datacolumn='corrected')&lt;br /&gt;
                slfcalms=slfcalms_&lt;br /&gt;
        else:&lt;br /&gt;
            if os.path.exists(slfcaledms):&lt;br /&gt;
                os.system('rm -rf '+slfcaledms)&lt;br /&gt;
            split(slfcalms,slfcaledms,datacolumn='corrected')&lt;br /&gt;
            print 'Final calibrated ms is {0:s}'.format(slfcaledms)&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
# =========== Step 4: Apply self-calibration tables =========&lt;br /&gt;
if doapply:&lt;br /&gt;
    import glob&lt;br /&gt;
    os.chdir(workdir)&lt;br /&gt;
    clearcal(ms_in)&lt;br /&gt;
    clearcal(slfcalms)&lt;br /&gt;
    applycal(vis=slfcalms,gaintable=slftbs,spw=','.join(spws),selectdata=True,\&lt;br /&gt;
             antenna=antennas,interp='linear',flagbackup=False,applymode='calonly',calwt=False)&lt;br /&gt;
    applycal(vis=ms_in,gaintable=slftbs,spw=','.join(spws),selectdata=True,\&lt;br /&gt;
             antenna=antennas,interp='linear',flagbackup=False,applymode='calonly',calwt=False)&lt;br /&gt;
    if os.path.exists(ms_slfcaled):&lt;br /&gt;
        os.system('rm -rf '+ms_slfcaled)&lt;br /&gt;
    split(ms_in, ms_slfcaled,datacolumn='corrected')&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
[[file:Slf.G0.png|thumb|center|500px|Figure 1: Phase before self-calibration]]&lt;br /&gt;
[[file:Slf.G1.png|thumb|center|500px|Figure 2: Phase after self-calibration]]&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
# =========== Step 5: Generate final self-calibrated images (optional) =========&lt;br /&gt;
if doclean_slfcaled:&lt;br /&gt;
    import glob&lt;br /&gt;
    pol='XX'&lt;br /&gt;
    print('Processing ' + trange)&lt;br /&gt;
    img_final=imagedir_slfcaled+'/slf_final_{0}_t0'.format(pol)&lt;br /&gt;
    vis = ms_slfcaled&lt;br /&gt;
    spws=[str(s+1) for s in range(30)]&lt;br /&gt;
    tb.open(vis+'/SPECTRAL_WINDOW')&lt;br /&gt;
    reffreqs=tb.getcol('REF_FREQUENCY')&lt;br /&gt;
    bdwds=tb.getcol('TOTAL_BANDWIDTH')&lt;br /&gt;
    cfreqs=reffreqs+bdwds/2.&lt;br /&gt;
    tb.close()&lt;br /&gt;
    sbeam=30.&lt;br /&gt;
    from matplotlib import gridspec as gridspec&lt;br /&gt;
    from sunpy import map as smap&lt;br /&gt;
    from matplotlib import pyplot as plt&lt;br /&gt;
    fitsfiles=[]&lt;br /&gt;
    for s,sp in enumerate(spws):&lt;br /&gt;
        cfreq=cfreqs[int(sp)]&lt;br /&gt;
        bm=max(sbeam*cfreqs[1]/cfreq,4.)&lt;br /&gt;
        imname=img_final+'_s'+sp.zfill(2)&lt;br /&gt;
        fitsfile=imname+'.fits'&lt;br /&gt;
        if not os.path.exists(fitsfile):&lt;br /&gt;
            print 'cleaning spw {0:s} with beam size {1:.1f}&amp;quot;'.format(sp,bm)&lt;br /&gt;
            try:&lt;br /&gt;
                tclean(vis=vis,&lt;br /&gt;
                        antenna=antennas,&lt;br /&gt;
                        imagename=imname,&lt;br /&gt;
                        spw=sp,&lt;br /&gt;
                        specmode='mfs',&lt;br /&gt;
                        timerange=trange,&lt;br /&gt;
                        imsize=[npix],&lt;br /&gt;
                        cell=['1arcsec'],&lt;br /&gt;
                        niter=1000,&lt;br /&gt;
                        gain=0.05,&lt;br /&gt;
                        stokes=pol,&lt;br /&gt;
                        weighting='briggs',&lt;br /&gt;
                        robust=2.0,&lt;br /&gt;
                        restoringbeam=[str(bm)+'arcsec'],&lt;br /&gt;
                        phasecenter=phasecenter,&lt;br /&gt;
                        mask='',&lt;br /&gt;
                        pbcor=True,&lt;br /&gt;
                        interactive=False)&lt;br /&gt;
            except:&lt;br /&gt;
                print 'cleaning spw '+sp+' unsuccessful. Proceed to next spw'&lt;br /&gt;
                continue&lt;br /&gt;
            if os.path.exists(imname+'.image.pbcor'):&lt;br /&gt;
                imn = imname+'.image.pbcor'&lt;br /&gt;
                hf.imreg(vis=vis,imagefile=imn,fitsfile=fitsfile,&lt;br /&gt;
                         timerange=trange,usephacenter=False,toTb=True,verbose=False)&lt;br /&gt;
            fitsfiles.append(fitsfile)&lt;br /&gt;
            junks=['.flux','.model','.psf','.residual','.mask','.image','.pb','.image.pbcor','.sumwt']&lt;br /&gt;
            for junk in junks:&lt;br /&gt;
                if os.path.exists(imname+junk):&lt;br /&gt;
                    os.system('rm -rf '+imname+junk)&lt;br /&gt;
        else:&lt;br /&gt;
            print('fits file '+fitsfile+' already exists, skip clean...')&lt;br /&gt;
            fitsfiles.append(fitsfile)&lt;br /&gt;
&lt;br /&gt;
    fig = plt.figure(figsize=(8.4,7.))&lt;br /&gt;
    gs = gridspec.GridSpec(5, 6)&lt;br /&gt;
    for s,sp in enumerate(spws):&lt;br /&gt;
        cfreq=cfreqs[int(sp)]&lt;br /&gt;
        ax = fig.add_subplot(gs[s])&lt;br /&gt;
        eomap=smap.Map(fitsfiles[s])&lt;br /&gt;
        eomap.plot_settings['cmap'] = plt.get_cmap('jet')&lt;br /&gt;
        eomap.plot(axes = ax)&lt;br /&gt;
        eomap.draw_limb()&lt;br /&gt;
        ax.set_title(' ')&lt;br /&gt;
        ax.get_xaxis().set_visible(False)&lt;br /&gt;
        ax.get_yaxis().set_visible(False)&lt;br /&gt;
        ax.set_xlim(xran)&lt;br /&gt;
        ax.set_ylim(yran)&lt;br /&gt;
        plt.text(0.98,0.85,'{0:.1f} GHz'.format(cfreq/1e9),transform=ax.transAxes,ha='right',color='w',fontweight='bold')&lt;br /&gt;
    plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)&lt;br /&gt;
    plt.show()&lt;br /&gt;
&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;br /&gt;
&lt;br /&gt;
[[file:Slf_t0_n0.png|thumb|center|500px|Figure 3: Images before self-calibration]]&lt;br /&gt;
[[file:Slf_t0_n1.png|thumb|center|500px|Figure 4: Images after self-calibration]]&lt;br /&gt;
&lt;br /&gt;
====Step 5: Quick-look imaging ==== &amp;lt;!-- ====Step 6: Making the full set of images====  --&amp;gt;&lt;br /&gt;
For further analysis of the event, follow the link given below.&lt;br /&gt;
https://colab.research.google.com/drive/1lSLLxgOG6b8kgu9Sk6kSKvrViyubnXG6?usp=sharing#scrollTo=xbXyyLmCFCGL&lt;br /&gt;
&lt;br /&gt;
&amp;lt;pre style=&amp;quot;background-color: #FCEBD9&amp;quot;&amp;gt;&lt;br /&gt;
from suncasa.utils import qlookplot as ql    ## (Optional) Supply the npz file of the dynamic spectrum from previous step.&lt;br /&gt;
				            ## If not provided, the program will generate a new one from the visibility.&lt;br /&gt;
                                            ## set the time interval&lt;br /&gt;
from suncasa import dspec as ds&lt;br /&gt;
import time&lt;br /&gt;
&lt;br /&gt;
specfile = visibility_data + '.dspec.npz'&lt;br /&gt;
d = ds.Dspec(visibility_data, bl='4&amp;amp;9', specfile=specfile)&lt;br /&gt;
&lt;br /&gt;
timerange = '19:02:00~19:02:10'             ## select frequency range from 2.5 GHz to 3.5 GHz&lt;br /&gt;
spw = '2~5'                                 ## select stokes XX&lt;br /&gt;
stokes = 'XX'                               ## turn off AIA image plotting, default is True&lt;br /&gt;
plotaia = False&lt;br /&gt;
ql.qlookplot(vis=visibility_data, specfile=specfile, timerange=timerange, spw=spw, \&lt;br /&gt;
stokes=stokes, plotaia=plotaia, restoringbeam=['60arcsec'], robust=0.5)&lt;br /&gt;
&amp;lt;/pre&amp;gt;&lt;/div&gt;</summary>
		<author><name>Sshaik</name></author>
	</entry>
</feed>