%--------------------------------------------------------------------------
% [OF] = FT_zernike_optical_features(A,R,n,m)
%
% A function for Matlab which returns the optical features (OF) structure
% for a given Zernike polynomial.  This structure contains information
% about the zernike polynomial in terms of the optical feature it
% represents and, in some cases, converts the Zernike amplitude into an
% equivalent measurement for this feature.  For example, the polynomial Z20
% refers to the optical feature curvature.  This function converts the
% amplitude for Z20 into the equivalent radius of curvature.  For 
% polynomials with n>2 the measurement is just the Zernike amplitude.
%
% A:    Amplitude of Zernike polynomial [nm]
% R:    Radius of Zernike polynomial, or equivalent surface [m]
% n,m:  Indices for the Zernike polynomial
% 
% Speacial case for astigmatism: Can get maximum and minimum values for
% radii of curvature of a surface, i.e. defining the astigmatism, by
% amplitude values for Z2-2, Z20 and Z22.  Need all three of these
% amplitudes.
% i.e. OF = FT_zernike_optical_features([1,1000,1],0.17,[2,2,2],[-2,0,2])
%
% OF:   Returned 'optical feature' structure.
%       OF.name:        Name of the feature.
%       OF.measurement: Description of OF.value (i.e. Rc for Z20)
%       OF.units:       Units of OF.value
%       OF.value:       Value of OF.measurement
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    05.05.2012
%--------------------------------------------------------------------------
%

function [OF] = FT_zernike_optical_features(A,R,n,m)

    baseid = 'zernike_optical_features';
    
    % Error messages
    % A, n and m must be the same length
    if (length(A)~=length(n) || length(n)~=length(m))
        result='Invalid input: A, n and m must be the same length';
        msgid=[baseid,':checkarguments']
        error(msgid,result);
    end
    
    if length(A)==1
        % Indices must be integers
        if (round(n)~=n || round(m)~=m)
            result='Invalid zernike indices: n and m must be integers';
            msgid=[baseid,':checkarguments']
            error(msgid,result);
        end
        % n must be positive or 0
        if (n<0)
            result='Invalid zernike indices: n must be 0 or positive';
            msgid=[baseid,':checkarguments']
            error(msgid,result);
        end
        % Integer difference between n and m must be even (otherwise polynomial is 
        % 0, but we don't want to include these).
        if (abs(m)>n || round((m-n)/2)~=(m-n)/2)
            result='Invalid zernike indices: m must be between -n and +n, in steps of 2';
            msgid=[baseid,':checkarguments']
            error(msgid,result);
        end
    elseif length(A)==3
        % Only valid for astigmatism (i.e. n must be 2)
        if (n(1)~=2 || n(2)~=2 || n(3)~=2)
            result='Invalid zernike indices: n should be [2,2,2]';
            msgid=[baseid,':checkarguments']
            error(msgid,result)
        end
        % Astigmatism needs m = -2,0 and 2;
        if ((m(1)~=-2 && m(2)~=-2 && m(3)~=-2) || (m(1)~=0 && m(2)~=0 && m(3)~=0) || (m(1)~=2 && m(2)~=2 && m(3)~=2))
            result='Invalid zernike indices: m should be [-2,0,2] (any order)';
            msgid=[baseid,':checkarguments']
            error(msgid,result)
        end
    else
        % Invalid amp size
        result='Invalid input: Valid size(A) is 1 or 3 (for astigmatism)';
        msgid=[baseid,':checkarguments']
        error(msgid,result)
    end
              
    % Amplitude of Zernikes is in nm
    scaling = 1e-9;
    amp = A*scaling;
    
    % Default values
    OF.name = sprintf('Z%g%g',n,m);
    OF.measurement = 'Zernike amplitude';
    OF.units = 'nm';
    OF.value = A;
    
    % Offset
    if n(1)==0 && m(1)==0
        
        OF.name = 'offset';
        OF.measurement = 'amplitude';
        OF.units = 'nm';
        OF.value = A;
        
    % Tilt
    elseif n(1)==1 && abs(m(1))==1
        
        OF.name = 'tilt';
        % Calculate tilt angle from Zernike amplitude
        OF.value = atan(amp/R);
        OF.units = 'rad';
        
        % Direction of tilt
        if m(1)==-1
            OF.measurement = 'ybeta';
        else
            OF.measurement = 'xbeta';
        end
    
    % Curvature
    elseif n(1)==2 
        
        if length(amp)==1            
            if m==0        
                OF.name = 'curvature';
        
                % Calculate Rc from Zernike amplitude (only valid for amp<R/2).
                % Otherwise OF returns defalt values for Zernike amplitude.
                if (amp<R/2)
                    OF.measurement = 'Rc';
                    OF.units = 'm';
                    OF.value = ((4*amp^2)+R^2)/(4*amp);
                end
            elseif abs(m)==2
                OF.name = 'astigmatism';    
            end
            
        elseif length(amp)==3    
            OF.name = 'astigmatism';
        
            i1 = find(m==-2);
            i2 = find(m==0);
            i3 = find(m==2);

            % Amplitude of curvature (Z20) and astigmatism (Z22) terms
            A20 = amp(i2);
            A22 = sqrt(amp(i1)^2+amp(i3)^2);
        
            % Calculate Rc from Zernike amplitudes if we are within limits
            % of this approximation.
            if ((2*A20-A22)<R)        
                OF.measurement = 'Rcs (maximum and minimum)';
                OF.units = 'm';
            
                % Calculate minimum and maximum radii of cutvature 
                Rc_min = ((2*A20+A22)^2+R^2)/(2*(2*A20+A22));
                Rc_max = ((2*A20-A22)^2+R^2)/(2*(2*A20-A22));
            
                OF.value = [Rc_max, Rc_min];            
            else 
                OF.value = A22/scaling;
            end        
        end
    end           

end

