How to use the ehtim.RADPERUAS function in ehtim

tstop_hr = 24
bw_hz = 4e9
obs = mod.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, ampcal=True, phasecal=False)

# Fit the observation using visibility amplitudes and closure phase
mod_init  = eh.model.Model()
mod_init  = mod_init.add_mring(1.5, 40.*eh.RADPERUAS, beta_list=[0.1+0.1j]).add_circ_gauss(1., 20.*eh.RADPERUAS)
mod_prior = mod_init.default_prior()
mod_prior[0]['x0'] = {'prior_type':'fixed'}
mod_prior[0]['y0'] = {'prior_type':'fixed'}
mod_prior[0]['F0'] = {'prior_type':'fixed'}
mod_prior[1]['F0'] = {'prior_type':'gauss','mean':1.0,'std':0.1}
mod_fit   = eh.modeler_func(obs, mod_init, mod_prior, d1='amp', d2='cphase')

eh.comp_plots.plotall_obs_im_compare(obs,mod_fit,'uvdist','amp') 
mod_fit.blur_circ(5.*eh.RADPERUAS).display()

# 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

import ehtim.modeling.modeling_utils as mu

# Define the ground-truth model
mod = eh.model.Model()
mod.add_ring(1., 50.*eh.RADPERUAS)
mod.add_ring(0.5, 30.*eh.RADPERUAS, 5.*eh.RADPERUAS, 5.*eh.RADPERUAS) #For testing gradients
mod.make_image(100.*eh.RADPERUAS, 128).blur_circ(5.*eh.RADPERUAS).display()

# Create an observation
obs = eh.obsdata.load_uvfits('/home/michael/Dropbox/ER5_polarization_calibration/ER5/postproc-hops-lo/3.+netcal/3601/hops_3601_M87+netcal.uvfits')
obs.add_scans()
obs = obs.avg_coherent(0.,scan_avg=True)
obs = mod.observe_same(obs,ampcal=True,phasecal=True)

# Testing the chi^2
dtypes = ['vis','amp','bs','cphase','logcamp','camp','logamp', 'logcamp_diag', 'cphase_diag'] #, 'bs', 'amp', 'cphase', 'cphase_diag', 'camp', 'logcamp', 'logcamp_diag']
for dtype in dtypes:
    print('\nTesting chi^2 dtype:',dtype)
    chisqdata = eh.modeling.modeling_utils.chisqdata(obs, dtype)
    chisq = eh.modeling.modeling_utils.chisq(mod, chisqdata[2], chisqdata[0], chisqdata[1], dtype)
    print("chisq: %f" % chisq)
    print('\nTesting gradient')
    for x in [['F0',1e-10,0],['d',1e-4*eh.RADPERUAS,1],['x0',1e-4*eh.RADPERUAS,2],['y0',1e-4*eh.RADPERUAS,3]]:

# This examples shows an image reconstruction using stochastic optics scattering mitigation at 86 GHz.
# 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

# 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

mod_fit = eh.modeler_func(obs, mod_fit, prior, d1='vis',alpha_d1=1,minimizer_method='L-BFGS-B')

# Fit holding some parameters fixed
mod_init = eh.model.Model()
mod_init.add_ring(1.5, 35.*eh.RADPERUAS)
mod_init.add_ring(0.5, 5.*eh.RADPERUAS)
prior = mod_init.default_prior()
prior[0]['x0'] = {'prior_type':'fixed'}
prior[0]['y0'] = {'prior_type':'fixed'}
mod_fit = eh.modeler_func(obs, mod_init, prior, d1='amp',d2='cphase',alpha_d1=1,stop=1e-50,minimizer_method='Powell')
mod_fit.make_image(100.*eh.RADPERUAS, 128).blur_circ(5.*eh.RADPERUAS).display()

# Fit a multi-component model
mod_init = eh.model.Model()
mod_init.add_circ_gauss(0.2, 5.*eh.RADPERUAS, 10.*eh.RADPERUAS, 0)
mod_init.add_circ_gauss(0.2, 5.*eh.RADPERUAS, 0., 10.*eh.RADPERUAS)
prior = mod_init.default_prior()
for j in range(len(prior)):
    prior[j]['F0']['prior_type'] = 'fixed'

mod_fit2 = eh.modeler_func(obs, mod_init, prior, d1='vis',alpha_d1=1,stop=1e-50,maxit=1000)

import ehtim.image as image
from ehtim.imaging import starwarps as sw
import matplotlib.pyplot as plt
import sys, os, copy
import scipy
import scipy.optimize as opt


############## parameters ##############

# data file
obsname = 'sample_movie.uvfits'

# image parameters
flux = 2.0
fwhm = 50 * eh.RADPERUAS
fov = 100 * eh.RADPERUAS
NPIX = 30
npixels = NPIX**2

# StarWarps optimization parameters
warp_method = 'phase'
measurement = {'vis':1 } # {'amp':1, 'cphase':1}
interiorPriors = True
reassign_apxImgs = False
numLinIters = 5
variance_img_diff = 1e-7

# parameters associated with EM 
nIters = 30
NHIST = 5000
stop=1e-10

# Add non-closing systematic noise to the observation
obs_sc = obs_sc.add_fractional_noise(sys_noise)

