Visit ASIAA Homepage Registration Deadline: August 17, 2008
Second Asian Radio Astronomy School
August 18-22, 2008
Taipei, Taiwan

SMA Lab Instructions

INTRODUCTION   PART I   PART II   PART III   PART IV   PART V   PART VI   PART VII   PART VIII   REFERENCES

Welcome!

The purpose of this series of lab sessions is to familiar you with fundamental data and imaging aspects in interferometric observations. By following the instruction set, you will get some hands-on experiences with real data. At each step, a few questions are posted for you to think through in order to understand the concept. Do ask the session instructor or teaching assistants if you run into any question or problem.

INTRODUCTION


MIR is the software package to reduce data taken with the SMA.  The MIR package is originally developed by N. Scoville at Caltech. The MIR is written in IDL (which is the complete computing environment for the interactive analysis and visualization of data, widely used in the astronomical community) based on the OVRO calibration package MMA. While MIRIAD stands for Multichannel Image Reconstruction, Image Analysis and Display and was initially developed for data reduction and imaging for the Berkely Illinois Maryland Association (BIMA) Array. MIRIAD is now frequently used in the millimeter astronomy community.

For demonstration purpose, we will in this lab exercise of SMA use the MIR package to calibrate the SMA data (Wednesday afternoon) and then use MIRAID to make maps (Thursday afternoon).

Before you type idl, change your path to ~/SMA_DATA. All the files saved will be in this directory unless you change your default directory within idl. Type
> cd SMA_DATA
> idl
to get in idl environment.  You should see a message similar to the one blow prompting out,
IDL Version 6.2 (linux x86_64 m64). (c) 2005, Research Systems, Inc.
Installation number: 60400.
Licensed for use by: ASIAA
% Compiled module: PRO_MIR_INIT.
% Compiled module: STRSPLIT.
% Compiled module: PRO_SETUP_ENV.
% Compiled module: PRO_CONSTANTS.
idl_startup completed
IDL>
In following sections, you will need to inspect data in graphics. For IDL on linux platforms, the color table is used differently from the solaris ones. To get the correct color displayed for linux users, try
     IDL> device, decomposed=0, retain=2
before you plot.

pro_hlp
prints help comments from beginning of .pro files.
pro_lkf
looks for a routine with a keyword in its first comment line.
mirhelp
provides brief online description of several routines and wrappers related to the SMA data-reduction.

For example:
     IDL> pro_hlp,variable=<'variable_name'>
 to get help on a variable.
     IDL> pro_hlp,structure=<'structure_name'>
 to get help on a structure.
     IDL> pro_hlp,/program
 to get help on a program and print out the basic
introduction and usage of the function.
     IDL> mirhelp,'<task name>'
 to get help on a wrapper program, e.g. gain_cal etc.

And you can find the MIR Cookbook here.

PART I:  Data Inspection & Editing


GOAL : visibilities of continuum data, uv-coverage, and system temperature

First check if you are at the SMA_DATA directory.

In this part, we inspect a real observational dataset taken with the Submillimeter Array (SMA) to get you familiar with the SMA data structure and the concept of system temperature.

The first step is to load the dataset into idl.
IDL> mir_restore,'data_hh212'
IDL> select,/reset,/pos_wt
IDL> device, decomposed=0, retain=2
Then we can inspect the whole continuum data using the command "plot_continuum". 
IDL> plot_continuum
you should see a message similar to the one blow prompting out,
IDL> plot_continuum
*** mouse sounds :
 left button => select panel to zoom
middle button => keyboard options (print, etc)
right button => unzoom or go to next page
******************
squeak !
Click left button to zoom the panel
Click the middle button for more options
The above information tells you how to control the graphic displays by using the mouse. If you run into any problem, please ask TAs. 

Calling plot_continuum with no arguments will default to plot amp,pha vs hours for each integration w/ separate panels for combinations of baseline, receiver, sideband and continuum band, and a maximum of 4 plots per page. The sources will appear in different colors.  Several options can be definded in "plot_continuum" :

IDL> plot_continuum,x_var='hours',y_vars='amp,pha',
frame_vars='blcd,rec,sb,band', color_vars='source',
symbol_vars=0,frames_per_page=4

For example, we can plot all continuum visibilities together in a page:
 IDL> plot_continuum, frames_per_page='15'

[Question] How many antennas and baselines are present in the dataset? How many sources are there in the dataset ? 

There are five source in the dataset:  3c454.3,  Uranus,  0530+135,  HH212_N,  and HH212_S.  HH212_N and HH212_S are the target sources. While 3c454.3 is a strong quasar, and 0530+135 is a quasar close to the target source HH212.

Actually you can aslo use the command "select" to know some basic  observational information.
IDL> select,/reset,/pos_wt,/display

The command "select,/reset,/pos_wt" is for us to select all good data with positive weights. We can also use this command to select part of data in the dataset. Note that putting the keyword /reset in the command will reset the filter, i.e., ignoring all the previous selection results.    

