Protoplanetary Disk Simulation - VLA-CASA5.3.0: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Line 244: Line 244:
</source>
</source>


When expressed in wavelengths, the uv-coverage and synthesized beam are somewhere in between that of Q-band A-configuration and C-configuration as expected (Fig. 13).


{|
{|
Line 254: Line 255:
</source>
</source>


The ring shape is still visible in the results (Fig. 14) although quite blurred. We note that the noise is considerably lower as expected from the lower band noise properties of the receiver and atmosphere. We adjusted the clean threshold accordingly.


{|
{|
|[[Image:VLAsim-psim8-ana.png|600px|thumb|left|Fig. 13: X-band simulation results.]]
|[[Image:VLAsim-psim8-ana.png|600px|thumb|left|Fig. 14: X-band simulation results.]]
|}
|}

Revision as of 19:08, 10 May 2018


The following tutorial shows how to adopt typical parameters for simulating Karl G. Jansky VLA data. We use the same image as the ALMA tutorial Protoplanetary Disk Simulation. Follow this link to obtain the protoplanetary disk model. Models are in units of Jy/pixel. Other simulation options, e.g. using component lists, or how to use the toolkit are available on the Simulations in CASA section of CASAguides.

In Fig. 1 we show the model that we use for this simulation tutorial.

Fig. 1: Model.

The ALMA version of the tutorial shows a way on how to derive the center of the image. We will use their results and specify direction='J2000 18h00m00.031s -22d59m59.6s' for all of our simulations. The image center can also be determined with the CASA viewer. Given that the VLA primary beams at the VLA frequencies are much larger than the image, the precise pointing direction center is less important.

We will use the simobserve and simanalyze tasks as in the ALMA tutorials. The model is derived for 672GHz, we will adapt the simulation to work at VLA frequencies.

Note that simobserve has a few limitations. E.g. it cannot simulate different spectral windows. If this is desired, one needs to set up the simulation for each spw separately, and then use concat to merge all MeasurementSets. simobserve also has no option to add pointing errors to the simulated data. All VLA configurations and the VLA receiver temperatures are, however, accessible in simobserve.


Q-band, 128MHz bandwidth, noiseless image, 1h integration time, A-configuration, no deconvolution

Let's start with a simulation at Q-band at 44GHz, with a bandwidth of 128MHz. We will use the VLA A-configuration as it provides the resolution that is needed for the disk to be well resolved. To start with we do not apply any noise to the data:

# In CASA
simobserve(project='psimvla1', skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel', incenter='44GHz', inwidth='128MHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='', graphics='both', overwrite=True)

project: The name of our project is psimalma1

skymodel: The input model image in Jy/pixel units. We adjust the peak to a lower [math]\displaystyle{ 3\times10^{-5} }[/math]Jy/pixel value with the inbright parameter, as expected at a lower frequency. We also overwrite the fits header to assume that the model is valid for 44GHz with incenter and the bandwidth to 128MHz with inwidth

setpointings: will allow simobserve to derive the pointign s by its own algorithm. Given that the primary beam at Q-band is about 1arcminutes (see the VLA [observational status summary], and the size of the model is less than an arcsecond, a single pointing will be more than enough.

integration: To avoid time smearing, we follow the guidance for data rates in the [OSS] and assume 2s per visibility integration.

direction: the center of the map. For a single pointing this is equivalent to the pointing center.

obsmode: int is used for interferometric data such as VLA observations.

antennalist: the VLA configuration antenna position file. The files are available in CASA via 'vla.x.cfg' where 'x' is the name of the array configuration. Here 'vla.a.cfg' is the VLA A configuration (the python command os.getenv("CASAPATH").split()[0]+"/data/alma/simmos/" shows the directory that contains all readily available configurations)

hourangle: is used to simulate observations at a specific hour angle. We use 'transit' for culmination.

totaltime: This is the time on source.

thermalnoise: We leave this parameter empoty for a this noise-less simulation.

graphics: 'both' will show graphics on the screen and save them as png files in the project directory.

overwrite: True will overwrite previous results: be careful when running multiple setups as the files may have different names and only the files with the same names will be overwritten.

The output of the simulation is shown in Figs. 2 and 3. The first image is the sky coverage which shows clearly that the primary beam exceeds the size of the model image by far.

Fig. 2: Sky Coverage.
Fig. 3: Output of simobserve. Upper left: Elevation vs time of the modeled source (blue). The vertical green line marks the transit time and the red bar is the time span of the simulated observations. Upper right: The VLA antenna in units of meters to the central VLA reference position. Lower left: The uv-coverage of the simulations (note that they do not account for the bandwidth), Lower right: the psf of the observation with fitted Gaussian major and minor beam sizes.


The task simanalyze can now procvess the nely created MeasurementSet and image the data:

# In CASA
simanalyze(project='psimvla1', image=True, vis='psimvla1.vla.a.ms',  imsize=[192, 192], interactive=False, niter=0, weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

project: The project name, same as in our call for simobserve.

image: True will image the visibilities with the following sub-parameters.

vis: The input MeasurementSet. Note that simobserve's notation is to use the project name appended by the antenna configuration file name.

imsize: the number of pixels for the image dimensions. We use 192 for each axis. By default simobserve will use the pixel size of the model image as the cell size for each pixel (3.11 milli-arcseconds in our case).

interactive: If set to True, the cleaning will be interactive to allow the step by step setting of clean boxes and clean depths (see the [VLA imaging topical guide]). We will turn this feature off.

niter: number of clean iterations. We start with a dirty map and set niter=0.

weighting: [Image weighting schemes]. We will use briggs which defaults to a robust = 0.5.

pbcor: If set to True the image will be primary beam corrected. In our case of a very large primary beam (see Fig. 1) primary beam correction will introduce very small corrections.

stokes: The Stokes polarization to be produced, we want Stokes I.

analyze: When set to True simobserve will perform basic analysis of the produced images and produce hard copies of the following displays.

showuv: We set this to False as we already show the uv-coverage in the {simobserve}} output (Fig. 3).

showpsf: Although we already have a plot of that in Fig. 3, we show it here again to have the comparison to the other plots.

showconvolved: produces a plot of the model convolved with the synthesized clean beam.

showclean: shows the deconvolved image (in our case we do not perform deconvolution, but show the plot anyways).

showresidual: a plot of the residual image after deconvolution.

showdifference: displays the difference between output cleaned image and input model sky image convolved with output clean beam, i.e. it shows remaining clean artefacts.

showfidelity: a fidelity image, as defined by [math]\displaystyle{ \frac{I}{|I-T|} }[/math] where [math]\displaystyle{ I }[/math] is the observed image and [math]\displaystyle{ T }[/math] the sky model, see the simanalyze CASAdoc help.

The resulting plots of our first simulation are shown in Fig. 4. Since this is a dirty image only (niter=0), the image and residual are the same dirty image. SInce we also did not turn on any thermal noise, the image corruptions are all due to the incomplete uv-coverage of our observations.

Fig. 4: Simulated Images for noise-less Q-band observations. Top right: the psf of the observations with the clean beam parameters. Top center: the Sky model. Top right: The sky model convolved with the clean beam. Bottom left: the simulated image. Bottom center: the residual image from deconvolution. Bottom right: the fidelity image.

Q-band, 128MHz bandwidth, noiseless image, 1h integration time, A-configuration, cleaned

As a second step we will deconvolve this data set. This could be done with inside the same project, but for a cleaner separation of our use cases, we will create a new MS with simobserve with a modified call of simanalyze:

# In CASA
simobserve(project='psimvla2', skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel', incenter='44GHz', inwidth='128MHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='', graphics='both', overwrite=True)

simanalyze(project='psimvla2', image=True, vis='psimvla2.vla.a.ms',  imsize=[192, 192], interactive=False, niter=1000, threshold='1e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

The only difference to the earlier call of simobserve is our setting of niter=1000, threshold='1e-4Jy' which corresponds to a maximum number of 1000 clean iterations and a clean threshold of 0.1mJy. The threshold was chosen based on the previous image to be a few times above the rms of the image. The residual image becomes flatter as expected and the bulk of the sidelobes disappear.

Fig. 5: Same as Figure 4, this time, however, with a cleaned image.

Q-band, 128MHz bandwidth, 4mm pwv, 1h integration time, A-configuration

We now add a precipitable water vapor of 4mm.

# In CASA
simobserve(project='psimvla3', skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel', incenter='44GHz', inwidth='128MHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='tsys-atm', user_pwv=4, t_ground=270.0, graphics='both', overwrite=True)

The simobserve output will be identical to Fig. 3 since we do not change any of the displayed parameters.

# In CASA
simanalyze(project='psimvla3', image=True, vis='psimvla3.vla.a.noisy.ms', imsize=[192, 192], interactive=False, niter=1000, threshold='4e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

Note that we adjusted the cleaning threshold to account for the higher noise. The results are displayed in Fig. 6. The noise is dominating the image now, which is expected.

Fig. 6: Output of simanalyze after injecting noise.

Q-band, 8GHz bandwidth, 4mm pwv, 1h integration time, A-configuration

To compensate for noise, we will extend our bandwidth from 128MHz to 8GHz. At this time, simobserve has no chanelization so the simulation will not have the better uv-coverage that is obtained in multi-frequency imaging. Consecutive simulations with different spectral windows, and combined MSs, however, can achieve that. We also need to be careful to not introduce [bandwidth smearing]. Our source is at the center of the image and bandwidth smearing will be minimal.

# In CASA
simobserve(project='psimvla4',   skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel', incenter='44GHz', inwidth='8GHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='tsys-atm', user_pwv=4, t_ground=270.0, graphics='both', overwrite=True)

The simobserve output again is identical with Fig. 3. In real, multi-channel data, however, the uv-coverage tracks would spread out and the psf would have reduced sidelobes.

# In CASA
simanalyze(project='psimvla4', image=True, vis='psimvla4.vla.a.noisy.ms',  imsize=[192, 192], interactive=False, niter=1000, threshold='2e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

The image is clearly better defined, as expected from a broader band. We also decreased the cleaning threshold to 0.2mJy.

Fig. 7: Imaging with 8GHz bandwidth.

Q-band, 8GHz bandwidth, 4mm pwv, 4h integration time, A-configuration

Let's improve the depth of the observations even further by going to 4h on-source integration, via totaltime='14400s' .

# In CASA
simobserve(project='psimvla5',   skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel',  incenter='44GHz', inwidth='8GHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='14400s',  thermalnoise='tsys-atm', user_pwv=4, t_ground=270.0, graphics='both', overwrite=True)

Fig. 8 shows a better uv-coverage (lower left image) and the 4h observation is also marked by a wider red bar in the top left image. This also improved the shape of the psf.

Fig. 8: Output of simobserve.
# In CASA
simanalyze(project='psimvla5', image=True, vis='psimvla5.vla.a.noisy.ms',  imsize=[192, 192], interactive=False, niter=1000, threshold='2e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

Given the expected deeper image we reduced the cleaning threshold to 0.1mJy, and the result is shown in Fig. 9. The image now is an almost perfect representation of the true sky model convolved with the clean beam.

Fig. 9: Images after increasing the integration time to 4h.

Q-band, 8GHz bandwidth, 18mm pwv, 1h integration time, A-configuration

We will now go back to 1h integration, but assume worse conditions of the VLA site with a precipitable water vapor of 18mm.

# In CASA
simobserve(project='psimvla6',   skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel',  incenter='44GHz', inwidth='8GHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='tsys-atm', user_pwv=18, t_ground=270.0, graphics='both', overwrite=True)
# In CASA
simanalyze(project='psimvla6', image=True, vis='psimvla6.vla.a.noisy.ms', imsize=[192, 192], interactive=False, niter=1000, threshold='2e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

We see that the resulting image is almost identical to the pwv of 4mm (Fig. 7), with almost identical noise figures. The observations seem to be very little affected by the variations in the VLA site pwv. Note that the simulations do not include pointing errors due to stronger winds.

Fig. 10: Assuming less than ideal weather conditions with pwv of 18mm.

Q-band, 8GHz bandwidth, 4mm pwv, 1h integration time, C-configuration

Our last Q-band simulation will use a precipitable water vapor of 4mm again, but to we will simulate the observations with the VLA C-configuration now.

# In CASA
simobserve(project='psimvla7',   skymodel='ppdisk672_GHz_50pc.fits', inbright='3e-5Jy/pixel',  incenter='44GHz', inwidth='8GHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.c.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='tsys-atm', user_pwv=4, t_ground=270.0, graphics='both', overwrite=True)

As expected, the uv-coverage decreases by a factor of about 10 and and the synthesize beam increases by the same factor, as seen in Fig. 11.

Fig. 11: Output of simobserve using the C-congifuration.
# In CASA
simanalyze(project='psimvla7', image=True, vis='psimvla7.vla.c.noisy.ms',  imsize=[1728, 1728], interactive=False, niter=1000, threshold='1e-4Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

Note that the output filenames will now change and reflect C-configuration in their names. The resulting images are displayed in Fig. 12, where the images now look like a point source with the degraded resolution. The surface brightness sensitivity of the observations, however, increase and we can clean to a lower threshold of 0.1mJy.

Fig. 12: C-configuration results.

X-band, 4GHz bandwidth, 4mm pwv, 1h integration time, A-configuration

Finally, we will go for a X-band, A-configuration simulation. Since X-band has a maximum bandwidth of 4GHz, we will reduce it to this value. We will also lower the flux of the model to 0.01mJy given that a the flux of the model will likely be lower, too.

# In CASA
simobserve(project='psimvla8',   skymodel='ppdisk672_GHz_50pc.fits', inbright='1e-5Jy/pixel',  incenter='10GHz', inwidth='4GHz' , setpointings=True, integration='2s',  direction='J2000 18h00m00.031s -22d59m59.6s',  mapsize= '0.78arcsec', obsmode='int', antennalist='vla.a.cfg', hourangle='transit', totaltime='3600s',  thermalnoise='tsys-atm', user_pwv=4, t_ground=270.0, graphics='both', overwrite=True)

When expressed in wavelengths, the uv-coverage and synthesized beam are somewhere in between that of Q-band A-configuration and C-configuration as expected (Fig. 13).

Fig. 13: Output of simobserve for X-band observations.
# In CASA
simanalyze(project='psimvla8', image=True, vis='psimvla8.vla.a.noisy.ms',  imsize=[192, 192], interactive=False, niter=1000, threshold='3e-5Jy', weighting='briggs',  pbcor=True, stokes='I', analyze=True, showuv=False, showpsf=True, showmodel=True, showconvolved=True, showclean=True, showresidual=True, showdifference=True, showfidelity=True, graphics='both', overwrite=True)

The ring shape is still visible in the results (Fig. 14) although quite blurred. We note that the noise is considerably lower as expected from the lower band noise properties of the receiver and atmosphere. We adjusted the clean threshold accordingly.

Fig. 14: X-band simulation results.