How to use the ehtim.image.make_square function in ehtim

def cleanbeam(self, npix, fov, pulse=PULSE_DEFAULT):

        """Make an image of the observation clean beam.

           Args:
               npix (int): The pixel size of the square output image.
               fov (float): The field of view of the square output image in radians.
               pulse (function): The function convolved with the pixel values for continuous image.

           Returns:
               (Image): an Image object with the clean beam.
        """

        im = ehtim.image.make_square(self, npix, fov, pulse=pulse)
        beamparams = self.fit_beam()
        im = im.add_gauss(1.0, beamparams)
        return im

# 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',
                alpha_s1=1, alpha_d1=50,

sgrscat=False, ampcal=True, phasecal=True)
obslist_nospec = [obs230,obs345_nospec]

# Resolution
beamparams230 = obs230.fit_beam() # fitted beam parameters (fwhm_maj, fwhm_min, theta) in radians
res230 = obs230.res() # nominal array resolution, 1/longest baseline
beamparams345 = obs345.fit_beam() # fitted beam parameters (fwhm_maj, fwhm_min, theta) in radians
res345 = obs345.res() # nominal array resolution, 1/longest baseline
print("Nominal Resolution: " ,res230,res345)

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


## Image both frequencies independently with complex visibilities
flux230 = zbl
imgr230  = eh.imager.Imager(obs230, gaussprior, gaussprior, flux230,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=200, ttype='nfft')
imgr230.make_image_I(show_updates=True)
out230 = imgr230.out_last()


flux345 = im345.total_flux()
imgr345  = eh.imager.Imager(obs345, gaussprior, gaussprior, flux345,

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

outbs = dd_clean_amp_cphase(obs2, prior, niter=50, loop_gain=0.1, loop_gain_init=.01,phaseweight=2, weighting='uniform', bscount="min")

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)

#-------------------------------------------------------------------------------
# Reconstruct an image
#-------------------------------------------------------------------------------


# First  Round of Imaging

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
npix = 128
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
imgr  = eh.imager.Imager(obs, gaussprior, gaussprior, flux,
                          data_term={'bs':100}, show_updates=False,
                          reg_term={'simple':1,'flux':100,'cm':50},
                          maxit=200, ttype='nfft')
imgr.make_image_I()

# Blur the image with a circular beam and image again to help convergance
out = imgr.out_last()
imgr.init_next = out.blur_circ(res)
imgr.prior_next = imgr.init_next
imgr.dat_term_next = {'amp':50, 'cphase':100}

# 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,
                      d1='bs', s1='simple',