How about taking a closer look at the data of 0530+135, HH212_N, and HH212_S from the baseline only between antenna "2" and antenna "3"?
IDL> select,/reset,/pos_wt,source=['0530+135','HH212_N','HH212_S'],baseline=['2-3']
IDL> plot_continuum
How about taking a closer look at the data of uranus and 3c454.3 from the same baseline (i.e., 2-3)?
IDL> select,/reset,/pos_wt,source=['3c454.3','uranus'],baseline=['2-3']
IDL> plot_continuum

[Question] Why do we need to observe 3c454.3, uranus, and 0530+135 in addition to the two target sources HH212_N and HH212_S? Why do we observe 3c453.3 and uranus only at the beginning of the observations, but observed 0530+135 much often (i.e., every 20 mins or so) ? 
 
Flux calibrator : Uranus
Bandpass calibrator : 3c454.3      
Gain calibrator : 0530+135

Now we would like to inspect the uv-coverage of the quasar "0530+135" and the target "HH212_N". We  will use the task  "plot_var".
IDL> select,/reset,/pos_wt,source=['0530+135']
IDL> plot_var,x='u',y='v',frame_vars='sb'
IDL> select,/reset,/pos_wt,source=['HH212_N']
IDL> plot_var,x='u',y='v',frame_vars='sb'

[Question] Given the data taken with the SMA array, can you guess what the source location is (positive, zero, or negative in declination angle)?


[Question] What do you see in the plots? Why the system temperatures vary in the way you see?

You can try the same command but on a single source source only.  Here are two examples, the target source HH212_N and the gain calibrator 0530+135, with x_axis of elevation  and hour angle :
IDL> select,/reset,/pos_wt,source='HH212_N'
IDL> plot_var,x='el',y='tssb',frames_per_page='15'
IDL> plot_var,x='ha',y='tssb',frames_per_page='15'
IDL> select,/reset,/pos_wt,source='0530+135'
IDL> plot_var,x='el',y='tssb',frames_per_page='15'
IDL> plot_var,x='ha',y='tssb',frames_per_page='15'

System temperature is a good measure of overall system performance. Tsys values are important in SMA data calibration since they are directly related to the visibility weights which scale as 1/Tsys2. The SMA does not automatically scale the data by Tsys as would be done at most other radio telescopes. The amplitude out from correlator is not corrected by the attenuation of the Earth atmosphere. Thus, the visibility amplitude needs to be multiplied by system temperatures. After you make sure the tsys values look fine, you can weight the data by Tsys:

IDL> select,/pos_wt,/reset
IDL> apply_tsys

Here we save the Tsys-corrected data to a new file.

IDL> select,/reset,/pos_wt
IDL> mir_save,'data_hh212_tsys',/nowait

PART II :  Bandpass Calibration


GOAL : SMA spectral-line data format, bandpass calibrators and techniques of bandpass calibration

Previously, we have investigated the continuum data of a real observational dataset. It was shown that observed visibilities in fact are not really a mere Fourier transform of the sky intensity distribution, as they are theoretically formulated. In Part II and III, we will further investigate these behavior and try to eliminate the "undesired distortion" of the data by applying certain calibrations.

Let us first inspect the visibilities of the spectral line data in the dataset.  We will use the command "plot_spectra" to plot the amplitude and phase of the data as a function of channel number (i.e., frequency).

IDL> select,/reset,/pos_wt,source='3c454.3'
IDL> plot_spectra,frames_per_page='15'
We can also plot the spectra, with x_corod in sky frequency :
IDL> plot_spectra,frames_per_page='15',x_var='fsky'

How about taking a closer look at the data from the baseline only between antenna "2" and antenna "3"?

IDL> select,/reset,/pos_wt,source='3c454.3',baseline=['2-3']
IDL> plot_spectra,color_vars='band'
IDL> plot_spectra,x_var='fsky',color_vars='band'
How about taking a further look at the data from only the first two "windows" in the spectra domain of the baseline 2-3?
IDL> select,/reset,/pos_wt,source='3c454.3',baseline=['2-3'],band=['s01','s02']
IDL> plot_spectra,color_vars='band'
IDL> plot_spectra,x_var='fsky',color_vars='band'

The SMA correlator processes two IF sidebands separated by 10 GHz, each of 1.968 GHz width (see figure below). The upper sideband (USB) and the lower sideband (LSB) are each divided into 6 "blocks" of 328 MHz width, and each block is further divided into 4 slightly overlapping "chunks" of 104 MHz width, where the usable bandwidth spans the central 82 MHz of each chunk (the sensitivity rolls off sharply at the chunk edges). The chunks are labeled s01, s02,...,s24. Note that the same correlator configuration applies to both the upper sideband and the lower sideband. The blocks are movable within the IF band, though there is very little experience with this feature so far.

[Question] How many chunks in each sideband ? What is the bandwidth of each chunk ? Do all chunks have the same channels ?

  Answer : 24, ~100 MHz, no

[Question] Can you examine the channel number in each chunks ? Why do we need large channel number in these particular chunks ?

  Answer :  256 channels for s02,s06,s12,s15,s16,s20,s21, and 64 channels for others

