Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# Note: must import ehtim outside the ehtim directory
# either in parent eht-imaging directory or after installing with setuptools
import numpy as np
import ehtim as eh
import ehtim.scattering as so
import ehtim.plotting as pl
# Load the image and the array
im = eh.image.load_txt('../models/avery_sgra_eofn.txt')
eht = eh.array.load_txt('../arrays/EHT2017.txt')
# If the image has an even number of pixels, make it odd
#if im.xdim%2 == 0:
# newim = im.imvec.reshape(im.ydim, im.xdim)
# newim = newim[:-1,:-1]
# im = eh.image.Image(newim, im.psize, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
# Create a ScatteringModel; this defaults to parameters appropriate for Sgr A*
sm = so.ScatteringModel()
# Scatter the image; display the result
im_scatt = sm.Scatter(im,DisplayImage=True)
# Compare visibilities of the unscattered, scattered, and deblurred images
def LoadSeriesImages(sourceName):
#loads a list of images saved according to the BU Blazar library convention
dirName = sourceName+'_READ_CLEAN_files'
fileList = os.listdir(dirName)
filenameList = [dirName+'/'+filename for filename in fileList]
im_List = [image.load_txt(filename) for filename in filenameList]
return im_List
#example script for running data domain clean
import ehtim as eh
from ehtim.imaging.clean import *
#################################################################
# Data domain clean with complex visibilities
im = eh.image.load_txt('./models/avery_sgra_eofn.txt')
arr = eh.array.load_txt('./arrays/EHT2017.txt')
#arr = eh.array.load_txt('./arrays/EHT2025.txt')
obs = im.observe(arr, 1000, 600, 0, 24., 4.e10, add_th_noise=False, phasecal=True)
prior = eh.image.make_square(obs, 128, 1.5*im.fovx())
# data domain clean with visibilities
outvis = dd_clean_vis(obs, prior, niter=100, loop_gain=0.1, method='min_chisq',weighting='uniform')
#################################################################
# Data domain clean directly with bispectrum
# trial image 2 -- 2 Gaussians
im2 = eh.image.make_square(obs, 256, 3*im.fovx())
im2 = im2.add_gauss(1., (1*eh.RADPERUAS, 1*eh.RADPERUAS, 0, 0, 0))
im2 = im2.add_gauss(.7, (1*eh.RADPERUAS, 1*eh.RADPERUAS, 0, -75*eh.RADPERUAS, 30*eh.RADPERUAS))
obs2 = im2.observe(arr, 600, 600, 0, 24., 4.e9, add_th_noise=False, phasecal=False)
prior = eh.image.make_square(obs, 50, 3*im.fovx())
# 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
# 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
# This is a rough script to verify the consistency of Fourier transform types and to check gradients of the various image regularization options
from __future__ import division
from __future__ import print_function
import numpy as np
import ehtim as eh
im = eh.image.load_txt('../models/avery_sgra_eofn.txt')
eht = eh.array.load_txt('../arrays/EHT2017.txt')
tint_sec = 5
tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 4e9
obs_dft = im.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True, ttype='direct', add_th_noise=False)
obs_nfft = im.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True, ttype='nfft', add_th_noise=False)
prior = im.copy()
im2 = im.copy() # This is our test image
# Add random noise to the image
for j in range(len(im2.imvec)):
im2.imvec[j] *= (1.0 + (np.random.rand()-0.5)/10.0)
# 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)