<jupyter_start><jupyter_text>How close is the solution to a shaped pupil? Take the sum of the absolute differences from rounded integer version<jupyter_code>nb_resid_frac = 100*np.sum(np.abs(A - np.round(A)))/np.sum(TelAp)
print("Non-binary residuals, as a percentage of clear telescope aperture area: {:.2}%".format(nb_resid_frac))
plt.figure(figsize=(12,12))
#plt.imshow(A*TelAp - A)
plt.imshow(A*TelAp)
plt.figure(figsize=(12,12))
plt.imshow(LS,vmin=0.5,vmax=1)<jupyter_output><empty_output><jupyter_text>Monochromatic PSF evaluation  Set dimensions and coordinates<jupyter_code>D = 1.
M_fp1 = FPM.shape[0]
fpm_rad = 4.
rho1 = 10.
fp2res = 8.
M_fp2 = int(np.ceil(rho1*fp2res))

# pupil plane
dx = (D/2)/N
dy = dx
xs = np.matrix(np.linspace(0.5,N-0.5,N)*dx)
ys = xs.copy()

# FPM
dmx = fpm_rad/M_fp1
dmy = dmx
mxs = np.matrix(np.linspace(0.5,M_fp1-0.5,M_fp1)*dmx)
mys = mxs.copy()

# FP2
dxi = 1/fp2res
xis = np.matrix(np.linspace(0.5,M_fp2-0.5,M_fp2)*dxi)
etas = xis<jupyter_output><empty_output><jupyter_text>Semi-analytical cosine propagation<jupyter_code>Psi_B = 4*dx*dy*np.dot(np.dot(np.cos(2*np.pi*np.dot(mys.T, ys)), A*TelAp ), np.cos(2*np.pi*np.dot(xs.T, mxs)))
Psi_B_stop = np.multiply(Psi_B, FPM)
Psi_C = A*TelAp - 4*dmx*dmy*np.dot(np.cos(2*np.pi*np.dot(ys.T, mys)), Psi_B_stop).dot(np.cos(2*np.pi*np.dot(mxs.T, xs)))
Psi_C_stop = np.multiply(Psi_C, LS)
Psi_D = 4*dx*dy*np.dot(np.cos(2*np.pi*np.dot(etas.T, ys)), Psi_C_stop ).dot(np.cos(2*np.pi*np.dot(xs.T, xis)))
Psi_D_0 = 4*dx*dy*np.dot(np.cos(2*np.pi*np.dot(etas.T, ys)), A*TelAp*LS ).dot(np.cos(2*np.pi*np.dot(xs.T, xis)))
Psi_D_0_peak = 4*np.sum(A*TelAp*LS)*dx*dy
Intens_B = np.power(np.absolute(Psi_B_stop), 2)
Intens_C = np.power(np.absolute(Psi_C_stop), 2)
Intens_D = np.power(np.absolute(Psi_D/Psi_D_0_peak), 2)
Intens_D_0 = np.power(np.absolute(Psi_D_0), 2)
Intens_D_0_peak = Psi_D_0_peak**2<jupyter_output><empty_output><jupyter_text>Final image<jupyter_code>plt.figure(figsize=(16,12))
plt.imshow(np.log10(Intens_D), vmin=-12, vmax=-8)
#plt.imshow(np.log10(Intens_D))
plt.colorbar()
ax = plt.gca()
ax.set_xticks(np.arange(0,M_fp2,4))
ax.set_xticklabels((np.arange(0,M_fp2,4)+0.5)*dxi)
ax.set_yticks(np.arange(0,M_fp2,4))
ax.set_yticklabels((np.arange(0,M_fp2,4)+0.5)*dxi)<jupyter_output><empty_output><jupyter_text>Airy throughput<jupyter_code>Psi_TelAp = 4*dx*dy*np.dot(np.cos(2*np.pi*np.dot(etas.T, ys)), TelAp ).dot(np.cos(2*np.pi*np.dot(xs.T, xis)))
Intens_TelAp = np.power(np.absolute(Psi_TelAp), 2)
Intens_TelAp_peak = (4*np.sum(TelAp)*dx*dy)**2
fwhm_ind_TelAp = np.greater_equal(Intens_TelAp, Intens_TelAp_peak/2)
fwhm_ind_APLC = np.greater_equal(Intens_D_0, Intens_D_0_peak/2)

fwhm_sum_TelAp = np.sum(Intens_TelAp[fwhm_ind_TelAp])
fwhm_sum_APLC = np.sum(Intens_D_0[fwhm_ind_APLC])

airy_thrupt = fwhm_sum_APLC/fwhm_sum_TelAp
print("Airy throughput = {:.1f}%".format(100*airy_thrupt))