[Question] Do those visibilities (amplitudes and phases) look flat (independent of frequency)? Why?

The sources "3c454.3" and "uranus" in the observation run is in fact to be used as "bandpass calibrator" sources. That is because we actually "know" how their visibilities should look like, given we know it should have no spectral feature. Therefore, any deviation in the measured data from the ideal case can be attributed to the "band pass" function. By differencing the observed and "theoretical" values, we can derive the "passband" solution as a function of frequency (channel).

As you already inspected, the spectral-line data contain two "kinds" of chunks, with channel number of 64 and 256.  When deriving bandpass solutions, usually a few channels are smoothed to ehance the S/N.  To do bandpass calibration, therefore, it is better to separate the data into two parts, one conatning 64-channel chunks and the other 256-channel chunks, to do the bandpass calibration. Furthermore, since uranus is also a qualified bandpass calibrator, you can use not only the quasar 3c454.3 but also uranus to derive the bandpass solutions. 

To derive the bandpass solutions for 256-channel chunks:

IDL> select,/pos_wt,/reset,band=['s21','s20','s16','s15','s12','s06','s02']
IDL> pass_cal,smoothing=10,frames_per_page=60,ntrim=10
then you should see a message similar to the one blow prompting out,

Set which sources will be used as calibrators and set the calibrator fluxes
The flux should be in Jy at the frequency of the observation
These are the sources and their current passband codes
name: uranus passband cal: NO
name: 3c454.3 passband cal: NO
name: 0530+135 passband cal: NO
name: HH212_N passband cal: NO
name: HH212_S passband cal: NO
Enter source and new cal code. eg: 3C273 YES
or hit Return if all the sources are correctly specified
: uranus yes    <- enter "uranus yes" here

These are the sources and their current passband codes
name:  uranus             passband cal: YES
name:  3c454.3            passband cal: NO
name:  0530+135           passband cal: NO
name:  HH212_N            passband cal: NO
name:  HH212_S            passband cal: NO
Enter source and new cal code. eg: 3C273 YES
or hit Return if all the sources are correctly specified
: 3c454.3 yes     <- enter "3c454.3 yes" here

These are the sources and their current passband codes
name:  uranus             passband cal: YES
name:  3c454.3            passband cal: YES
name:  0530+135           passband cal: NO
name:  HH212_N            passband cal: NO
name:  HH212_S            passband cal: NO
Enter source and new cal code. eg: 3C273 YES
or hit Return if all the sources are correctly specified
:       <-  just hit Return here

Afterward you should see two plots  (here and here) of the derived bandpass curves. If you are satisfied with the solutions, you can apply the solutions.

Stacking gain curves to the ca structure
Apply passband solution? [NO <YES>]: yes <-  enter "yes" here

Similar to that for 256-channel chunks, we now do bandpass calibration for the chunks with 64 channels. Here we average 4 channels and ignore the information of the edge 4 channels (due to their poor responses).
IDL> select,/pos_wt,/reset,band=['-s21','-s20','-s16','-s15','-s12','-s06','-s02']
IDL> pass_cal,smoothing=4,frames_per_page=48,ntrim=4
again, '3c454.3' and 'uranus' are chose as bandpass calibrators.

After applying the bandpass solutions, we can inspect the bandpass-calibrated dataset:
IDL> select,/reset,/pos_wt,source='3c454.3'
IDL> plot_spectra,frames_per_page='15'
IDL> select,/reset,/pos_wt,source='3c454.3',baseline=['2-3']
IDL> plot_spectra,color_vars='band'
IDL> plot_spectra,x_var='fsky',color_vars='band'
IDL> select,/reset,/pos_wt,source='3c454.3',baseline=['2-3'],band=['s01','s02']
IDL> plot_spectra,color_vars='band'
IDL> plot_spectra,x_var='fsky',color_vars='band'
We can save the bandpass-calibrated data to a new file.

IDL> select,/reset,/pos_wt
IDL> mir_save,'data_hh212_bp',/nowait

Finally, please compare the line features seen in the target source data before and after the bandpass calibration.
IDL> select,/reset,/pos_wt,source='HH212_N',baseline=['2-3'],band=['s02']
IDL> plot_spectra,x_var='fsky',color_vars='band'
To plot the spectra before bandpass calibration, you have to open another xterm, run idl, and reload the data wo bandpass calibration. Please run the following commnads in a new terminal.
IDL> mir_restore,'data_hh212_tsys'
IDL> select,/reset,/pos_wt
IDL> device, decomposed=0, retain=2
IDL> select,/reset,/pos_wt,source='HH212_N',baseline=['2-3'],band=['s02']
IDL> plot_spectra,x_var='fsky',color_vars='band'
[Question] What features do you see in the target source data HH212_N before and after the bandpass calibration?

PART III :  Flux and Gain Calibration


GOAL : flux calibrators, flux calibration, gain calibrators and techniques of gain calibration

Take the "bandpass-calibrated" dataset, Let us further investigate the phase and amplitude behavior verses time during an observation run. Recall how to plot out the amplitudes and phases of a visibility dataset?  Let us inspect the amplitudes and phases of visibilities of the gain calibrator 0530+135.