# Make a copy of the initial data (before any self-calibration but after the taper)
obs_sc_init = obs_sc.copy()

# Self-calibrate the LMT to a Gaussian model
# (Refer to Section 4's "Pre-Imaging Considerations")
print("Self-calibrating the LMT to a Gaussian model for LMT-SMT...")

obs_LMT = obs_sc_init.flag_uvdist(uv_max=2e9) # only consider the short baselines (LMT-SMT)
if reverse_taper_uas > 0:
    # start with original data that had no reverse taper applied.
    # Re-taper, if necessary
    obs_LMT = obs_LMT.taper(reverse_taper_uas*eh.RADPERUAS)

# Make a Gaussian image that would result in the LMT-SMT baseline visibility amplitude
# as estimated in Section 4's "Pre-Imaging Considerations".
# This is achieved with a Gaussian of size 60 microarcseconds and total flux of 0.6 Jy
gausspriorLMT = eh.image.make_square(obs, npix, fov)
gausspriorLMT = gausspriorLMT.add_gauss(0.6, (60.0*eh.RADPERUAS, 60.0*eh.RADPERUAS, 0, 0, 0))

# Self-calibrate the LMT visibilities to the gausspriorLMT image
# to enforce the estimated LMT-SMT visibility amplitude
caltab = eh.selfcal(obs_LMT, gausspriorLMT, sites=['LM'], gain_tol=1.0,
                    method='both', ttype=ttype, caltable=True)

# Spply the calibration solution to the full (and potentially tapered) dataset
obs_sc = caltab.applycal(obs_sc, interp='nearest', extrapolate=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)

# Export the visibility data to uvfits/text
#obs.save_txt('obs.txt') # exports a text file with the visibilities
#obs.save_uvfits('obs.uvp') # exports a UVFITS file modeled on template.UVP

# Generate an image prior
npix = 256
fov = 1*im.fovx()
zbl = im.total_flux() # total flux
prior_fwhm = 200*eh.RADPERUAS # Gaussian size in microarcssec
emptyprior = eh.image.make_square(obs, npix, fov)
flatprior = emptyprior.add_flat(zbl)
gaussprior = emptyprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0))

# Image total flux with bispectrum
flux = zbl
out  = eh.imager_func(obs, gaussprior, gaussprior, flux,
                      d1='bs', s1='simple',
                      alpha_s1=1, alpha_d1=100,
                      alpha_flux=100, alpha_cm=50,
                      maxit=100, ttype='nfft')

# Blur the image with a circular beam and image again to help convergance
out = out.blur_circ(res)
out = eh.imager_func(obs, out, out, flux,
                d1='bs', s1='tv',

def padNewFOV(im, fov_arcseconds):
    
    oldfov = im.psize * im.xdim
    newfov = fov_arcseconds * ehtim.RADPERUAS
    tnpixels = np.ceil(im.xdim * newfov/oldfov).astype('int')
    
    origimg = np.reshape(im.imvec, [im.xdim, im.xdim])
    padimg = np.pad(origimg, ((0,tnpixels-im.xdim), (0,tnpixels-im.xdim)), 'constant')
    padimg = np.roll(padimg, np.floor((tnpixels-im.xdim)/2.).astype('int'), axis=0)
    padimg = np.roll(padimg, np.floor((tnpixels-im.xdim)/2.).astype('int'), axis=1)
    
    return image.Image(padimg.reshape((tnpixels, tnpixels)), im.psize, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd, pulse=im.pulse)

ax.set_xlim(0, 22)
    ax.set_ylim(-10, 22)

    ax.set_xticks(x_loc)
    ax.set_xticklabels(fracsteps)


    ax.set_yticks([])
    ax.set_yticklabels([])

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.spines['bottom'].set_visible(False)

    ax.set_title('Blurred comparison of all images; cutoff={0}, fov (uas)={1}'.format(str(cutoff), str(im_clique_fraclevels[-1][-1].fovx()/eh.RADPERUAS)))

#     for item in [fig, ax]:
#         item.patch.set_visible(False)
#     fig.patch.set_visible(False)
#     ax.axis('off')
    if show == True:
        plt.show()

# 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)

# Export the visibility data to uvfits/text
#obs.save_txt('obs.txt') # exports a text file with the visibilities
#obs.save_uvfits('obs.uvp') # exports a UVFITS file modeled on template.UVP

# Generate an image prior
npix = 256
fov = 1*im.fovx()
zbl = im.total_flux() # total flux
prior_fwhm = 200*eh.RADPERUAS # Gaussian size in microarcssec
emptyprior = eh.image.make_square(obs, npix, fov)
flatprior = emptyprior.add_flat(zbl)
gaussprior = emptyprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0))

# Average the closure quantities and add them to the obsdata object

avg_time = 600
obs.add_bispec(avg_time=avg_time)
obs.add_amp(avg_time=avg_time)
obs.add_cphase(avg_time=avg_time)
obs.add_camp(avg_time=avg_time)
obs.add_logcamp(avg_time=avg_time)

# Image directly with these averaged data
flux = zbl
out  = eh.imager_func(obs, gaussprior, gaussprior, flux,