Skip to content

Commit

Permalink
Merge pull request achael#79 from achael/polrep
Browse files Browse the repository at this point in the history 
Polrep updates through February 1, 2019
  • Loading branch information
@achael
achael authored Feb 2, 2019
2 parents 4ba0f11 + e5d0de3 commit b359373
Show file tree
Hide file tree
Showing 9 changed files with 307 additions and 76 deletions.
1 change: 0 additions & 1 deletion README.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ Python modules for simulating and manipulating VLBI data and producing images wi

The package contains several primary classes for loading, simulating, and manipulating VLBI data. The main classes are the ``Image``, ``Array``, ``Obsdata``, ``Imager``, and ``Caltable`` classes, which provide tools for loading images and data, producing simulated data from realistic u-v tracks, calibrating, inspecting, and plotting data, and producing images from data sets in various polariazations using various data terms and regularizers.

Note that this is a pre-release of ehtim. If you have a problem please submit a pull request on the git repository.

Installation
------------
Expand Down
7 changes: 5 additions & 2 deletions ehtim/caltable.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -399,7 +399,8 @@ def pad_scans(self, maxdiff=60, padtype='median'):
return Caltable(self.ra, self.dec, self.rf, self.bw, outdict, self.tarr,
source=self.source, mjd=self.mjd, timetype=self.timetype)

def applycal(self, obs, interp='linear', extrapolate=None, force_singlepol=False, copy_closure_tables=True):
def applycal(self, obs, interp='linear', extrapolate=None,
force_singlepol=False, copy_closure_tables=True):
"""Apply the calibration table to an observation.
Args:
Expand Down Expand Up @@ -493,7 +494,9 @@ def applycal(self, obs, interp='linear', extrapolate=None, force_singlepol=False
calobs = ehtim.obsdata.Obsdata(obs.ra, obs.dec, obs.rf, obs.bw, np.array(datatable), obs.tarr,
polrep=obs.polrep, scantable=obs.scans, source=obs.source, mjd=obs.mjd,
ampcal=obs.ampcal, phasecal=obs.phasecal, opacitycal=obs.opacitycal,
dcal=obs.dcal, frcal=obs.frcal,timetype=obs.timetype,).switch_polrep(orig_polrep)
dcal=obs.dcal, frcal=obs.frcal,timetype=obs.timetype)
calobs = calobs.switch_polrep(orig_polrep)

if copy_closure_tables:
calobs.camp = obs_orig.camp
calobs.logcamp = obs_orig.logcamp
Expand Down
2 changes: 1 addition & 1 deletion ehtim/const_def.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -324,7 +324,7 @@ def prog_msg(nscan, totscans, msgtype='bar',nscan_last=0):
' ..,**,,.......,,*(((#%&&%#%&@@@@@@@@&&&&%(*... @ @ @ @ # # # # .. ',
' .,,*,,,,.......,/(((#%&&%##&@@&@@@&&%%%##(*.. @@ @ @ @@ ## # # # . ',
' ..,,,,,,,.......,/((((%&&%#(%&&&@@%%%###((/*.. . ',
' ...,,.,,,....,,,*/((/(#%&(#%&@@%%%#(///**,.. v 1.0 , 2018 . ',
' ...,,.,,,....,,,*/((/(#%&(#%&@@%%%#(///**,.. v 0.1 , 2019 . ',
' ........,.....,,*////((%&&%###&@&%###(/*,,,,... .. ',
' .............,,**///((#%&&%%#@&%%#((((/*....... ,. ',
' ............,,,**////((#%%&%@@#####((/***,.... .,. ',
Expand Down
143 changes: 107 additions & 36 deletions ehtim/image.py
View file

Large diffs are not rendered by default.

12 changes: 10 additions & 2 deletions ehtim/io/load.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,11 +218,19 @@ def load_im_fits(filename, aipscc=False, pulse=PULSE_DEFAULT,

# Update the image using the AIPS CC table
if aipscc:
# Get AIPS CC Table
try:
aipscctab = hdulist["AIPS CC"]
except:
raise ValueError("Input FITS file does not have an AIPS CC table.")
print("Input FITS file does not have an AIPS CC table. Loading image instead.")
aipscc = False
#raise ValueError("Input FITS file does not have an AIPS CC table.")

if aipscc:
# Get AIPS CC Table
# try:
# aipscctab = hdulist["AIPS CC"]
# except:
# raise ValueError("Input FITS file does not have an AIPS CC table.")

print("loading the AIPS CC table.")
print("force the pulse function to be the delta function.")
Expand Down
18 changes: 15 additions & 3 deletions ehtim/obsdata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1251,11 +1251,13 @@ def cleanbeam(self, npix, fov, pulse=PULSE_DEFAULT):

return im

def fit_beam(self):
def fit_beam(self, weighting='uniform', units='rad'):

"""Fit a Gaussian to the dirty beam and return the parameters (fwhm_maj, fwhm_min, theta).
Args:
weighting (str): 'uniform' or 'natural'.
units (string): 'rad' returns values in radians, 'natural' returns FWHMs in uas and theta in degrees
Returns:
(tuple): a tuple (fwhm_maj, fwhm_min, theta) of the dirty beam parameters in radians.
Expand All @@ -1279,8 +1281,13 @@ def fit_chisq(beamparams, db_coeff):
# For a point (x,y) in the image plane, the dirty beam expansion is 1-ax^2-by^2-cxy
u = self.unpack('u')['u']
v = self.unpack('v')['v']
sigma = self.unpack('sigma')['sigma']
n = float(len(u))
abc = (2.*np.pi**2/n) * np.array([np.sum(u**2), np.sum(v**2), 2*np.sum(u*v)])
weights = np.ones(u.shape)
if weighting == 'natural':
weights = 1./sigma**2