IDL> select,/reset,/pos_wt,source=['0530+135']
IDL> plot_continuum,frames_per_page='15'
How about taking a closer look at the data from the baseline only between antenna "2" and antenna "3"?
IDL> select,/reset,/pos_wt,source=['0530+135'],baseline=['2-3']  
IDL> plot_continuum

[Question] Do the amplitude data show flat features? What about the phases? 
How should the visibilities look like for this source in an ideal case? Do you agree that the data look consistent with the source being a point source? Where does these inconsistency come from?

IDL> select,/reset,/pos_wt,source=['0530+135','HH212_N'],baseline=['2-3']
IDL> plot_continuum

[Question] How do the two sources get observed? Do the phases of two sources follow each other? What do you expect? Why?

Similar to the previous bandpass calibration case, the source 0530+135 in this observation run is in fact to be used as a "gain calibrator" source. That is because we actually "know" how its visibisities should look like, given we know the source structure (and what that is?). Therefore, any deviation in the measured data from the ideal case can be attributed to the "gain" function. Again, by differencing the observed and "theoretical" measurements, we can derive the gain solutions as a function of time.

Here we will use the command "gaincal" to derive and apply the gain solutions. Before we do this, however, we need to determine the absolute flux scale of the gain calibrator 0530+135 first by using the commands "sma_flux_cal" and "flux_measure".

IDL> select,/reset,/pos_wt
IDL> sma_flux_cal

you should see a message similar to the one blow prompting out,

IDL> sma_flux_cal
% Compiled module: SMA_FLUX_CAL.
% Compiled module: SMA_FLUX_CAL_INI.
       15165 passed in list
% Compiled module: UTI_DISTINCT.
% Compiled module: UNIQ.
% Compiled module: FLUX_PRIMARY.
% Compiled module: UTI_INTERP.
% Compiled module: FLUX_SECONDARY.
#   Source   Flags   Nscans      Amp        Flux(Jy)
    uranus             94         21.        65.18
   3c454.3             90         5.3        0.000 
0530+135            139         2.0        0.000
   HH212_N            347        0.96        0.000
   HH212_S            341        0.93        0.000
Enter flux calibrator source, and if needed, flux in Jy, eg: 3c279 18.1
: urnaus 65.18 <- enter "urnaus 65.18" here


then you should see a table of the derived scaling fator like,

  all        sidebands           baselines

% Compiled module: UTI_MEANVAR.
  1.33E+00   l    1.33E+00   l  1-2      9.64E-01
                             l  1-3      1.24E+00
                             l  1-5      1.41E+00
                             l  2-3      1.21E+00
                             l  2-5      1.26E+00
                             l  2-6      1.13E+00
                             l  2-7      1.46E+00
                             l  3-5      1.22E+00
                             l  3-6      1.26E+00
                             l  3-7      1.66E+00
                             l  5-7      1.92E+00
                             l  6-7      1.42E+00
select which scale factors to scale data:
1. sidebands; 2. baselines; other. single(all)
: 2 <- enter "2" here (i.e., using different scaling factor for different baselines)

then you have to decide whether applying the soluitions or not,

% Compiled module: CAL_STORE.
% Compiled module: DAT_COMB_SEP.
Stacking gain curves to the ca structure
Apply Flux Calibration? [NO <YES>]:  yes <- enter "yes" here
           The next step is to measure the flux denisty of the gain calibrator 0530+135 using the command "flux_measure"
IDL> select,/reset,/pos_wt,source=['uranus','0530+135']
IDL> flux_measure

Either scalar or vector average can be used to determine the flux of the quasar, here we choose scalar average
(The vector average of a series of complex numbers is obtained by averaging the real and imaginary part separately. For scalar average, amplitudes and phases are averaged)

% Compiled module: FLUX_MEASURE.
% Compiled module: SMA_FLUX_MEASURE.
      379125 passed in list
% Compiled module: MEAN.
% Compiled module: MOMENT.
Scalar Average or Vector Average ? [S <V>]: s <- enter "s" here
          then you can see the flux of the gain calibrator 0530+135:
Scalar average: 
#   Source   Flags   Nscans  Flux(Jy)   SNR    meantime    REAL       IMAG
  0530+135            139    2.5910     127     11.20      0.0545      0.2932

Okay, now we can use the command "gaincal" to derive the gain solutions. We choose ant 2 as reference antenna to derive antenna-based solutions 

IDL>select,/reset,/pos_wt,source=['0530+135','HH212_N','HH212_S']
IDL>gain_cal,x='hours',poly=4,/preavg,tel_bsl='telescope',refant=2,frames_per_page=6
% Compiled module: GAIN_CAL.
-------------------------------------------
      Single Band Calibration Mode        
          Main Band Calibration           
-------------------------------------------
% Compiled module: GAIN_INI.
Set which sources will be used as calibrators and set the calibrator fluxes
The flux should be in Jy at the frequency of the observation
These are the sources and their current gain codes
name:  0530+135                gain code: NO   flux (Jy):       0.0000000
name:  HH212_N                 gain code: NO   flux (Jy):       0.0000000
Enter source, cal code, and if cal, flux in Jy, eg: 3C273 YES 3.1
or hit Return if all the sources are correctly specified
: 0530+135 yes 2.6 <- enter "0530+135 yes 2.6" here


These are the sources and their current gain codes
name: 0530+135 gain code: YES flux (Jy): 2.5999999
name: HH212_N gain code: NO flux (Jy): 0.0000000
name: HH212_S gain code: NO flux (Jy): 0.0000000
Enter source, cal code, and if cal, flux in Jy, eg: 3C273 YES 3.1
or hit Return if all the sources are correctly specified
:
<- just hit Return here

Afterward you should see a plot of derived gain curves. If you are satisfied with the solutions, you can apply the solutions.

EXITING FROM THE INTERACTIVE MODE
% Compiled module: GAIN_STORE.
% Compiled module: CAL_STORE.
% Compiled module: DAT_COMB_SEP.
Stacking gain curves to the ca structure
% Compiled module: PLO_PLID_REL.
Apply gain solution? [NO <YES>]: yes <- enter "yes" here

Then we will use "mir_save" to store the well-calibrated data. 

IDL> select,/reset,/pos_wt
IDL> mir_save,filename='data-hh212_final',/nowait
Now, looks at the amplitudes and phases of the calibrator source again,
IDL> select,/reset,/pos_wt,source=['0530+135']
IDL> plot_continuum,frames_per_page='15'
IDL> select,/reset,/pos_wt,source=['0530+135'],baseline=['2-3']
IDL> plot_continuum
[Question] Do they look more like expected values for a point source
IDL> select,/reset,/pos_wt,source=['0530+135','HH212_N']  
IDL> plot_continuum,frames_per_page='15'
IDL> select,/reset,/pos_wt,source=['0530+135','HH212_N'],baseline=['2-3']
IDL> plot_continuum
[Question] How do the amplitudes and phases for the two sources look like now?

At this point, our dataset is calibrated and ready to be imaged. We need transfer the data from MIR format to Miriad format.

IDL> select,/reset,/pos_wt
IDL> idl2miriad,source='0530+135',sideband='l',dir='0530+135'
IDL> idl2miriad,source='HH212_N',sideband='l',dir='HH212N'
IDL> idl2miriad,source='HH212_S',sideband='l',dir='HH212S'

Tomorrow afternoon you will learn how to make maps with the calibrated visibilities.


PART IV: Checking Visibility with MIRIAD


GOAL: Inspecting visibility data with MIRIAD

Yesterday we learned how to calibrate the visibility data. Today we will follow with imaging process using MIRIAD.

Before we start imaging, let us take the calibrated data set and have a few more looks on the visibilities with MIRIAD. Let us first investigate the phase and amplitude behavior of 0530+135 verses time during an observation run.

  > uvplt vis=0530+135 device=/xw nxy=5,3 line=wide,1,1 axis=time,amp
  > uvplt vis=0530+135 device=/xw nxy=5,3 line=wide,1,1 axis=time,pha yrange=-180,180
[Question] How should the visibilities look like for this source in an ideal case?

The source 0530+135 in this observation run is in fact to be used as a "gain calibrator." That is because we actually "know" how its visibilities should look like, given we know the source structure (and what that is?). Therefore, any deviation in the measured data from the ideal case can be attributed to the "gain" function.

Now look at the amplitudes and phases of the calibrator source again as a function of the uv-distance,

  > uvplt vis=0530+135 device=/xw line=wide,1,1 axis=uvdis,amp options=nobase 
  > uvplt vis=0530+135 device=/xw line=wide,1,1 axis=uvdis,pha yrange=-180,180 options=nobase
[Question] Do they look more like expected values for a point source?

Let us further investigate the frequency dependence of phase and amplitude. Plot out the amplitudes and phases of a visibility dataset in the frequency domain for the shortest baseline, 2-3. We expect the visibilities to be in fact flat as quasars normally have no spectral feature.

  > smauvspec device=/xw select='ant(2)(3)' interval=10000 nxy=1,1 \
    options=avall hann=3 vis=0530+135 axis=freq,amp
  > smauvspec device=/xw select='ant(2)(3)' interval=10000 nxy=1,1 \
    options=avall hann=3 vis=0530+135 axis=freq,pha
[Question] Do they look more like expected values for a quasar?

Let us also plot the freqency dependence of phase and amplitude for the source HH212.

  > smauvspec device=/xw select='ant(2)(3)' interval=10000 nxy=1,1 \
    options=avall hann=3 vis=HH212N axis=freq,amp
  > smauvspec device=/xw select='ant(2)(3)' interval=10000 nxy=1,1 \
    options=avall hann=3 vis=HH212N axis=vel,amp
[Question] What features do you see in the target source data HH212?

The source data has its rest freqeuency set to 345.78500 GHz, which is different from that of our target line SiO (8-7) at 347.330631 GHz. Such a "wrong" rest frequency translates to incorrect velocity for our target line. Let us just correct the rest frequency to 347.330631 GHz for SiO (8-7).

  > puthd in=HH212N/restfreq value=347.330631
  > puthd in=HH212S/restfreq value=347.330631

