In [1]:
import numpy as np
#import skimage.transform
import astropy.io.fits as fits
import matplotlib.pyplot as plt
import os
%pylab inline --no-import-all
matplotlib.rcParams['image.origin'] = 'lower'
matplotlib.rcParams['image.interpolation'] = 'nearest'
matplotlib.rcParams['image.cmap'] = 'gray'
In [2]:
telap_key_map = {'hex1':'hex1', 'hex2':'hex2', 'hex3':'hex3', 'hex4':'hex4',
'key24':'keystone24', 'pie08':'piewedge8', 'pie12':'piewedge12'}
secobs_key_map = {'Cross':'cross', 'X':'x'}
In [336]:
overwrite = False
N = 250 # pupil array quadrant width after binning
iD = 30
oD = 76
aligntol = 2 # units of thousandths of pupil diameter
pad = 0 # units of thousandths of pupil diameter
L = 1000 # starting array width
D = 1000 # starting pupil diameter
symm = 'quart' # set to either 'quart' or 'half'
shape = 'ann' # set to 'ann' or 'hex'
edge = 'round' # set to 'gray' or 'round'
secobs_key = "X"
In [337]:
#LS_dir = os.path.abspath("../InputMasks/LS")
LS_dir = os.path.abspath("/astro/opticslab1/SCDA/Apertures/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 [338]:
#telap_dir = os.path.abspath('../Apertures/JPL/offset_masks')
telap_dir = os.path.abspath('/astro/opticslab1/SCDA/Apertures/JPL/offset_masks')
secobs_fname = os.path.join(telap_dir, "{0:s}_spiders_{1:04d}pix_2.5cm_offset.fits".format(
secobs_key_map[secobs_key], D))
secobs_hdulist = fits.open(secobs_fname, "readonly")
secobs = secobs_hdulist[0].data
secobs_hdulist.close()
In [339]:
plt.figure(figsize=(10,10))
plt.imshow(secobs)
Out[339]:
In [340]:
max_shift = int(round(D*float(pad)/1000))
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)
XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])
len(XYshifts_allowed)
Out[340]:
In [341]:
padded_secobs_accum = np.ones(secobs.shape)
for (xshift,yshift) in XYshifts_allowed:
secobs_shifted = np.roll(np.roll(secobs, yshift, 0), xshift, 1)
padded_secobs_accum *= secobs_shifted
In [342]:
#plt.figure(figsize=(10,10))
#plt.imshow(padded_secobs_accum)
In [343]:
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)
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_secobs_accum)
Lyot_stop[inside_ann] = 0.
Lyot_stop[outside_ann] = 0.
In [344]:
plt.figure(figsize=(12,12))
plt.imshow(Lyot_stop)
Out[344]:
In [345]:
#scalefac = float(N)/(D/2)
#Lyot_stop_binned = skimage.transform.rescale(Lyot_stop, scalefac, order=1)
In [346]:
N_orig = D/2
scalefac = int(N_orig/N)
print("Binning the original aperture 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)
In [347]:
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_crop).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
In [348]:
plt.figure(figsize=(10,10))
plt.imshow(Lyot_stop_binned_crop_binary)
#plt.imshow(np.floor(Lyot_stop_binned_quad))
Out[348]:
In [349]:
if symm is 'half':
LS_quart_dat_fname_tail = "LS_half_ann{0:02d}D{1:02d}_{2:s}025Pad{3:02d}_N{4:04d}.dat".format(
int(round(100*inD)), int(round(100*outD)), secobs_key, pad, N)
else:
LS_quart_dat_fname_tail = "LS_quart_ann{0:02d}D{1:02d}_{2:s}025Pad{3:02d}_N{4:04d}.dat".format(
int(round(100*inD)), int(round(100*outD)), secobs_key, pad, N)
LS_quart_dat_fname = os.path.join(LS_dir, LS_quart_dat_fname_tail)
#np.savetxt(LS_quart_dat_fname, Lyot_stop_binned_quad, fmt='%.6f', delimiter=" ")
np.savetxt(LS_quart_dat_fname, Lyot_stop_binned_crop_binary, fmt='%d', delimiter=" ")
print("Wrote binned, cropped LS array to {0:s}".format(LS_quart_dat_fname))
In [350]:
os.listdir(LS_dir)
Out[350]:
In [351]:
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 [352]:
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 [353]:
allowed_shifts
Out[353]:
In [354]:
XYshifts_allowed = zip(Xshifts[allowed_shifts], Yshifts[allowed_shifts])
len(XYshifts_allowed)
Out[354]:
In [355]:
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 [356]:
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:]
In [357]:
plt.figure(figsize=(16,16))
plt.imshow(Lyot_dz)
Out[357]:
In [358]:
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))
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 [359]:
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}025Pad{4:02d}_Tol{5:02d}_N{6:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), secobs_key, pad, aligntol, N)
else:
LDZ_dat_fname_tail = "LDZ_quart_{0:s}{1:02d}D{2:02d}_{3:s}025Pad{4:02d}_Tol{5:02d}_N{6:04d}.dat".format(
shape, int(round(100*inD)), int(round(100*outD)), 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 [360]:
os.listdir(LS_dir)
Out[360]:
In [ ]: