function [SphT, ComaS, AstTS, Ptz, Dist]=aber3_data_opt_sys(opt_sys, flag_xy, i_c)
%ABER3_DATA_OPT_SYS 3rd order aberration contributions
%	ABER3_DATA_OPT_SYS returns the 3rd order aberration contributions of each surface of an optical system
%		(also with non-spherical surfaces) 
%		(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 (OPT_SYS.N_SURF-element array of the curvatures of the surfaces including the detector plane)
%		OPT_SYS.D (OPT_SYS.N_SURF-element array of the distances between the surfaces; 
%			the first distance D(1) is the distance from the obejct point 
%			to the first refracting surface
%		OPT_SYS.N (OPT_SYS.N_SURF-element array of the refraction indices of the media)
%		OPT_SYS.L_0 (height of the principal ray at the first (non-refracting) surface)
%		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, PAR_TRACE_OPT_SYS, EFL, BFL
%
%   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=[opt_sys.n(1,i_c)*opt_sys.h_EP/opt_sys.d_EP; opt_sys.h_EP/opt_sys.d_EP*opt_sys.d(1)];
[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.l_0/opt_sys.d_EP; opt_sys.l_0*(1-opt_sys.d(1)/opt_sys.d_EP)];
[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   

m=nuh_r(1,1)/nuh_r(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)';
w=nuh_h(1,1:opt_sys.n_surf-1)'./n;
i_b=w+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=-(m*opt_sys.l_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 + (m*opt_sys.l_0)/2*(u_b_s.^2-u_b.^2);