Check the rest frequency value and plot the spectra again.

  > gethd in=HH212N/restfreq 
  > smauvspec device=/xw select='ant(2)(3)' interval=10000 nxy=1,1 \
    options=avall hann=3 vis=HH212N axis=vel,amp
[Question] Can you find the chunk with the SiO line?

At this point, we are ready for imaging.

PART V: Imaging a Point Source


GOAL: Various imaging aspects including Fourier transform, weighting schemes, de-convolution, reconstruction, etc. Case study of the calibrator 0530+135.

Following the previous parts, we finally have "calibrated" visibilities to make images of our observed sources. We first use INVERT to demonstrate how Fourier transform is accomplished, than use CLEAN and RESTOR for de-convolution and restoration.

Let us start with imaging the continuum emission of the calibrator to build some sense of the whole process. First, we perform an "inverse" Fourier transform with INVERT, which converts the visibilities to a "dirty" map as well as the sampling function to a synthesized beam.

 > invert vis=0530+135 map=0530+135.map beam=0530+135.beam \
   line=wide,1,1 sup=0 cell=0.35 imsize=97 options=systemp  

To get the best signal-to-noise ratio, natural weighting is applied with keyword "sup=0." At this step, two images are created, one is the source brightness distribution, often called "dirty map", and another one is the PSF (point-spread-function), often called "synthesized beam". To look at the two images,

  > cgdisp in=0530+135.map device=/xw options=full labtyp=arcsec
  > cgdisp in=0530+135.beam device=/xw options=full labtyp=arcsec
[Question] Does the target source really look like a point source? What is the maximum value in the dirty map? Can you explain all that you see?

We have to play with the "dirty" map to make it "cleaner" and reflect better what the real source brightness distribution look like. CLEAN tries to remove the PSF response in the dirty map by finding out where the source is likely to be, assuming it is a point source with certain intensity, 1) registering that point source intensity into the so called clean component or model image, and 2) subtracting the contribution of this component from the dirty map. This contribution, as you should realize, is the clean component convolved with the PSF. After the subtraction, the remaining image is called a "residual" image, we can then repeat this "clean" process with the new "residual" image until at some point we think there is no more "source" needs to be extracted from the residual map.

First try clean "one step (niters=1"),

   > clean map=0530+135.map beam=0530+135.beam out=0530+135.mod1 niters=1

Take a look at the clean component image and find out what the maximum value is

  > cgdisp in=0530+135.mod1 device=/xw options=full labtyp=arcsec 

Now we create the residual image and take a look at it.

  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod1 \
    mode=residual out=0530+135.res1
  > cgdisp in=0530+135.res1 device=/xw options=full labtyp=arcsec
[Question] How does the residual map look like? What is the maximum value in the residual map?

As mentioned in the class, a "restored" map, which is often called "clean" map, is created in the following way, 1) The "clean" beam is created by taking the inner (main) peak of the synthesized beam (PSF) and fitting it with a 2D Gaussian. 2) The "restored" image is generated by convolving the clean component with the "clean" beam (an ideal image which hopefully represent the source brightness distribution) and adding back the residual map (to the former ideal image).

  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod1 \
    out=0530+135.cm1
[Question] What do you notice in the terminal output?

As you can see, the Gaussian beam FWHM and position angle are the results from Gaussian fitting of the synthesized beam. Look at the "restored" image

  > cgdisp in=0530+135.cm1 device=/xw options=full,beambl labtyp=arcsec

or try

  > cgdisp in=0530+135.cm1 device=/xw type=c options=full,beambl labtyp=arcsec
[Question] Do you notice a ellipse at the bottom-left corner of the image? Can you identify what that means?

You should find that this kind of "restored" image is far from satisfactory. There are significant effects from the PSF (these defects are called "sidelobes", as the conventional terminology would refer to). That means, we did not clean "deep" enough. In the following two sequences, we will clean "10 iterations" and "50 iterations".

SEQUENCE A:

  > clean map=0530+135.map beam=0530+135.beam out=0530+135.mod10 niters=10
  > cgdisp in=0530+135.mod10 device=/xw options=full labtyp=arcsec
  >
  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod10 \
    mode=residual out=0530+135.res10
  > cgdisp in=0530+135.res10 device=/xw options=full labtyp=arcsec
  >
  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod10 out=0530+135.cm10
  > cgdisp in=0530+135.cm10 device=/xw options=full
  > cgdisp in=0530+135.cm10 device=/xw type=c options=full

SEQUENCE B:

  > clean map=0530+135.map beam=0530+135.beam out=0530+135.mod50 niters=50
  > cgdisp in=0530+135.mod50 device=/xw options=full labtyp=arcsec
  >
  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod50 \
    mode=residual out=0530+135.res50
  > cgdisp in=0530+135.res50 device=/xw options=full labtyp=arcsec
  >
  > restor map=0530+135.map beam=0530+135.beam model=0530+135.mod50 out=0350+135.cm50
  > cgdisp in=0530+135.cm50 device=/xw options=full
  > cgdisp in=0530+135.cm50 device=/xw type=c options=full
