"""Spherical-Sun optical benchmark, not a mission or NIAC image reproduction. Run: python reproduce.py. Requires numpy, scipy, matplotlib. Equations: Turyshev & Toth (2019), PRD100 084018, Eqs9,81-83. All files are written beside this script. No external image/data is required. """ from pathlib import Path import csv import json import platform import numpy as np import scipy from scipy.special import j0, j1, jn_zeros import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt OUT = Path(__file__).resolve().parent AU = 149597870700.0 RG = 2953.25 WAVELENGTH = 1e-6 Z = 650 * AU APERTURE = 1.0 ALPHA = 2 * np.pi / WAVELENGTH * np.sqrt(2 * RG / Z) GAIN = 4 * np.pi**2 * RG / WAVELENGTH N = 24 WIDTH = 6.0 DX = WIDTH / N def pupil_average(rho, radial=96, angular=192, aperture=APERTURE): """Disk-area average of relative J0^2; no gain or power units attached.""" g, w = np.polynomial.legendre.leggauss(radial) a = aperture / 2 r = (g + 1) * a / 2 rw = w * a / 2 * r * 2 / a**2 theta = 2 * np.pi * np.arange(angular) / angular result = [] for chunk in np.array_split(np.asarray(rho), max(1, int(np.ceil(len(rho) / 32)))): distance = np.sqrt(np.maximum(0, chunk[:, None, None]**2 + r[None, :, None]**2 + 2 * chunk[:, None, None] * r[None, :, None] * np.cos(theta))) result.extend((j0(ALPHA * distance)**2).mean(axis=2) @ rw) return np.asarray(result) def cell_kernel(radial_grid, averaged_psf, quadrature): """Integrate uniform source cells; retain physical cell-area weights.""" q, w = np.polynomial.legendre.leggauss(quadrature) q = q * DX / 2 weights = np.outer(w, w) * DX**2 / 4 k = np.zeros((N, N)) for iy in range(N): for ix in range(iy, N): rho = np.hypot(ix * DX + q[:, None], iy * DX + q[None, :]) k[iy, ix] = k[ix, iy] = np.sum(np.interp(rho, radial_grid, averaged_psf) * weights) return k def operator(k): yy, xx = np.indices((N, N)) return k[np.abs(yy.ravel()[:, None] - yy.ravel()), np.abs(xx.ravel()[:, None] - xx.ravel())] def write_csv(name, header, rows): with (OUT / name).open('w', newline='', encoding='utf-8') as f: writer = csv.writer(f) writer.writerow(header) writer.writerows(rows) def main(): plt.rcParams.update({'font.family': 'sans-serif', 'font.sans-serif': ['Arial', 'DejaVu Sans', 'sans-serif'], 'font.size': 12, 'axes.spines.top': False, 'axes.spines.right': False, 'svg.fonttype': 'none', 'savefig.facecolor': 'white'}) rho = np.linspace(0, 9, 9001) avg = pupil_average(rho) check_rho = np.array([0, .013, .049, .125, .25, .499, .5, .501, .75, 1, 2.345, 4, 6, 8.5, 9]) pupil_low = pupil_average(check_rho) pupil_high = pupil_average(check_rho, 128, 256) pupil_error = float(np.max(np.abs(pupil_low - pupil_high))) exact_axis = float(j0(ALPHA * APERTURE / 2)**2 + j1(ALPHA * APERTURE / 2)**2) k8 = cell_kernel(rho, avg, 8) k16 = cell_kernel(rho, avg, 16) k32 = cell_kernel(rho, avg, 32) coarse = cell_kernel(rho[::2], avg[::2], 32) A = operator(k32) source_y, source_x = np.meshgrid((np.arange(N) + .5) * DX - WIDTH / 2, (np.arange(N) + .5) * DX - WIDTH / 2, indexing='ij') source = .12 + .75 * np.exp(-((source_x + 1.1)**2 + (source_y - .6)**2) / .55) source += .5 * ((source_x > .5) & (source_x < 1.5) & (source_y > -1.5) & (source_y < 1.5)) source[(np.abs(source_y + .9) < .2) & (source_x > -.5)] = .03 source /= source.max() truth = source.ravel() clean = A @ truth U, singular, Vt = np.linalg.svd(A, full_matrices=False) rng = np.random.default_rng(271828) normal = rng.normal(size=N*N) cases = [] for noise in [0, .001, .01]: measured = clean + noise * clean.max() * normal for regularization in [0, .01, .1]: lam = regularization * singular[0] recovered = Vt.T @ ((singular / (singular**2 + lam**2)) * (U.T @ measured)) error = recovered - truth cases.append({'id': f'n{noise:g}-r{regularization:g}', 'noise_fraction': noise, 'regularization_fraction': regularization, 'source_rmse': float(np.sqrt(np.mean(error**2))), 'source_relative_l2_error': float(np.linalg.norm(error)/np.linalg.norm(truth)), 'realized_noise_rms': float(np.sqrt(np.mean((measured-clean)**2))), 'measurement_relative_residual': float(np.linalg.norm(A@recovered-measured)/np.linalg.norm(measured)), 'recovered_min': float(recovered.min()), 'recovered_max': float(recovered.max()), 'negative_cells': int(np.sum(recovered < 0)), 'above_one_cells': int(np.sum(recovered > 1)), 'measured': measured.tolist(), 'recovered': recovered.tolist(), 'error': error.tolist()}) selected = next(c for c in cases if c['noise_fraction']==.01 and c['regularization_fraction']==.01) first_zero = float(jn_zeros(0, 1)[0] / ALPHA) checks = { 'analytic_on_axis_relative_average': exact_axis, 'quadrature_on_axis_absolute_error': float(abs(avg[0] - exact_axis)), 'pupil_refinement_max_absolute_error_15_radii': pupil_error, 'cell_8_to_16_relative_frobenius': float(np.linalg.norm(operator(k8)-operator(k16))/np.linalg.norm(A)), 'cell_16_to_32_relative_frobenius': float(np.linalg.norm(operator(k16)-A)/np.linalg.norm(A)), 'radial_interpolation_2mm_to_1mm_relative_frobenius': float(np.linalg.norm(operator(coarse)-A)/np.linalg.norm(A)), 'operator_symmetry_max_absolute': float(np.max(np.abs(A-A.T))), 'operator_min_entry': float(A.min()), 'same_model_no_noise_unregularized_rmse': cases[0]['source_rmse'], 'condition_number': float(singular[0]/singular[-1]), } alternate_forward = operator(k16) @ truth alternate_recovery = Vt.T @ ((U.T @ alternate_forward) / singular) checks['different_quadrature_no_noise_unregularized_rmse'] = float(np.sqrt(np.mean((alternate_recovery-truth)**2))) mismatch_average = pupil_average(rho, aperture=1.02) mismatch_signal = operator(cell_kernel(rho, mismatch_average, 32)) @ truth mismatch_recovered = Vt.T @ ((singular/(singular**2+(.01*singular[0])**2))*(U.T@mismatch_signal)) checks['aperture_mismatch_1p02m_generated_1m_inverse_no_noise_r0p01_rmse'] = float(np.sqrt(np.mean((mismatch_recovered-truth)**2))) checks['cell_refinement_forward_pattern_relative_l2'] = float(np.linalg.norm(alternate_forward-clean)/np.linalg.norm(clean)) assert checks['quadrature_on_axis_absolute_error'] < 1e-12 assert pupil_error < 1e-10 assert checks['cell_16_to_32_relative_frobenius'] < 2e-4 assert checks['radial_interpolation_2mm_to_1mm_relative_frobenius'] < 2e-4 assert checks['same_model_no_noise_unregularized_rmse'] < 1e-8 methods = {'classification': 'Original finite-source numerical benchmark using published optical equations; not a mission forecast or reproduction of published planet images', 'sources': [{'doi': '10.1103/PhysRevD.100.084018', 'url': 'https://arxiv.org/pdf/1908.01948', 'locators': 'Eq9 PDF p5; Eqs81-82 PDF p18; Eq83 PDF p19'}, {'doi': '10.1103/PhysRevD.103.124038', 'url': 'https://arxiv.org/abs/2012.05477', 'locators': 'Eqs1,4-6; forward and inverse image models'}], 'parameters': {'wavelength_m': WAVELENGTH,'z_AU':650,'rg_m':RG,'AU_m':AU,'aperture_diameter_m':APERTURE, 'alpha_per_m':float(ALPHA),'peak_point_gain':float(GAIN),'first_zero_radius_m':first_zero, 'on_axis_aperture_gain':float(exact_axis*GAIN),'grid_side':N,'field_width_m':WIDTH,'cell_width_m':DX, 'radial_kernel_step_m':.001,'radial_kernel_max_m':9,'cell_quadrature_order':32, 'pupil_radial_order':96,'pupil_angular_samples':192,'noise_seed':271828}, 'operator': 'A_ij = integral over square source cell j of aperture-area-average J0(alpha|x_i-u+r|)^2 d2u. Fixed observation centers; zero source outside 6m square; no circular convolution or row normalization. Source cell brightness is dimensionless; clean signal has relative brightness times m2 units, with common gain omitted.', 'inverse': 'Unconstrained Tikhonov: V diag(s/(s2+lambda2)) U^T y; lambda = control times largest singular value. No positivity clipping in computations or downloadable data. Exact same operator generates and reconstructs measurements; zero-noise recovery is a code consistency check, not physical validation.', 'noise': 'Fixed Gaussian sample shared across cases, sigma = selected fraction times maximum clean signal. It is analyst-selected noise, not a photon count, corona calculation or the published output SNR7.', 'mismatch_check': 'Additional no-noise signal generated with a 1.02m aperture and reconstructed using the nominal 1m aperture, regularization0.01. This tests one deliberately chosen PSF error, not a measured uncertainty or comprehensive robustness analysis.', 'display': 'True and recovered brightness use fixed [0,1] scale, error [-1,1]. Out-of-range values saturate only for display; extrema/counts and unmodified data are supplied. Measurement color scale fixed [0,max clean].', 'limitations': ['Artificial small patch, zero exterior: omits contamination from the rest of a planet.', 'Spherical static Sun and monochromatic source; no solar multipoles, pointing error or PSF mismatch.', 'No telescope focal-plane coronagraph/detector model, spectral response, photon budget, real corona noise or exposure prediction.', 'No rotating planet, cloud evolution, scan timing or navigation dynamics.', 'No source-matched NIAC image reproduction or achieved planetary resolution.'], 'checks': checks,'software': {'python':platform.python_version(),'numpy':np.__version__,'scipy':scipy.__version__,'matplotlib':matplotlib.__version__}} (OUT/'methods.json').write_text(json.dumps(methods,indent=2)+'\n',encoding='utf-8') bundle = {'n':N,'width_m':WIDTH,'source':truth.tolist(),'clean':clean.tolist(),'measurement_scale':float(clean.max()),'cases':cases,'checks':checks} (OUT/'cases.json').write_text(json.dumps(bundle,separators=(',',':'))+'\n',encoding='utf-8') write_csv('point-and-aperture-response.csv',['radius_m','point_relative_response','aperture_mean_relative_response'],zip(rho,j0(ALPHA*rho)**2,avg)) write_csv('source-and-clean.csv',['row','column','x_m','y_m','source_relative_brightness','clean_relative_signal_m2'], ((i//N,i%N,source_x.ravel()[i],source_y.ravel()[i],truth[i],clean[i]) for i in range(N*N))) write_csv('operator-cell-kernel.csv',['absolute_row_offset','absolute_column_offset','integrated_relative_kernel_m2'], ((y,x,k32[y,x]) for y in range(N) for x in range(N))) write_csv('case-metrics.csv',['noise_fraction','regularization_fraction','source_rmse','measurement_relative_residual','recovered_min','recovered_max','negative_cells','above_one_cells'], ([c[k] for k in ['noise_fraction','regularization_fraction','source_rmse','measurement_relative_residual','recovered_min','recovered_max','negative_cells','above_one_cells']] for c in cases)) # Publication-style curves: no raster traced from the source papers. for narrow in [False,True]: plt.rcParams.update({'font.size':22 if narrow else 12}) fig, axs = plt.subplots(2,1,figsize=(6,12) if narrow else (10,9),layout='constrained') small=rho<=.3 axs[0].plot(rho[small]*100,j0(ALPHA*rho[small])**2,color='#9b288f',lw=2) axs[0].axvline(first_zero*100,color='#186b8c',ls='--') axs[0].set(xlabel='Offset (cm)',ylabel='Relative point response',title='A point spreads\ninto rings' if narrow else 'A point spreads into rings') axs[0].text(.96,.85,f'First zero: {first_zero*100:.2f} cm',transform=axs[0].transAxes,ha='right') large=rho<=2 axs[1].plot(rho[large],j0(ALPHA*rho[large])**2,color='#9b288f',alpha=.6,label='Point sample') axs[1].plot(rho[large],avg[large],color='#186b8c',lw=2,label='1 m aperture mean') axs[1].set(xlabel='Offset (m)',ylabel='Relative response',yscale='log',ylim=(1e-5,1.1),title='The telescope averages\nmany rings' if narrow else 'The telescope averages many rings') if narrow: axs[0].set_xticks([0,10,20,30]);axs[1].set_xticks([0,1,2]) axs[1].legend(loc='upper right'); [ax.grid(alpha=.2) for ax in axs] fig.suptitle('Calculated optical response\n1 µm · 650 AU' if narrow else 'Calculated optical response · 1 µm · 650 AU',fontsize=22 if narrow else 14) stem='point-response-mobile' if narrow else 'point-response' fig.savefig(OUT/(stem+'.svg'));fig.savefig(OUT/(stem+'.png'),dpi=130) if not narrow:fig.savefig(OUT/'point-response.pdf') plt.close(fig) for narrow in [False,True]: plt.rcParams.update({'font.size':22 if narrow else 12}) fig,axs=plt.subplots(4,1,figsize=(6,23),layout='constrained') if narrow else plt.subplots(2,2,figsize=(10,9),layout='constrained') arrays=[truth,selected['measured'],selected['recovered'],selected['error']] titles=['Known synthetic source','Computed measurements + 1% noise','Recovered · regularization 0.01','Recovered minus known source'] if narrow:titles=['Known synthetic\nsource','Computed measurements\n+ 1% noise','Recovered\nregularization 0.01','Recovered minus\nknown source'] for i,(ax,arr,title) in enumerate(zip(np.ravel(axs),arrays,titles)): image=ax.imshow(np.reshape(arr,(N,N)),origin='lower',extent=(-3,3,-3,3),cmap='coolwarm' if i==3 else 'viridis', vmin=-1 if i==3 else 0,vmax=1 if i!=1 else clean.max(),interpolation='nearest') ax.set(title=title,xlabel='Projected x (m)',ylabel='Projected y (m)') if narrow:ax.set_xticks([-3,0,3]);ax.set_yticks([-3,0,3]) fig.colorbar(image,ax=ax,shrink=.75,label='Brightness error' if i==3 else ('Relative signal × m²' if i==1 else 'Relative brightness')) fig.suptitle((f'Finite optical benchmark\nSource RMSE {selected["source_rmse"]:.3f}\n6 m patch; zero outside' if narrow else f'Finite optical benchmark · source RMSE {selected["source_rmse"]:.3f}\nArtificial 6 m patch; source is zero outside'),fontsize=22 if narrow else 14) stem='reconstruction-mobile' if narrow else 'reconstruction' fig.savefig(OUT/(stem+'.svg'));fig.savefig(OUT/(stem+'.png'),dpi=130) if not narrow:fig.savefig(OUT/'reconstruction.pdf') plt.close(fig) print(json.dumps(checks,indent=2)) if __name__=='__main__': main()