aber3_data_opt_sys_inf.m

用matlab实现光学系统像差计算

function [SphT, ComaS, AstTS, Ptz, Dist]=aber3_data_opt_sys_inf(opt_sys, flag_xy, i_c)
%ABER3_DATA_OPT_SYS_INF 3rd order aberration contributions for an object at infinity
%	ABER3_DATA_OPT_SYS_INF returns the 3rd order aberration contributions of each surface of an optical system
%		(also with non-spherical surfaces) for an object point at infinity 
%		(including entrance pupil).
%		See Malacara, D., Handbook of Lens-Design, Marcel Dekker, 1994, pp 132, 169, 175, 178, 196)
%		SPHT: 3rd order spherical, transversal (SPHL=SPHT*U_K')
%		PTZ: Petzval (1/R_PTZ=2*SUM(PTZ)/H_K'^2)
%		COMAS: coma sagittal (COMAT=3*COMAS)
%		ASTTS: astigmatism, transversal, sagittal (ASTTT=3*ASTTS, ASTLS=ASTTS*U_K', ASTLT=3*ASTLS))
%		DIST: distortion
%
%	The optical system is given by:
%		OPT_SYS.N_SURF (number of optical surfaces including the detector plane)
%		OPT_SYS.C (N_SURF-element array of the curvatures of the surfaces including the detector plane)
%		OPT_SYS.D (N_SURF-element array of the distances between the surfaces; 
%			the first distance OPT_SYS.D(1) is the distance from a starting surface 
%			to the first refracting surface, i.e. OPT_SYS.D(1)-OPT_SYS.D_EP is the distance between the 
%			entrance pupil and the first refracting surface
%		OPT_SYS.N (OPT_SYS.N_SURF-element array of the refraction indices of the media)
%		OPT_SYS.U_0 (angle of the principal ray in rad)
%		OPT_SYS.H_EP (height (radius) of the entrance pupil)
%		OPT_SYS.D_EP (distance of the entrance pupil from the first (non-refracting) surface)
%
%	An additional optional input parameter for the ray trace
%		FLAG_XY (specifies the direction, i.e. FLAG_XY=1 the x-direction, FLAG_XY=0 the y-direction; default FLAG_XY=1)
%		I_C (an integer between 1 and N_C which specifies the wavelength, i.e. which column of refraction indices
%			is used for the calculation of the efl; default I_C=1)
%
%   See also PAR_DATA_INF, EFL_OPT_SYS, BFL_OPT_SYS
%
%   Peter Ott 30.11.2000
%   Copyright (c) by Peter Ott

if nargin < 3, i_c = 1; end
if nargin < 2, flag_xy = 1; end

%marginal ray (Randstrahl)
nuh_1=[0; opt_sys.h_EP];
[nuh_r_x,nuh_r_y]=par_trace_opt_sys(opt_sys, nuh_1, nuh_1, i_c);
if flag_xy==1, nuh_r=nuh_r_x; else, nuh_r=nuh_r_y; end   
%principal ray (Hauptstrahl)
nuh_1=[opt_sys.n(1,i_c)*opt_sys.u_0; (opt_sys.d(1)-opt_sys.d_EP)*opt_sys.u_0];
[nuh_h_x,nuh_h_y]=par_trace_opt_sys(opt_sys, nuh_1, nuh_1, i_c);
if flag_xy==1, nuh_h=nuh_h_x; else, nuh_h=nuh_h_y; end   

lag=nuh_h(2,1)*nuh_r(1,1)-nuh_r(2,1)*nuh_h(1,1);
efl=lag*opt_sys.n(opt_sys.n_surf, i_c)/(nuh_r(1,opt_sys.n_surf)*nuh_h(1,1)-nuh_r(1,1)*nuh_h(1,opt_sys.n_surf));

y=nuh_r(2,1:opt_sys.n_surf-1)';
n=opt_sys.n(1:opt_sys.n_surf-1,i_c);
n_s=opt_sys.n(2:opt_sys.n_surf,i_c);
u=nuh_r(1,1:opt_sys.n_surf-1)'./n;
u_s=nuh_r(1,2:opt_sys.n_surf)'./n_s;
i=u+y.*opt_sys.c(1:opt_sys.n_surf-1);

y_b=nuh_h(2,1:opt_sys.n_surf-1)';
u_b=nuh_h(1,1:opt_sys.n_surf-1)'./n;
i_b=u_b+y_b.*opt_sys.c(1:opt_sys.n_surf-1);

sig=y.*(n./n_s).*(n-n_s).*(i+u_s)/(2*n_s(opt_sys.n_surf-1)*u_s(opt_sys.n_surf-1));

SphT=sig.*i.^2;

Ptz=-(efl*opt_sys.u_0)^2*n_s(opt_sys.n_surf-1)/2*((n_s-n).*opt_sys.c(1:opt_sys.n_surf-1)./(n.*n_s));

ComaS=sig.*i_b.*i;

AstTS=sig.*i_b.^2;

u_b_s=nuh_h(1,2:opt_sys.n_surf)'./n_s;
sig_b=y_b.*(n./n_s).*(n-n_s).*(i_b+u_b_s)/(2*n_s(opt_sys.n_surf-1)*u_s(opt_sys.n_surf-1));

Dist=sig_b.*i_b.*i + (efl*opt_sys.u_0)/2*(u_b_s.^2-u_b.^2);