[Question] What did you notice in these two sequence of outputs?

As in the first trial sequence, you should pay attention to how the maximum values in those "clean component", "residual", and "restored" images change, as well as how the levels of those "sidelobes" change.

[Question] Think about how to decide how deep or when to stop "clean"ing?

The proper iteration number can be quite different from map to map. One way to get a handle with CLEAN is to relax the niters number and set a "cutoff" limit, where deconvolution stops when the extreme values of the residual map goes below the cutoff value. In practive, a cutoff value of roughly 2-3 times as large as the theoretical rms value is usually good enough to produce reasonable clean components. To speed up the rest of this exercise, we will just follow this simple guideline.

Let us read out the rms value from the dirty map and set a cutoff 3 times as large as the theoretical rms value for CLEAN.

  > gethd in=0530+135.map/rms
  > clean map=0530+135.map beam=0530+135.beam model=0530+135.model \
    niters=1000000 cutoff=0.07
  
  > restor map=0530+135.map beam=0530+135.beam model=0530+135.model \
    mode=residual out=0530+135.res 
  > cgdisp device=/xw in=0530+135.res options=full

  > restor map=0530+135.map beam=0530+135.beam model=0530+135.model \
    out=0530+135.cm 
  > cgdisp device=/xw in=0530+135.cm options=full

PART VI: Single-Channel Mosaic


GOAL: Combining multiple fields (pointings) into a mosaicked image. Case study of the continuum emission of HH212.

Some astronomical phenomena appear on the sky with angular extents larger than the observed field of view (FoV), which in our case is set by the diffraction limit of one single element. When this happens, we will need to observe the desired region with more than one pointings and linearly "add" multiple fields to make one mosaicked image. The SiO protostallr jet of HH212 is one of these cases.

Before we start the mosaicking, let us try to image one single field and examine whether CLEAN has been properly "configured."

  > invert vis=HH212N map=HH212Nd.map beam=HH212Nd.beam line=wide,1,1 \
    cell=0.35 imsize=97 sup=0 options=systemp
  > gethd in=HH212Nd.map/rms
  > clean out=HH212Nd.model map=HH212Nd.map beam=HH212Nd.beam niters=100000 \
    cutoff=0.043
  > restor out=HH212Nd.res map=HH212Nd.map beam=HH212Nd.beam model=HH212Nd.model \
    mode=residual
  > restor out=HH212Nd.cm map=HH212Nd.map beam=HH212Nd.beam model=HH212Nd.model 

Now let us compare the image sizes of the dirty map, beam, and clean component (model).

  > cgdisp device=/xw labtyp=arcsec options=full in=HH212Nd.map
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212Nd.beam
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212Nd.model
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212Nd.res
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212Nd.cm

[Question] Can you see any difference in the sizes of the maps? Why are they different? Is it reasonable to keep them the way they are?

For a point source at the center of the map, such as our calibrator, it is fine to have a clean model that is much smaller than the dirty map. However, in the case of the source, the emission can occur beyond the cnetral quarter, and CLEAN has to work across the entire mapping region. Let us make a synthesized beam twice larger than the dirty map and repeat the process again.

  > invert vis=HH212N map=HH212N.map beam=HH212N.beam line=wide,1,1 \
    cell=0.35 imsize=97 sup=0 options=systemp,double
  > gethd in=HH212N.map/rms
  > clean out=HH212N.model map=HH212N.map beam=HH212N.beam niters=100000 \
    cutoff=0.043
  > restor out=HH212N.res map=HH212N.map beam=HH212N.beam model=HH212N.model \
    mode=residual
  > restor out=HH212N.cm map=HH212N.map beam=HH212N.beam model=HH212N.model 

Check again the map sizes,

  > cgdisp device=/xw labtyp=arcsec options=full in=HH212N.map
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212N.beam
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212N.model
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212N.res
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212N.cm

[Question] How does every map compare with the one in the previous sequence in terms of quality and size?

To build any decent "clean" maps, one has to involve the keyword "options=double" all the time, escpeically when extended emission is present.

Now we are ready to make a mosaicked imaging, which can be easily accomplished by using "options=mosaic" along with "offset=05:43:51.404,-01:02:53.10" to register the continuum peak at the map center. Note that, MOSSDI, an improved de-convolution alogrithm, should be used for the "clean"ing of a mosaicked image.

  > invert vis=HH212N,HH212S map=HH212.map beam=HH212.beam line=wide,1,1 \
    cell=0.35 imsize=97 sup=0 options=systemp,double,mosaic \
    offset=05:43:51.404,-01:02:53.10
  > mossdi out=HH212.model map=HH212.map beam=HH212.beam niters=100000 \
    cutoff=0.043
  > restor out=HH212.res map=HH212.map beam=HH212.beam model=HH212.model \
    mode=residual
  > restor out=HH212.cm map=HH212.map beam=HH212.beam model=HH212.model 