abc = (2.*np.pi**2/np.sum(weights)) * np.array([np.sum(weights*u**2), np.sum(weights*v**2), 2*np.sum(weights*u*v)])
abc = 1e-20 * abc # Decrease size of coefficients

# Fit the beam
Expand All @@ -1299,6 +1306,11 @@ def fit_chisq(beamparams, db_coeff):

gparams = np.array((fwhm_maj, fwhm_min, theta))

if units == 'natural':
gparams[0] /= RADPERUAS
gparams[1] /= RADPERUAS
gparams[2] *= 180./np.pi

return gparams


Expand Down Expand Up @@ -1340,7 +1352,7 @@ def dirtybeam(self, npix, fov, pulse=PULSE_DEFAULT, weighting='uniform'):
# TODO right normalization?

src = self.source + "_DB"
outim2 = ehtim.image.Image(im, pdim, self.ra, self.dec, rf=self.rf, source=src, mjd=self.mjd, pulse=pulse)
outim = ehtim.image.Image(im, pdim, self.ra, self.dec, rf=self.rf, source=src, mjd=self.mjd, pulse=pulse)

return outim

Expand Down
54 changes: 37 additions & 17 deletions ehtim/observing/obs_simulate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -973,13 +973,33 @@ def add_noise(obs, add_th_noise=True, opacitycal=True, ampcal=True, phasecal=Tru
uv = recarr_to_ndarr(obsdata[['u','v']],'f8')
taus = np.abs(recarr_to_ndarr(obsdata[['tau1','tau2']],'f8'))
elevs = recarr_to_ndarr(obs.unpack(['el1','el2'],ang_unit='deg'),'f8')
time = obsdata['time']
times = obsdata['time']
tint = obsdata['tint']
vis1 = obsdata[obs.poldict['vis1']]
vis2 = obsdata[obs.poldict['vis2']]
vis3 = obsdata[obs.poldict['vis3']]
vis4 = obsdata[obs.poldict['vis4']]

times_stable_phase = times.copy()
times_stable_amp = times.copy()
times_stable = times.copy()

if stabilize_scan_phase==True or stabilize_scan_amp==True:
scans = obs.scans
if np.all(scans) == None or len(scans) == 0:
print('ADDING SCANS')
obs_scans = obs.copy()
obs_scans.add_scans()
scans = obs_scans.scans
for j in range(len(times_stable)):
for scan in scans:
if scan[0] <= times_stable[j] and scan[1] >= times_stable[j]:
times_stable[j] = scan[0]
break

if stabilize_scan_phase==True: times_stable_phase = times_stable.copy()
if stabilize_scan_amp==True: times_stable_amp = times_stable.copy()

# Recompute perfect sigmas from SEFDs
bw = obs.bw
if np.any(obs.tarr['sefdr'] <= 0):
Expand Down Expand Up @@ -1018,26 +1038,26 @@ def add_noise(obs, add_th_noise=True, opacitycal=True, ampcal=True, phasecal=Tru
if not ampcal:
# Amplitude gain
if type(gain_offset) == dict:
goff1 = np.fromiter((gain_offset[sites[i,0]] for i in range(len(time))), float)
goff2 = np.fromiter((gain_offset[sites[i,1]] for i in range(len(time))), float)
goff1 = np.fromiter((gain_offset[sites[i,0]] for i in range(len(times))), float)
goff2 = np.fromiter((gain_offset[sites[i,1]] for i in range(len(times))), float)
else:
goff1 = np.fromiter((gain_offset for i in range(len(time))), float)
goff2 = np.fromiter((gain_offset for i in range(len(time))), float)
goff1 = np.fromiter((gain_offset for i in range(len(times))), float)
goff2 = np.fromiter((gain_offset for i in range(len(times))), float)

if type(gainp) is dict:
gain_mult_1 = np.fromiter((gainp[sites[i,0]] for i in range(len(time))), float)
gain_mult_2 = np.fromiter((gainp[sites[i,1]] for i in range(len(time))), float)
gain_mult_1 = np.fromiter((gainp[sites[i,0]] for i in range(len(times))), float)
gain_mult_2 = np.fromiter((gainp[sites[i,1]] for i in range(len(times))), float)
else:
gain_mult_1 = np.fromiter((gainp for i in range(len(time))), float)
gain_mult_2 = np.fromiter((gainp for i in range(len(time))), float)
gain_mult_1 = np.fromiter((gainp for i in range(len(times))), float)
gain_mult_2 = np.fromiter((gainp for i in range(len(times))), float)

gain1 = np.abs(np.fromiter(
((1.0 + goff1[i] * np.abs(hashrandn(sites[i,0], 'gain', seed)))*(1.0 + gain_mult_1[i] * hashrandn(sites[i,0], 'gain', time[i], seed))
for i in range(len(time))
((1.0 + goff1[i] * np.abs(hashrandn(sites[i,0], 'gain', seed)))*(1.0 + gain_mult_1[i] * hashrandn(sites[i,0], 'gain', times_stable_amp[i], seed))
for i in range(len(times))
), float))
gain2 = np.abs(np.fromiter(
((1.0 + goff2[i] * np.abs(hashrandn(sites[i,1], 'gain', seed)))*(1.0 + gain_mult_2[i] * hashrandn(sites[i,1], 'gain', time[i], seed))
for i in range(len(time))
((1.0 + goff2[i] * np.abs(hashrandn(sites[i,1], 'gain', seed)))*(1.0 + gain_mult_2[i] * hashrandn(sites[i,1], 'gain', times_stable_amp[i], seed))
for i in range(len(times))
), float))
gain_true = np.sqrt(gain1 * gain2)
else:
Expand All @@ -1049,8 +1069,8 @@ def add_noise(obs, add_th_noise=True, opacitycal=True, ampcal=True, phasecal=Tru
tau_est = np.sqrt(np.exp(taus[:,0]/(EP+np.sin(elevs[:,0]*DEGREE)) + taus[:,1]/(EP+np.sin(elevs[:,1]*DEGREE))))

# Opacity Errors
tau1 = np.abs(np.fromiter((taus[i,0]* (1.0 + taup * hashrandn(sites[i,0], 'tau', time[i], seed)) for i in range(len(time))),float))
tau2 = np.abs(np.fromiter((taus[i,1]* (1.0 + taup * hashrandn(sites[i,1], 'tau', time[i], seed)) for i in range(len(time))),float))
tau1 = np.abs(np.fromiter((taus[i,0]* (1.0 + taup * hashrandn(sites[i,0], 'tau', times_stable_amp[i], seed)) for i in range(len(times))),float))
tau2 = np.abs(np.fromiter((taus[i,1]* (1.0 + taup * hashrandn(sites[i,1], 'tau', times_stable_amp[i], seed)) for i in range(len(times))),float))

# Correct noise RMS for opacity
tau_true = np.sqrt(np.exp(tau1/(EP+np.sin(elevs[:,0]*DEGREE)) + tau2/(EP+np.sin(elevs[:,1]*DEGREE))))
Expand Down Expand Up @@ -1082,8 +1102,8 @@ def add_noise(obs, add_th_noise=True, opacitycal=True, ampcal=True, phasecal=Tru

# Add random atmospheric phases
if not phasecal:
phase1 = np.fromiter((2 * np.pi * hashrand(sites[i,0], 'phase', time[i], seed) for i in range(len(time))),float)
phase2 = np.fromiter((2 * np.pi * hashrand(sites[i,1], 'phase', time[i], seed) for i in range(len(time))),float)
phase1 = np.fromiter((2 * np.pi * hashrand(sites[i,0], 'phase', times_stable_phase[i], seed) for i in range(len(times))),float)
phase2 = np.fromiter((2 * np.pi * hashrand(sites[i,1], 'phase', times_stable_phase[i], seed) for i in range(len(times))),float)

vis1 *= np.exp(1j * (phase2-phase1))
vis2 *= np.exp(1j * (phase2-phase1))
Expand Down
5 changes: 3 additions & 2 deletions ehtim/plotting/comp_plots.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@

from builtins import range
import numpy as np
import numpy.matlib as matlib
import matplotlib as mpl
import matplotlib.pyplot as plt
import itertools as it
Expand Down Expand Up @@ -214,7 +215,7 @@ def plot_cphase_compare(obslist, imlist, site1, site2, site3,
except TypeError: obslist = [obslist]

if len(cphases)==0:
cphases = np.matlib.repmat([],len(obslist),1)
cphases = matlib.repmat([],len(obslist),1)

if len(cphases) != len(obslist):
raise Exception("cphases list must be same length as obslist!")
Expand Down Expand Up @@ -305,7 +306,7 @@ def plot_camp_compare(obslist, imlist, site1, site2, site3, site4,
except TypeError: obslist = [obslist]

if len(camps)==0:
camps = np.matlib.repmat([],len(obslist),1)
camps = matlib.repmat([],len(obslist),1)

if len(camps) != len(obslist):
raise Exception("camps list must be same length as obslist!")
Expand Down
Loading

0 comments on commit b359373

Please sign in to comment.