**1. 収差モデリングとシミュレーションフレームワーク**
光学収差のシミュレーションは幾何光学と波動光学の原理を統合して実現します。以下に主要なモジュールを示します。
**1.1 Zernike多項式による波面表現**
Zernike多項式は光学収差を記述する基礎的な関数系です。極座標系での定義は以下の通りです：
Z_n^m(ρ,θ) = R_n^m(ρ) × e^(imθ)
主要な収差に対応する項：
- 球面収差：Z_4^0項
- コマ収差：Z_3^1とZ_3^−1項
- 非点収差：Z_4^±2項
- 像面湾曲：Z_4^0の高次成分
MATLAB実装コード**：
function wavefront = calc_zernike(order, freq, radial_coord, angular_coord)
% 入力検証
if abs(freq) > order || mod(order - abs(freq), 2) ~= 0
wavefront = zeros(size(radial_coord));
return;
end
% 規格化係数の計算
norm_coeff = sqrt(2 * (order + 1) / (1 + (freq == 0)));
% 径向多項式の構築
radial_term = zeros(size(radial_coord));
for k = 0:(order - abs(freq))/2
num = (-1)^k * factorial(order - k);
den = factorial(k) * factorial((order + abs(freq))/2 - k) * factorial((order - abs(freq))/2 - k);
radial_term = radial_term + num/den * radial_coord.^(order - 2*k);
end
% 方位角成分との結合
if freq >= 0
wavefront = norm_coeff * radial_term .* cos(freq * angular_coord);
else
wavefront = norm_coeff * radial_term .* sin(abs(freq) * angular_coord);
end
end
**1.2 光線追跡アルゴリズム**
近軸光学公式に基づく光線追跡により、光学系通過後の波面歪みを計算します。
function [xf, yf] = trace_ray_surface(x0, y0, inc_angle, n_medium, n_surface, radius)
% 球面における光線の屈折計算
% 面法線ベクトルの角度
normal_angle = atan2(y0, radius - x0);
% 実効入射角
effective_inc = inc_angle - normal_angle;
% スネルの法則による屈折角
sin_refr = (n_medium/n_surface) * sin(effective_inc);
refr_angle = asin(sin_refr);
% 出射方向の決定
exit_direction = refr_angle + normal_angle;
% 新しい交点座標
xf = x0 + radius * (1 - cos(exit_direction));
yf = y0 + radius * sin(exit_direction);
end
**2. 収差シミュレーションの実装**
**2.1 単色収差のシミュレーション**
球面収差の可視化**：
% パラメータ設定
order = 4; freq = 0; % Z4^0は球面収差
r_vec = linspace(0, 1, 200);
phi_vec = linspace(0, 2*pi, 200);
[R_grid, Phi_grid] = meshgrid(r_vec, phi_vec);
% 波面計算
spherical_wave = calc_zernike(order, freq, R_grid, Phi_grid);
% 3Dプロット
figure;
surf(R_grid.*cos(Phi_grid), R_grid.*sin(Phi_grid), spherical_wave*5);
title('球面収差の波面分布');
xlabel('X座標(mm)'); ylabel('Y座標(mm)'); zlabel('波面収差(λ)');
colormap jet; colorbar;
コマ収差の解析**：
% Z3^1とZ3^−1の合成
zernike_3_1 = calc_zernike(3, 1, R_grid, Phi_grid);
zernike_3_minus1 = calc_zernike(3, -1, R_grid, Phi_grid);
coma_pattern = zernike_3_1 + zernike_3_minus1;
% ベクトル場表示
[x_mesh, y_mesh] = meshgrid(linspace(-1,1,50), linspace(-1,1,50));
figure;
quiver(x_mesh, y_mesh, gradient(coma_pattern));
title('コマ収差の光路ベクトル場');
axis equal; grid on;
**2.2 色収差の解析**
複数波長の重ね合わせによる軸上色収差のシミュレーション：
wavelengths = [0.486, 0.546, 0.656]; % 青、緑、赤の波長(μm)
focal_len = 100; % 焦点距離(mm)
color_map = lines(length(wavelengths));
figure; hold on;
for idx = 1:length(wavelengths)
% 各波長の焦点位置計算
focal_shift = focal_len * (wavelengths(idx) - wavelengths(2))/wavelengths(2);
spot_size = 0.01 * idx;
theta = linspace(0, 2*pi, 100);
plot(focal_len + focal_shift + spot_size*cos(theta), ...
spot_size*sin(theta), 'Color', color_map(idx,:), 'LineWidth', 2);
end
title('色収差による焦点ずれ');
xlabel('光軸方向(mm)'); ylabel('径方向(mm)');
legend('F線(486nm)', 'e線(546nm)', 'C線(656nm)'); hold off;
**3. 高度な機能モジュール**
**3.1 乱流位相スクリーン**
大気乱流を模擬するためのランダム位相スクリーン生成：
function turbulence_screen = generate_turbulence_screen(grid_size, aperture_dia, fried_param)
% 空間周波数領域の設定
delta_f = 1 / aperture_dia;
freq_range = (-grid_size/2:grid_size/2-1) * delta_f;
[fx, fy] = meshgrid(freq_range);
freq_magnitude = sqrt(fx.^2 + fy.^2);
% ボン/モナジノ位相スペクトル
power_spectrum = 0.023 * fried_param^(-5/3) * freq_magnitude.^(-11/3);
power_spectrum(grid_size/2+1, grid_size/2+1) = 0; % DC成分除去
% ランダム位相生成
random_phase = exp(1i * 2*pi * rand(size(freq_magnitude)));
turbulence_screen = ifft2(sqrt(power_spectrum) .* random_phase);
turbulence_screen = turbulence_screen / std(turbulence_screen(:));
end
**3.2 収差補正シミュレーション**
複合レンズによる球面収差の補正効果の比較：
% 二重膠合レンズ設計パラメータ
curvatures = [100, -50, 80]; % 各面の曲率半径(mm)
refractive_indices = [1.5, 1.6, 1.5]; % 対応する屈折率
% 波面収差の計算
[wavefront_before, wavefront_after] = evaluate_correction(curvatures, refractive_indices, 25);
% 結果の可視化
figure;
subplot(1,2,1); imagesc(wavefront_before); colorbar;
title('補正前の波面収差'); xlabel('瞳面X(mm)'); ylabel('瞳面Y(mm)');
subplot(1,2,2); imagesc(wavefront_after); colorbar;
title('補正後の波面収差'); xlabel('瞳面X(mm)'); ylabel('瞳面Y(mm)');
**4. 可視化と解析ツール**
**4.1 三次元波面表示**
% 収差のある波面の生成
[axis_x, axis_y] = meshgrid(linspace(-1,1,400));
base_sphere = sqrt(max(0, 1 - axis_x.^2 - axis_y.^2));
aberration_term = 0.15 * sin(4*pi*axis_x) .* cos(2*pi*axis_y);
distorted_wavefront = base_sphere + aberration_term;
% 高品質サーファスプロット
figure;
surf(axis_x, axis_y, distorted_wavefront, 'EdgeColor', 'none');
title('収差を含む波面形状');
xlabel('X軸(mm)'); ylabel('Y軸(mm)'); zlabel('光軸方向変位(mm)');
colormap parula; camlight; lighting gouraud;
**4.2 スポットダイアグラムとRay Fanプロット**
% 光線束の追跡
incident_heights = linspace(-12, 12, 80);
focal_plane = 150; % 像面位置(mm)
[x_spots, y_spots] = trace_ray_bundle(incident_heights, focal_plane);
% スポットダイアグラム
figure; plot(x_spots, y_spots, 'ko', 'MarkerSize', 4);
title('スポットダイアグラム'); axis equal; grid on;
xlabel('X方向像面位置(mm)'); ylabel('Y方向像面位置(mm)');
% Ray Fan図
figure; plot(incident_heights, x_spots*1000, 'b-', 'LineWidth', 1.5); hold on;
plot(incident_heights, y_spots*1000, 'r--', 'LineWidth', 1.5);
title('Ray Fan図'); xlabel('入射高さ(mm)'); ylabel('像面ずれ(μm)');
legend('X成分', 'Y成分'); grid on;
**5. 完全なシミュレーションワークフロー例**
% ステップ1: 基本パラメータ設定
lambda_design = 0.55e-3; % 設計波長(mm)
aperture_size = 25; % 有効口径(mm)
efl_target = 100; % 目標焦点距離(mm)
% ステップ2: 収差成分の生成
[rho_mesh, theta_mesh] = meshgrid(linspace(0,1,256), linspace(0,2*pi,256));
zernike_coeff = 0.08; % 収差係数
wavefront_error = zernike_coeff * calc_zernike(4, 0, rho_mesh, theta_mesh);
% ステップ3: 光学系による光線追跡
num_trace_rays = 500;
[final_x, final_y, final_u, final_v] = optical_system_trace(wavefront_error, num_trace_rays);
% ステップ4: 性能評価指標の算出
rms_wavefront = sqrt(mean(wavefront_error(:).^2));
strehl_estimate = exp(-(2*pi * rms_wavefront)^2);
fprintf('RMS波面収差: %.4fλ\n', rms_wavefront);
fprintf('Strehl比: %.3f\n', strehl_estimate);
**6. エンジニアリング応用の拡張**
**6.1 適応光学的補正**
% 測定波面のZernike係数分解
measured_wavefront = acquire_wavefront_sensor_data();
[max_order, ~] = size(measured_wavefront);
zernike_orders = 6; % 補正次数
% 係数ベクトルの計算
z_coeffs = zeros(zernike_orders*(zernike_orders+1)/2, 1);
idx = 0;
for n_val = 0:zernike_orders
for m_val = -n_val:2:n_val
if abs(m_val) <= n_val
idx = idx + 1;
basis = calc_zernike(n_val, m_val, rho_mesh, theta_mesh);
z_coeffs(idx) = sum(measured_wavefront(:) .* basis(:)) / sum(basis(:).^2);
end
end
end
% デフォーカス補正の適用
corrected_wavefront = measured_wavefront - z_coeffs(4) * calc_zernike(2, 0, rho_mesh, theta_mesh);
**6.2 マルチチャネル並列処理**
% 波長帯域の並列シミュレーション
spectral_bandwidth = linspace(0.4, 0.7, 10); % 可視光帯域(μm)
parpool(6); % 6ワーカー並列起動
% 分散型計算
spectral_results = cell(size(spectral_bandwidth));
parfor idx = 1:length(spectral_bandwidth)
spectral_results{idx} = polychromatic_simulation(spectral_bandwidth(idx), aperture_size);
end
% 結果の統合
polychromatic_psf = zeros(size(spectral_results{1}));
for idx = 1:length(spectral_results)
polychromatic_psf = polychromatic_psf + spectral_results{idx};
end
imagesc(polychromatic_psf); title('ポリクロマティックPSF'); colorbar;
delete(gcp); % プールの解放