Again, check your images with CGDISP,

  > cgdisp device=/xw labtyp=arcsec options=full in=HH212.map
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212.beam
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212.model
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212.res
  > cgdisp device=/xw labtyp=arcsec options=full in=HH212.cm,HH212.cm type=c,p 

PART VII: Multi-Channel Mosaic


GOAL: Generating images for SiO (8-7) line emission. Case study of the SiO protostellar jet from HH212.

Finally, we can apply the command sequence for mapping velocity channels of our target source data. As previosly mentioned, the datasets HH212N and HH212S include SiO (8-7) line emission in the second chunk (window). Let us make a series of mosaicked images with velocity step and width of 1 km/s.

  > invert vis=HH212N,HH212S map=sio.map beam=sio.beam cell=0.35 imsize=194 \
    sup=0 line=vel,39,-19.3,1,1 options=systemp,double,mosaic \
    offset=05:43:51.404,-01:02:53.10
  > mossdi out=sio.model map=sio.map beam=sio.beam niters=100000 cutoff=0.4 \
    region='arcsec,box(-10,-20,10,20)' 
  > restor out=sio.res map=sio.map beam=sio.beam model=sio.model mode=residual
  > restor out=sio.cm map=sio.map beam=sio.beam model=sio.model 
  > cgdisp in=sio.cm,sio.cm type=p,c device=/xw labtyp=arcsec \
    nxy=3,2 options=full,3val 3format='f10.1' csize=0,.8 olay=sio.olay 

[Question] The image size is given by the keyword imsize, which has a value twice larger than before. Can you identify the effect of this change and find the reason behind it?

[Question] Can you identify the blue-shifted (negative velocity) channels and the red-shifted (positive velocity) channels with respect to the systemic velocity of 1.7 km/s?

Now we can also make two velocity-integrated images, one for the blue-shifted component and the other for the red-shifted component. The blue-shifted componenet is approximately from -17.3 to 1.7 km/s, corresponding to a line width of 19 km/s centered at -7.8 km/s.

  > invert vis=HH212N,HH212S map=blue.map beam=blue.beam cell=0.35 imsize=194 \
    sup=0 line=vel,1,-7.8,19,19 options=systemp,double,mosaic \
    offset=05:43:51.404,-01:02:53.10
  > mossdi out=blue.model map=blue.map beam=blue.beam niters=100000 cutoff=0.4 \
    region='arcsec,box(-10,-20,10,20)' 
  > restor map=blue.map beam=blue.beam model=blue.model out=blue.cm
  > cgdisp device=/xw labtyp=arcsec in=blue.cm,blue.cm type=c,p

The red-shifted component is roughly from 1.7 to 14.7 km/s, corresponding to a line width of 13 km/s centered at 8.2 km/s.

  > invert vis=HH212N,HH212S map=red.map beam=red.beam cell=0.35 imsize=194 \
    sup=0 line=vel,1,8.2,13,13 options=systemp,double,mosaic \
    offset=05:43:51.404,-01:02:53.10
  > mossdi out=red.model map=red.map beam=red.beam niters=100000 cutoff=0.4 \
    region='arcsec,box(-10,-20,10,20)' 
  > restor map=red.map beam=red.beam model=red.model out=red.cm
  > cgdisp device=/xw labtyp=arcsec in=red.cm,red.cm type=c,p

Let us place the two velocity-integrated images together in one plot.

  > cgdisp device=/xw in=red.cm,blue.cm type=c,c options=full \
    beamtyp=b,l labtyp=hms,dms olay=sio.olay 

Alternatively, we may output the plot as a postscript file

  > cgdisp device=jet.ps/cps in=red.cm,blue.cm type=c,c options=full \
    beamtyp=b,l labtyp=hms,dms olay=sio.olay 

Overlaid with an H2 image (2.12 μm)

To compare obesrvations in different wavebands, we often plot several images together as a composite image. Here we will show an example of plotting our red-shifted and blue-shifted components on top of a near-infrared image showing the H2 emission around HH212. ▶download tar file

  > cgdisp in=HH212H2,red.cm,blue.cm type=p,c,c labtyp=arcsec options=full \
    beamtyp=b,l olay=sio.olay range=-0.1,13,lin,3 device=/xw 
  > cgdisp in=HH212H2,red.cm,blue.cm type=p,c,c labtyp=arcsec options=full \
    beamtyp=b,l olay=sio.olay range=-0.1,13,lin,3 device=brh2.ps/cps 

PART VIII: More Post Processing


GOAL: Analysizing the SiO data cube

You can run other MIRIAD tasks such as velplot or imspec to further investigate image cubes. For example, the MIRIAD command wrapper velplot can be used to interactively take spectra and position-velocity slices.

 velplot in=sio.cm

Beyond the SiO data cube

You may try to create a data cube for the CO (3-2) emission on your own. The rest frequency of CO (3-2) is at 345.795991 GHz.

REFERENCES


SMA MIR Website : http://cfa-www.harvard.edu/~cqi/mircook.html
MIRIAD Website : http://bima.astro.umd.edu/miriad
ASIAA will not contact participants for credit card information by phone. Privacy and Security Policy