Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Create a sample unscattered image -- a ring
npix = 200 #Linear pixel dimension
fov = 500*eh.RADPERUAS #Field of view
total_flux = 3.4 # Jy
im = eh.image.Image(np.zeros((npix,npix)), fov/npix, eh.RA_DEFAULT, eh.DEC_DEFAULT, rf=86e9)
im = im.add_crescent(total_flux, 105*eh.RADPERUAS, 95*eh.RADPERUAS, 0, 0)
im = im.blur_circ(10*eh.RADPERUAS)
im.display()
# Scatter the image
ep = so.MakeEpsilonScreen(im.xdim,im.ydim,rngseed=34)
sm = so.ScatteringModel()
scatt = sm.Scatter(im,ep,DisplayImage=True)
# Observe the average image
eht = eh.array.load_txt('../arrays/VLBA_86GHz.txt') # Load the observing array
#observing parameters:
tint_sec = 60
tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 0.5e9
#create the observation
obs = scatt.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True)
# Generate an image prior
npix = 55 #This must be odd
prior_fwhm = 300*eh.RADPERUAS # Gaussian size in microarcssec
gaussprior = eh.image.Image(np.zeros((npix,npix)), fov/npix, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
gaussprior = gaussprior.add_gauss(total_flux,(prior_fwhm, prior_fwhm, 0, 0, 0))
# Now try imaging
# Note: must import ehtim outside the ehtim directory
# either in parent eht-imaging directory or after installing with setuptools
from __future__ import division
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import ehtim as eh
from ehtim.calibrating import self_cal as sc
#from ehtim.plotting import self_cal as sc
# Load the image and the array
im = eh.image.load_txt('../models/avery_sgra_eofn.txt')
eht = eh.array.load_txt('../arrays/EHT2017.txt')
# Look at the image
im.display()
# Observe the image
# tint_sec is the integration time in seconds, and tadv_sec is the advance time between scans
# tstart_hr is the GMST time of the start of the observation and tstop_hr is the GMST time of the end
# bw_hz is the bandwidth in Hz
# sgrscat=True blurs the visibilities with the Sgr A* scattering kernel for the appropriate image frequency
# ampcal and phasecal determine if gain variations and phase errors are included
tint_sec = 5
tadv_sec = 30
tstart_hr = 0
tstop_hr = 24
bw_hz = 4e9
obs = im.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz,
# Note: this is an example sequence of commands to run in ipython
# The matplotlib windows may not open/close properly if you run this directly as a script
from __future__ import division
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import ehtim as eh
from ehtim.calibrating import self_cal as sc
#from ehtim.plotting import self_cal as sc
# Load the image and the array
im = eh.image.load_txt('../models/avery_sgra_eofn.txt')
eht = eh.array.load_txt('../arrays/EHT2017.txt')
# Observe the image
tint_sec = 5
tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 400e9
obs = im.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz,
sgrscat=False, ampcal=True, phasecal=True)
# Resolution
beamparams = obs.fit_beam() # fitted beam parameters (fwhm_maj, fwhm_min, theta) in radians
res = obs.res() # nominal array resolution, 1/longest baseline
print("Clean beam parameters: " , beamparams)
print("Nominal Resolution: " ,res)
# Generate an image prior
# Multifrequency imager example
# RIAF model with a constant spectral index
from __future__ import division
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
import ehtim as eh
from ehtim.calibrating import self_cal as sc
#from ehtim.plotting import self_cal as sc
# Load the image and the array
im = eh.image.load_txt('./models/avery_sgra_eofn.txt')
eht = eh.array.load_txt('./arrays/EHT2017.txt')
# Add a spectral index to the image
alpha = 2.5
im.rf = 230.e9
im = im.add_const_mf(alpha)
im230 = im.get_image_mf(230.e9)
im345 = im.get_image_mf(345.e9)
# Observe the image at two different frequencies
tint_sec = 120
tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 2e9
obs230 = im230.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz,