In [1]:
import numpy as np
import skimage.transform
import astropy.io.fits as fits
import matplotlib.pyplot as plt
import os
import matplotlib.patches
import PIL.ImageDraw
import PIL.Image
%pylab inline --no-import-all
matplotlib.rcParams['image.origin'] = 'lower'
matplotlib.rcParams['image.interpolation'] = 'nearest'
matplotlib.rcParams['image.cmap'] = 'gray'
In [2]:
prim_key_map = {'hex1':'hex1', 'hex2':'hex2', 'hex3':'hex3', 'hex4':'hex4',
'key24':'keystone24', 'pie08':'piewedge8', 'pie12':'piewedge12'}
secobs_key_map = {'Cross':'cross', 'X':'x'}
In [151]:
overwrite = False
N = 125 # pupil array quadrant width after binning
iD = 20
oD = 82
aligntol = 3 # units of thousandths of pupil diameter
pad = 5 # units of thousandths of pupil diameter
prim_key = "hex3"
centobs = True
secobs_key = "X"
gap = 2 # minimum segment gap width, in units of 0.2% diameter, or pixels at N=250
D = 1000
symm = 'quart' # set to either 'quart' or 'half'
shape = 'ann' # set to 'ann' or 'hex'
In [152]:
#LS_dir = os.path.abspath("../InputMasks/LS")
#LS_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks/LS")
#LS_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks_v2/LS")
#LS_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks_v3/LS")
LS_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks_v4/LS")
#LS_dir = os.path.normpath("/Users/neil/Box Sync/scda/InputMasks/LS")
if not os.path.exists(LS_dir):
os.mkdir(LS_dir)
print("created {:s} for LS arrays".format(LS_dir))
else:
print("Destination {:s} already exists".format(LS_dir))
In [153]:
#telap_dir = os.path.abspath('../Apertures/JPL/offset_masks')
telap_dir = os.path.normpath('/astro/opticslab1/SCDA/Apertures/JPL/offset_masks')
opt_telap_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks_v4/TelAp")
prim_fname = os.path.join(telap_dir, "{0:s}_{1:04d}pix_offset.fits".format(prim_key_map[prim_key], D))
secobs_fname = os.path.join(telap_dir, "{0:s}_spiders_{1:04d}pix_2.5cm_offset.fits".format(
secobs_key_map[secobs_key], D))
prim_hdulist = fits.open(prim_fname, "readonly")
prim = prim_hdulist[0].data
prim_hdulist.close()
secobs_hdulist = fits.open(secobs_fname, "readonly")
secobs = secobs_hdulist[0].data
secobs_hdulist.close()
telap = prim*secobs
In [154]:
plt.figure(figsize=(10,10))
plt.imshow(telap[500:700,500:700])
Out[154]:
In [155]:
max_shift = int(round(D*float(pad)/1000))
shift_range = range(-max_shift,max_shift+1,1)
In [156]:
[Xshifts, Yshifts] = np.meshgrid(shift_range, shift_range)
In [157]:
allowed_shifts = np.less_equal(Xshifts**2 + Yshifts**2, max_shift**2)
In [158]:
allowed_shifts
Out[158]:
In [159]:
allowed_shifts.shape
Out[159]:
In [160]:
XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])
In [161]:
len(XYshifts_allowed)
Out[161]:
In [162]:
padded_telap_accum = np.ones(telap.shape)
for (xshift,yshift) in XYshifts_allowed:
telap_shifted = np.roll(np.roll(telap, yshift, 0), xshift, 1)
padded_telap_accum *= telap_shifted
In [163]:
plt.figure(figsize=(10,10))
plt.imshow(padded_telap_accum)
Out[163]:
In [164]:
L = secobs.shape[0]
inD = float(iD)/100
outD = float(oD)/100
xs = np.linspace(-L/2 + 0.5, L/2 - 0.5, L)
[Xs, Ys] = np.meshgrid(xs, xs)
if shape is 'hex':
xycent = (L/2-0.5,L/2-0.5)
if prim_key is 'hex1':
hexagon_obj = matplotlib.patches.RegularPolygon(xycent,6,radius=outD*D/2/np.cos(np.pi/6), orientation=np.pi/6)
else:
hexagon_obj = matplotlib.patches.RegularPolygon(xycent,6,radius=outD*D/2/np.cos(np.pi/6))
# hexagon_obj = matplotlib.patches.RegularPolygon(xycent,6,radius=outD*D/2, orientation=np.pi/6)
hexagon_verts = hexagon_obj.get_verts()
hexagon_img = PIL.Image.new('L', secobs.shape, 0)
PIL.ImageDraw.Draw(hexagon_img).polygon(list(hexagon_verts.ravel()), outline=1, fill=1)
hexagon_mask = np.array(hexagon_img)
Lyot_stop = np.logical_and(np.round(padded_telap_accum), hexagon_mask)
inside_ann = np.less_equal(Xs**2 + Ys**2, (inD*D/2)**2)
Lyot_stop[inside_ann] = 0.
else:
inside_ann = np.less_equal(Xs**2 + Ys**2, (inD*D/2)**2)
outside_ann = np.greater_equal(Xs**2 + Ys**2, (outD*D/2)**2)
Lyot_stop = np.round(padded_telap_accum)
Lyot_stop[inside_ann] = 0.
Lyot_stop[outside_ann] = 0.
In [165]:
plt.figure(figsize=(12,12))
plt.imshow(Lyot_stop)
Out[165]:
In [166]:
#scalefac = float(N)/(D/2)
#Lyot_stop_binned = scipy.ndimage.zoom(Lyot_stop, scalefac, order=1)
N_orig = D/2
scalefac = int(N_orig/N)
print("Binning the original LS array {0:d}x".format(scalefac))
Lyot_stop_binned = np.reshape(Lyot_stop, (Lyot_stop.shape[0]/scalefac, scalefac,
Lyot_stop.shape[1]/scalefac, scalefac)).mean(1).mean(2)
telap_binned = np.reshape(telap, (telap.shape[0]/scalefac, scalefac,
telap.shape[1]/scalefac, scalefac)).mean(1).mean(2)
In [167]:
L_bin = Lyot_stop_binned.shape[0]
if symm is 'half':
Lyot_stop_binned_crop = Lyot_stop_binned[L_bin/2-N:L_bin/2+N,L_bin/2:L_bin/2+N]
Lyot_stop_binned_crop_binary = np.round(Lyot_stop_binned_half).astype(int)
print Lyot_stop_binned_crop.shape
else:
Lyot_stop_binned_crop = Lyot_stop_binned[L_bin/2:L_bin/2+N,L_bin/2:L_bin/2+N]
Lyot_stop_binned_crop_binary = np.round(Lyot_stop_binned_crop).astype(int)
print Lyot_stop_binned_crop.shape
Lyot_stop_full_binary = np.round(Lyot_stop_binned[L_bin/2-N:L_bin/2+N,L_bin/2-N:L_bin/2+N]).astype(int)
In [168]:
if symm is 'half':
plt.figure(figsize=(20,10))
plt.imshow(Lyot_stop_binned_crop_binary)
#plt.imshow(Lyot_stop_binned_crop)
else:
plt.figure(figsize=(10,10))
#plt.imshow(Lyot_stop_binned_crop)
plt.imshow(Lyot_stop_binned_crop_binary)
In [169]:
if symm is 'half':
LS_dat_fname_tail = "LS_half_{0:s}{1:02d}D{2:02d}_{3:s}{4:s}025cobs1Pad{5:02d}_N{6:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), prim_key, secobs_key, pad, N)
else:
LS_dat_fname_tail = "LS_quart_{0:s}{1:02d}D{2:02d}_{3:s}{4:s}025cobs1Pad{5:02d}_N{6:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), prim_key, secobs_key, pad, N)
LS_dat_fname = os.path.join(LS_dir, LS_dat_fname_tail)
if not os.path.exists(LS_dat_fname) or overwrite==True:
#np.savetxt(LS_dat_fname, Lyot_stop_binned_crop, fmt='%.6f', delimiter=" ")
np.savetxt(LS_dat_fname, Lyot_stop_binned_crop_binary, fmt='%d', delimiter=" ")
print("Wrote binned, cropped LS array to {0:s}".format(LS_dat_fname))
else:
print("LS array {0:s} already exists, will not overwrite".format(LS_dat_fname))
In [170]:
os.listdir(LS_dir)
Out[170]:
In [171]:
# os.remove(LS_dat_fname)
In [172]:
opt_telap_dir = os.path.normpath("/astro/opticslab1/SCDA/Apertures/InputMasks_v4/TelAp")
if symm is 'half':
if centobs is True:
telap_bin_dat_fname_tail = "TelAp_half_{0:s}{1:s}025cobs1gap{2:1d}_N{3:04d}.dat".format(
prim_key, secobs_key, gap, N)
telap_bin_dat_fname = os.path.join(opt_telap_dir, telap_bin_dat_fname_tail)
else:
telap_bin_dat_fname_tail = "TelAp_half_{0:s}{1:s}025cobs0gap{2:1d}_N{3:04d}.dat".format(
prim_key, secobs_key, gap, N)
telap_bin_dat_fname = os.path.join(opt_telap_dir, telap_bin_dat_fname_tail)
opt_telap_half = np.loadtxt(telap_bin_dat_fname)
opt_telap = np.concatenate((opt_telap_half[:,::-1], opt_telap_half),axis=1)
else:
if centobs is True:
telap_bin_dat_fname_tail = "TelAp_quart_{0:s}{1:s}025cobs1gap{2:1d}_N{3:04d}.dat".format(
prim_key, secobs_key, gap, N)
telap_bin_dat_fname = os.path.join(opt_telap_dir, telap_bin_dat_fname_tail)
else:
telap_bin_dat_fname_tail = "TelAp_quart_{0:s}{1:s}025cobs1gap{2:1d}_N{3:04d}.dat".format(
prim_key, secobs_key, gap, N)
telap_bin_dat_fname = os.path.join(opt_telap_dir, telap_bin_dat_fname_tail)
opt_telap_quart = np.loadtxt(telap_bin_dat_fname)
opt_telap = np.concatenate((np.concatenate((opt_telap_quart[::-1,::-1], opt_telap_quart[:,::-1]),axis=0),
np.concatenate((opt_telap_quart[::-1,:], opt_telap_quart),axis=0)), axis=1)
In [173]:
orig_LS = np.round(Lyot_stop_binned[L_bin/2-N:L_bin/2+N,L_bin/2-N:L_bin/2+N]).astype(int)
dz_width_fac = float(aligntol)/1000 # dark zone is +/- this fraction of pupil diameter
dz_width = np.ceil(2*N*dz_width_fac).astype(int) # dark zone is +/- this number of pixels in binned pupil array
print("Lyot plane dark zone width in binned array: +/- {:d} pixels".format(dz_width))
In [174]:
max_shift = dz_width
shift_range = range(-max_shift,max_shift+1,1)
[Xshifts, Yshifts] = np.meshgrid(shift_range, shift_range)
allowed_shifts = np.less_equal(Xshifts**2 + Yshifts**2, max_shift**2)
In [175]:
allowed_shifts
Out[175]:
In [176]:
XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])
len(XYshifts_allowed)
Out[176]:
In [177]:
fat_LS = np.ones(orig_LS.shape)
for (xshift,yshift) in XYshifts_allowed:
LS_shifted = np.roll(np.roll(orig_LS, yshift, 0), xshift, 1)
fat_LS *= LS_shifted
inv_thin_LS = np.ones(orig_LS.shape)
for (xshift,yshift) in XYshifts_allowed:
inv_LS_shifted = 1-np.roll(np.roll(orig_LS, yshift, 0), xshift, 1)
inv_thin_LS *= inv_LS_shifted
thin_LS = 1-inv_thin_LS
In [178]:
Lyot_dz = np.logical_xor(thin_LS, fat_LS)
L = Lyot_dz.shape[0]
if symm is 'half':
Lyot_dz_crop = Lyot_dz[:,N:]
else:
Lyot_dz_crop = Lyot_dz[N:,N:]
#plt.figure(figsize=(16,16))
#plt.imshow(Lyot_dz)
In [179]:
#plt.figure(figsize=(16,16))
#plt.imshow(1-np.floor(orig_telap))
#plt.imshow(np.logical_and(Lyot_dz,1-np.floor(orig_telap)))
print("None of the zero-valued telescope aperture points intersect with the Lyot constraint region? {:}".format(
~np.any(np.logical_and(Lyot_dz,1-np.round(opt_telap)))))
In [180]:
plt.figure(figsize=(16,16))
#plt.imshow(orig_telap)
#plt.imshow(np.floor(orig_telap).astype(int))
#plt.imshow(Lyot_dz-(1-np.floor(orig_telap)),cmap='gray')
plt.imshow(Lyot_dz-(1-np.round(opt_telap)),cmap='gray')
Out[180]:
In [181]:
max_shift_tol = 2*N*dz_width_fac
max_shift_tol_int = int(np.round(max_shift_tol))
test_shift = (max_shift_tol_int,0)
print("The LDZ accomomdates a translation +/-{0:.1f}% of D={1:d} pixels = +/-{2:.2f} pixels, up to +/-{3:d} whole pixels".format(
float(aligntol)/10, 2*N, max_shift_tol, max_shift_tol_int))
print("Testing an (x,y) translation of {0:} pixels. Within the design tolerance? {1:}".format(
test_shift, test_shift[0]**2 + test_shift[1]**2 <= max_shift_tol))
shift_LS = np.roll(np.roll(orig_LS, test_shift[0], axis=1), test_shift[1], axis=0)
LS_err_mask = np.ceil(np.abs(shift_LS - orig_LS)).astype(bool)
LDZ_valid = ~np.any(np.logical_and(LS_err_mask, ~Lyot_dz)) and ~np.any(np.logical_and(Lyot_dz,1-np.round(opt_telap)))
print("LDZ encompasses the LS transmission error region and does not overlap with tel pupil obscurations? {0:}".format(LDZ_valid))
print("Total unconstrained \"leak\" area after translation = {0:d} pixels".format(
int(np.sum(np.logical_and(LS_err_mask, ~Lyot_dz)))))
In [182]:
plt.figure(figsize=(16,6))
plt.subplot(131)
plt.imshow(LS_err_mask)
lims = plt.axis('off')
t=plt.title('Lyot transmission error region')
plt.subplot(132)
plt.imshow(~Lyot_dz)
lims = plt.axis('off')
t=plt.title('Inverse of LDZ')
plt.subplot(133)
plt.imshow(np.logical_and(LS_err_mask, ~Lyot_dz))
lims = plt.axis('off')
t=plt.title('Lyot leak region (black is good)')
In [183]:
if aligntol > 0 and LDZ_valid:
if symm is 'half':
LDZ_dat_fname_tail = "LDZ_half_{0:s}{1:02d}D{2:02d}_{3:s}{4:s}025cobs1Pad{5:02d}_Tol{6:02d}_N{7:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), prim_key,
secobs_key, pad, aligntol, N)
else:
LDZ_dat_fname_tail = "LDZ_quart_{0:s}{1:02d}D{2:02d}_{3:s}{4:s}025cobs1Pad{5:02d}_Tol{6:02d}_N{7:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), prim_key,
secobs_key, pad, aligntol, N)
LDZ_dat_fname = os.path.join(LS_dir, LDZ_dat_fname_tail)
if not os.path.exists(LDZ_dat_fname) or overwrite==True:
np.savetxt(LDZ_dat_fname, Lyot_dz_crop, fmt='%d', delimiter=" ")
print("Wrote binned, cropped LDZ array to {0:s}".format(LDZ_dat_fname))
else:
print("LDZ array {0:s} already exists, will not overwrite".format(LDZ_dat_fname))
In [184]:
Lyot_dz_crop.shape
Out[184]:
In [185]:
# os.remove(LDZ_dat_fname)
In [186]:
os.listdir(LS_dir)
Out[186]:
In [187]:
LS_full_fits_fname = LS_dat_fname_tail.replace("quart", "full")
LS_full_fits_fname = LS_full_fits_fname.replace("half", "full")
LS_full_fits_fname = LS_full_fits_fname.replace(".dat", ".fits")
LS_full_fits_fname = os.path.join(LS_dir, LS_full_fits_fname)
if not os.path.exists(LS_full_fits_fname) or overwrite==True:
LS_full_hdu = fits.PrimaryHDU(Lyot_stop_full_binary)
LS_full_hdu.writeto(LS_full_fits_fname, clobber=True)
print("Wrote full, rounded Lyot stop to {:s}".format(LS_full_fits_fname))
LDZ_full_fits_fname = LDZ_dat_fname_tail.replace("quart", "full")
LDZ_full_fits_fname = LDZ_full_fits_fname.replace("half", "full")
LDZ_full_fits_fname = LDZ_full_fits_fname.replace(".dat", ".fits")
LDZ_full_fits_fname = os.path.join(LS_dir, LDZ_full_fits_fname)
if not os.path.exists(LDZ_full_fits_fname) or overwrite==True:
LDZ_full_hdu = fits.PrimaryHDU(Lyot_dz.astype(int))
LDZ_full_hdu.writeto(LDZ_full_fits_fname, clobber=True)
print("Wrote full, rounded Lyot dark zone to {:s}".format(LDZ_full_fits_fname))