plt.subplot(121)
plt.imshow(Intens_D_0[:10,:10])
plt.subplot(122)
plt.imshow(Intens_TelAp[:10,:10])<jupyter_output><empty_output><jupyter_text>First image<jupyter_code>#plt.figure(figsize=(16,12))
#plt.imshow(np.log10(Intens_B),vmin=-6)
#plt.colorbar()<jupyter_output><empty_output><jupyter_text>Lyot plane<jupyter_code>plt.figure(figsize=(16,12))
plt.imshow(np.log10(Intens_C))
plt.colorbar()<jupyter_output><empty_output><jupyter_text>Polychromatic PSF evaluation<jupyter_code>D = 1.
bw = 0.10
Nlambda = 3
M_fp1 = FPM.shape[0]
fpm_rad = 4.
rho1 = 10.
fp2res = 16.
M_fp2 = int(np.ceil(rho1*fp2res))

# pupil plane
dx = (D/2)/N
dy = dx
xs = np.matrix(np.linspace(0.5,N-0.5,N)*dx)
ys = xs.copy()

# FPM
dmx = fpm_rad/M_fp1
dmy = dmx
mxs = np.matrix(np.linspace(0.5,M_fp1-0.5,M_fp1)*dmx)
mys = mxs.copy()

# FP2
dxi = 1/fp2res
xis = np.matrix(np.linspace(0.5,M_fp2-0.5,M_fp2)*dxi)
etas = xis

# wavelength-dependent variables
wrs = np.linspace(1.-bw/2, 1.+bw/2, Nlambda)
Intens_D_polychrom = []
airy_thrupt_polychrom = []
fwhm_area_polychrom = []
for wr in wrs:
    Psi_B = 4*dx*dy/wr*np.dot(np.dot(np.cos(2*np.pi/wr*np.dot(mys.T, ys)), A*TelAp ), np.cos(2*np.pi/wr*np.dot(xs.T, mxs)))
    Psi_B_stop = np.multiply(Psi_B, FPM)
    Psi_C = A*TelAp - 4*dmx*dmy/wr*np.dot(np.cos(2*np.pi/wr*np.dot(ys.T, mys)), Psi_B_stop).dot(np.cos(2*np.pi/wr*np.dot(mxs.T, xs)))
    Psi_C_stop = np.multiply(Psi_C, LS)
    Psi_D = 4*dx*dy/wr*np.dot(np.cos(2*np.pi/wr*np.dot(etas.T, ys)), Psi_C_stop ).dot(np.cos(2*np.pi/wr*np.dot(xs.T, xis)))
    Psi_D_0 = 4*dx*dy/wr*np.dot(np.cos(2*np.pi/wr*np.dot(etas.T, ys)), A*TelAp*LS ).dot(np.cos(2*np.pi/wr*np.dot(xs.T, xis)))
    Psi_D_0_peak = 4*np.sum(A*TelAp*LS)*dx*dy/wr
    Intens_D_0 = np.power(np.absolute(Psi_D_0), 2)
    Intens_D_0_peak = Psi_D_0_peak**2
    Intens_D_polychrom.append( np.power(np.absolute(Psi_D/Psi_D_0_peak), 2) )
    
    Psi_TelAp = 4*dx*dy/wr*np.dot(np.cos(2*np.pi/wr*np.dot(etas.T, ys)), TelAp ).dot(np.cos(2*np.pi/wr*np.dot(xs.T, xis)))
    Intens_TelAp = np.power(np.absolute(Psi_TelAp), 2)
    Intens_TelAp_peak = (4*np.sum(TelAp)*dx*dy/wr)**2
    fwhm_ind_TelAp = np.greater_equal(Intens_TelAp, Intens_TelAp_peak/2)
    fwhm_ind_APLC = np.greater_equal(Intens_D_0, Intens_D_0_peak/2)

    fwhm_sum_TelAp = np.sum(Intens_TelAp[fwhm_ind_TelAp])
    fwhm_sum_APLC = np.sum(Intens_D_0[fwhm_ind_APLC])
    
    fwhm_area_polychrom.append(4*np.sum(fwhm_ind_APLC)/(fp2res**2))
    airy_thrupt_polychrom.append(fwhm_sum_APLC/fwhm_sum_TelAp)
print(airy_thrupt_polychrom)
airy_thrupt_bandavg = np.mean(airy_thrupt_polychrom)
fwhm_area_bandavg = np.mean(fwhm_area_polychrom)
print(airy_thrupt_bandavg)
print(fwhm_area_bandavg)
plt.figure(figsize=(16,12))
plt.imshow(np.log10(Intens_D_polychrom[0]), vmin=-12, vmax=-10)
#plt.imshow(np.log10(Intens_D))
plt.colorbar()
ax = plt.gca()
ax.set_xticks(np.arange(0,M_fp2,4))
ax.set_xticklabels((np.arange(0,M_fp2,4)+0.5)*dxi)
ax.set_yticks(np.arange(0,M_fp2,4))
ax.set_yticklabels((np.arange(0,M_fp2,4)+0.5)*dxi)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18962887646161666, 'fraction_numerical': 0.04253516353160481, 'mean_word_length': 6.286419753086419, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '259598', 'license': 'MIT', 'n_tokens_mistral': 3201, 'n_tokens_neox': 3020, 'n_words': 315, 'path': 'old_notebooks/eval_qpaplc_mono.ipynb', 'repo_name': 'neilzim/SCDA', 'size': '5901'}