Code covered by the BSD License  

Highlights from
Zernike calculator

image thumbnail

Zernike calculator

by

 

07 Sep 2010 (Updated )

Graphical calculator showing Zernike polynomials for a variety of aperture shapes.

zernikes.m
% Zernike polynomial calculator
% Christina Dunn
% 7 Sept 2010
% christina.r.dunn@gmail.com

% This calculator displays Zernike polynomials for several types of
% apertures, including circular, annular, rectangular, hexagonal and
% elliptical apertures.

% The circular Zernikes displayed are the FRINGE set, as defined in the
% Code V reference Manual, R. C. Juergens, ed. (Optical Research
% Associates, Pasadena, CA, 1998), Vol 1, version 8.30.
% The annular Zernikes are from V. N. Mahajan, "Zernike annular polynomials
% for imaging systems with annular pupils," J. Opt. Soc. Am., Vol. 71, No.
% 1, pg 75-85 (1981). All other Zernikes are from V. N. Mahajan and G. M. Dai,
% "Orthonormal polynomials in wavefront analysis: analytical solution,"
% Vol. 24, No. 9, pg 2994-3016 (Sept. 2007). 

% This calculator makes use of the POLAR3D function written by J M De Freitas. 
% POLAR3D: A 3-Dimensional Polar Plot Function 
% in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, 
% Dorchester DT2 8XJ. UK. 2 June 2005. 

%%
function varargout = zernikes(varargin)
% ZERNIKES M-file for zernikes.fig
%      ZERNIKES, by itself, creates a new ZERNIKES or raises the existing
%      singleton*.
%
%      H = ZERNIKES returns the handle to a new ZERNIKES or the handle to
%      the existing singleton*.
%
%      ZERNIKES('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ZERNIKES.M with the given input arguments.
%
%      ZERNIKES('Property','Value',...) creates a new ZERNIKES or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before zernikes_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to zernikes_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help zernikes

% Last Modified by GUIDE v2.5 11-Aug-2010 11:25:57

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @zernikes_OpeningFcn, ...
                   'gui_OutputFcn',  @zernikes_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


%% --- Executes just before zernikes is made visible.
function zernikes_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to zernikes (see VARARGIN)

warning off all

% Choose default command line output for zernikes
handles.output = hObject;

% Add tool bar to figure
set(hObject, 'toolbar', 'figure');

% Colors for tabs
handles.inactiveColor = [0.9 0.9 0.9];
handles.activeColor = [0.80 0.80 0.95];

% Plot style
% Options are 'SURFACE', 'MESH', 'FRINGES', 'PROFILE'
handles.plotStyle = 'SURFACE';
set(handles.plotModePanel,'SelectionChangeFcn',@plotModePanel_SelectionChangeFcn);

% The application opens in FRINGE Zernike mode
% Possible modes are:
% 'FRINGE','ANNULAR', 'RECT', 'HEX', 'ELLIPSE'
set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.activeColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);
handles.currentZernikeMode = 'FRINGE';

% No special inputs for the FRINGE set, so these are set to be invisible
set(handles.obscurationSlider, 'Visible', 'off');
set(handles.obscurationLabel, 'Visible', 'off');
set(handles.obscurationMin, 'Visible', 'off');
set(handles.obscurationMax, 'Visible', 'off');
set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'off');
set(handles.rectHalfWidthBox, 'Visible', 'off');
set(handles.rectHalfWidthSlider, 'Visible', 'off');
set(handles.rectHalfWidthMin, 'Visible', 'off');
set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'off');
set(handles.ellipseAspectRatioLabel, 'Visible', 'off');
set(handles.ellipseAspectSlider, 'Visible', 'off');
set(handles.ellipseAspectMin, 'Visible', 'off');
set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'on');

% Default values for apertures
handles.resolution = 100;
handles.obscuration = 0.5;
handles.rectHalfWidth = 0.5;
handles.ellipseAspectRatio = 0.5;

% Coefficients of all the Zernike terms start out as zero
handles.table = zeros(5,6);
handles.table(:,1) = [1:5];
handles.table(:,3) = [6:10];
handles.table(:,5) = [11:15];
set(handles.uitable, 'data', handles.table);
handles.terms = zeros(1,15);

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes zernikes wait for user response (see UIRESUME)
% uiwait(handles.figure1);


%% --- Outputs from this function are returned to the command line.
function varargout = zernikes_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


%% --- Executes on button press in FRINGEZernikeButton.
function FRINGEZernikeButton_Callback(hObject, eventdata, handles)
% hObject    handle to FRINGEZernikeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.currentZernikeMode = 'FRINGE';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.activeColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

% No special inputs for the FRINGE set, so these are set to be invisible
set(handles.obscurationSlider, 'Visible', 'off');
set(handles.obscurationLabel, 'Visible', 'off');
set(handles.obscurationMin, 'Visible', 'off');
set(handles.obscurationMax, 'Visible', 'off');
set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'off');
set(handles.rectHalfWidthBox, 'Visible', 'off');
set(handles.rectHalfWidthSlider, 'Visible', 'off');
set(handles.rectHalfWidthMin, 'Visible', 'off');
set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'off');
set(handles.ellipseAspectRatioLabel, 'Visible', 'off');
set(handles.ellipseAspectSlider, 'Visible', 'off');
set(handles.ellipseAspectMin, 'Visible', 'off');
set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'on');

guidata(hObject, handles);

%% --- Executes on button press in annularZernikeButton.
function annularZernikeButton_Callback(hObject, eventdata, handles)
% hObject    handle to annularZernikeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.currentZernikeMode = 'ANNULAR';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.activeColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.obscurationSlider, 'Visible', 'on');
set(handles.obscurationLabel, 'Visible', 'on');
set(handles.obscurationMin, 'Visible', 'on');
set(handles.obscurationMax, 'Visible', 'on');
set(handles.obscurationRatioBox, 'Visible', 'on');

set(handles.rectHalfWidthLabel, 'Visible', 'off');
set(handles.rectHalfWidthBox, 'Visible', 'off');
set(handles.rectHalfWidthSlider, 'Visible', 'off');
set(handles.rectHalfWidthMin, 'Visible', 'off');
set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'off');
set(handles.ellipseAspectRatioLabel, 'Visible', 'off');
set(handles.ellipseAspectSlider, 'Visible', 'off');
set(handles.ellipseAspectMin, 'Visible', 'off');
set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on button press in rectZernikeButton.
% The Rectangular Zernikes tab sets the GUI to draw Zernikes with a 
% rectangular aperture.
function rectZernikeButton_Callback(hObject, ~, handles)
% hObject    handle to rectZernikeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.currentZernikeMode = 'RECT';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.activeColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.obscurationSlider, 'Visible', 'off');
set(handles.obscurationLabel, 'Visible', 'off');
set(handles.obscurationMin, 'Visible', 'off');
set(handles.obscurationMax, 'Visible', 'off');
set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'on');
set(handles.rectHalfWidthBox, 'Visible', 'on');
set(handles.rectHalfWidthSlider, 'Visible', 'on');
set(handles.rectHalfWidthMin, 'Visible', 'on');
set(handles.rectHalfWidthMax, 'Visible', 'on');

set(handles.ellipseAspectRatioBox, 'Visible', 'off');
set(handles.ellipseAspectRatioLabel, 'Visible', 'off');
set(handles.ellipseAspectSlider, 'Visible', 'off');
set(handles.ellipseAspectMin, 'Visible', 'off');
set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on button press in hexZernikeButton.
function hexZernikeButton_Callback(hObject, eventdata, handles)
% hObject    handle to hexZernikeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.currentZernikeMode = 'HEX';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.activeColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.inactiveColor);

set(handles.obscurationSlider, 'Visible', 'off');
set(handles.obscurationLabel, 'Visible', 'off');
set(handles.obscurationMin, 'Visible', 'off');
set(handles.obscurationMax, 'Visible', 'off');
set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'off');
set(handles.rectHalfWidthBox, 'Visible', 'off');
set(handles.rectHalfWidthSlider, 'Visible', 'off');
set(handles.rectHalfWidthMin, 'Visible', 'off');
set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'off');
set(handles.ellipseAspectRatioLabel, 'Visible', 'off');
set(handles.ellipseAspectSlider, 'Visible', 'off');
set(handles.ellipseAspectMin, 'Visible', 'off');
set(handles.ellipseAspectMax, 'Visible', 'off');

set(handles.noVariablesText, 'Visible', 'on');

guidata(hObject, handles);

%% --- Executes on button press in ellipticalZernikeButton.
function ellipticalZernikeButton_Callback(hObject, eventdata, handles)
% hObject    handle to ellipticalZernikeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles.currentZernikeMode = 'ELLIPSE';

set(handles.FRINGEZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.annularZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.rectZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.hexZernikeButton, 'BackgroundColor', handles.inactiveColor);
set(handles.ellipticalZernikeButton, 'BackgroundColor', handles.activeColor);

set(handles.obscurationSlider, 'Visible', 'off');
set(handles.obscurationLabel, 'Visible', 'off');
set(handles.obscurationMin, 'Visible', 'off');
set(handles.obscurationMax, 'Visible', 'off');
set(handles.obscurationRatioBox, 'Visible', 'off');

set(handles.rectHalfWidthLabel, 'Visible', 'off');
set(handles.rectHalfWidthBox, 'Visible', 'off');
set(handles.rectHalfWidthSlider, 'Visible', 'off');
set(handles.rectHalfWidthMin, 'Visible', 'off');
set(handles.rectHalfWidthMax, 'Visible', 'off');

set(handles.ellipseAspectRatioBox, 'Visible', 'on');
set(handles.ellipseAspectRatioLabel, 'Visible', 'on');
set(handles.ellipseAspectSlider, 'Visible', 'on');
set(handles.ellipseAspectMin, 'Visible', 'on');
set(handles.ellipseAspectMax, 'Visible', 'on');

set(handles.noVariablesText, 'Visible', 'off');

guidata(hObject, handles);

%% --- Executes on slider movement.
% Sets the obscuration ratio for the annular Zernike
function obscurationSlider_Callback(hObject, eventdata, handles)
% hObject    handle to obscurationSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.obscuration = get(hObject, 'Value');

%set the text box
set(handles.obscurationRatioBox, 'String', handles.obscuration);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function obscurationSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to obscurationSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


%%
function rectHalfWidthBox_Callback(hObject, eventdata, handles)
% hObject    handle to rectHalfWidthBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of rectHalfWidthBox as text
%        str2double(get(hObject,'String')) returns contents of rectHalfWidthBox as a double

a = str2double(get(hObject, 'String'));

if a <= 0.9 && a >= 0.1
    handles.rectHalfWidth = a;
    set(handles.rectHalfWidthSlider, 'Value', a);
else
    errordlg('The rectangle half width must be between 0 and 1.')
end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function rectHalfWidthBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to rectHalfWidthBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%%
function ellipseAspectRatioBox_Callback(hObject, eventdata, handles)
% hObject    handle to ellipseAspectRatioBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ellipseAspectRatioBox as text
%        str2double(get(hObject,'String')) returns contents of ellipseAspectRatioBox as a double
a = str2double(get(hObject, 'String'));

if a <= 1 && a >= 0
    handles.ellipseAspectRatio = a;
    set(handles.ellipseAspectSlider, 'Value', a);
else
    errordlg('The aspect ratio must be between 0 and 1.')
end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function ellipseAspectRatioBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ellipseAspectRatioBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%%
function obscurationRatioBox_Callback(hObject, eventdata, handles)
% hObject    handle to obscurationRatioBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of obscurationRatioBox as text
%        str2double(get(hObject,'String')) returns contents of obscurationRatioBox as a double

a = str2double(get(hObject, 'String'));

if a <= 1 && a >= 0
    handles.obscuration = a;
    set(handles.obscurationSlider, 'Value', a);
else
    errordlg('The obscuration ratio must be between 0 and 1.')
end

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function obscurationRatioBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to obscurationRatioBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


%% --- Executes on slider movement.
function rectHalfWidthSlider_Callback(hObject, eventdata, handles)
% hObject    handle to rectHalfWidthSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider


handles.rectHalfWidth = get(hObject, 'Value');

%set the text box
set(handles.rectHalfWidthBox, 'String', handles.rectHalfWidth);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function rectHalfWidthSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to rectHalfWidthSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


%% --- Executes on slider movement.
function ellipseAspectSlider_Callback(hObject, eventdata, handles)
% hObject    handle to ellipseAspectSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.ellipseAspectRatio = get(hObject, 'Value');

%set the text box
set(handles.ellipseAspectRatioBox, 'String', handles.ellipseAspectRatio);

guidata(hObject, handles);

%% --- Executes during object creation, after setting all properties.
function ellipseAspectSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ellipseAspectSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


% --- Executes on button press in displayButton.
function displayButton_Callback(hObject, eventdata, handles)
% hObject    handle to displayButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
cla reset

switch handles.currentZernikeMode
    case 'FRINGE'
        drawFringeZernikes(hObject, handles);
    case 'ANNULAR'
        drawAnnularZernikes(hObject, handles);
    case 'RECT'
        drawRectZernikes(hObject, handles);
    case 'HEX'
        drawHexZernikes(hObject, handles);
    case 'ELLIPSE'
        drawEllipseZernikes(hObject, handles);
end

function drawFringeZernikes(hObject, handles)

resolution = handles.resolution;

resolution = handles.resolution;
r = linspace(1, 0, resolution)';
a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

F = zeros(resolution, resolution, 15);

F(:,:,1) = ones(resolution, resolution);
F(:,:,2) = r * cos(a);
F(:,:,3) = r * sin(a);
F(:,:,4) = (2*r.^2 - ones(size(r))) * ones(size(a));
F(:,:,5) = r.^2 * cos(2*a);
F(:,:,6) = r.^2 * sin(2*a);
F(:,:,7) = (3*r.^3 - 2*r) * cos(a);
F(:,:,8) = (3*r.^3 - 2*r) * sin(a);
F(:,:,9) = (6*r.^4 - 6*r.^2 + 1) * ones(size(a));
F(:,:,10) = r.^3 * cos(3*a);
F(:,:,11) = r.^3 * sin(3*a);
F(:,:,12) = (4*r.^4 -3*r.^2) * cos(2*a);
F(:,:,13) = (4*r.^4 -3*r.^2) * sin(2*a);
F(:,:,14) = (10*r.^5 - 12*r.^3 + 3*r) * cos(a);
F(:,:,15) = (10*r.^5 - 12*r.^3 + 3*r) * sin(a);

result = zeros(resolution, resolution);
for n = 1:15
    result = result + handles.terms(n) * F(:,:,n);
end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1);

plotZernike(x, y, z, handles);

function drawAnnularZernikes(hObject, handles)

e = handles.obscuration;
resolution = handles.resolution;

r = linspace(1, e, resolution)';
a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

A = zeros(resolution, resolution, 15);

A(:,:,1) = ones(resolution, resolution);
A(:,:,2) = 2 * (r/sqrt(1 + e^2)) * cos(a);
A(:,:,3) = 2 * (r/sqrt(1 + e^2)) * sin(a);
A(:,:,4) = sqrt(3) *((2*r.^2 - 1 - e^2)/(1 - e^2))*ones(size(a));
A(:,:,5) = sqrt(6) * (r.^2 / sqrt(1 + e^2 + e^4)) * sin(2*a);
A(:,:,6) = sqrt(6) * (r.^2 / sqrt(1 + e^2 + e^4)) * cos(2*a);
A(:,:,7) = sqrt(8) * ((3*(1+e^2)*r.^3 - 2*(1+ e^2 + e^4)*r)/ ...
    ((1-e^2)*sqrt((1+4*e^2+e^4)))) * sin(a);
A(:,:,8) = sqrt(8) * ((3*(1+e^2)*r.^3 - 2*(1+ e^2 + e^4)*r)/...
    ((1-e^2)*sqrt((1+4*e^2+e^4)))) * cos(a);
A(:,:,9) = sqrt(8) * (r.^3 / sqrt(1 + e^2 + e^4 + e^6)) * sin(3*a);
A(:,:,10) = sqrt(8) * (r.^3 / sqrt(1 + e^2 + e^4 + e^6)) * cos(3*a);
A(:,:,11) = sqrt(5) * (6*r.^4 - 6*(1+e^2)*r.^2 + 1+ 4*e^2 + e^4)/(1-e^2)^2 ...
    * ones(size(a));
A(:,:,12) = sqrt(10) * (4 * r.^2 - 3*((1-e^8)/(1-e^6))*r.^2) / ...
    sqrt((1-e^2)^-1 * (16*(1-e^10) - 15*(1-e^8)^2/(1-e^6))) * cos(2*a);
A(:,:,13) = sqrt(10) * (4 * r.^2 - 3*((1-e^8)/(1-e^6))*r.^2) / ...
    sqrt((1-e^2)^-1 * (16*(1-e^10) - 15*(1-e^8)^2/(1-e^6))) * sin(2*a);
A(:,:,14) = sqrt(10) * r.^4/sqrt(1 + e^2 + e^4 + e^6 + e^8) * cos(4*a);
A(:,:,15) = sqrt(10) * r.^4/sqrt(1 + e^2 + e^4 + e^6 + e^8) * sin(4*a);

result = zeros(resolution, resolution);
for n = 1:15
    result = result + handles.terms(n) * A(:,:,n);
end

[x, y, z] = polar3d(result, 0, 2*pi, e, 1, 1);

plotZernike(x, y, z, handles);

%%
function drawRectZernikes(hObject, handles)

a = handles.rectHalfWidth;
resolution = handles.resolution;
xRes = linspace(-a, a, resolution);
yRes = linspace(-sqrt(1-a^2), sqrt(1-a^2), resolution);

[x, y] = meshgrid(xRes, yRes);
r = sqrt(x.^2 + y.^2);

u = sqrt(9 - 36*a^2 + 103*a^4 - 134*a^6 + 67*a^8);
v = sqrt(49 - 196*a^2 + 330*a^4 - 268*a^6 + 134*a^8);
t = 1/(128*v*a^4*(1-a^2)^2);
n = 9 - 45*a^2 + 139*a^4 - 237*a^6 + 210*a^8 - 67*a^10;

R = zeros(resolution, resolution,15);
R(:,:,1) = ones(size(x));
R(:,:,2) = (sqrt(3)/a) * x;
R(:,:,3) = sqrt(3/(1-a^2)) * y;
R(:,:,4) = (sqrt(5)/(2 * sqrt(1 - 2*a^2 + 2*a^4)))*(3*(x.^2 + y.^2) - 1);
R(:,:,5) = (3/(a*sqrt(1-a^2))) * x .* y;
R(:,:,6) = (sqrt(5)*(2*a^2*(1-a^2)*sqrt(1-2*a^2+2*a^4))) * ...
    (3*(1-a^2)^2*x.^2 - 3*a^4*y.^2 - a^2*(1 - 3*a^2 + 2*a^4));
R(:,:,7) = (sqrt(21)/(2*sqrt(27-81*a^2+116*a^4 - 62*a^6)))*(15*(x.^2 + y.^2) - 9 + 4*a^2) .* y;
R(:,:,8) = (sqrt(21)/(2*a*sqrt(35-70*a^2+62*a^4)))*(15*r.^2 - 5 - 4*a^2)*x;
R(:,:,9) = (sqrt(5)*sqrt((27-54*a^2+62*a^4)/(1-a^2))/(2*a^2*(27-81*a^2+116*a^4-62*a^6))) ...
    * (27*(1-a^2)^2*x.^2 - 35*a^4*y.^2 - a^2*(9 - 39*a^2 + 30*a^4)) .* y;
R(:,:,10) = (sqrt(5)/(2*a^3*(1-a^2)*sqrt(35-70*a^2+62*a^4)))*...
    (35*(1-a^2)^2*x.^2 - 27*a^4*y.^2 - a^2*(21 - 51*a^2 + 30*a^4))*x;
R(:,:,11) = (1/(8*u))*(315*r.^4 - 30*(7+2*a^2)*x.^2-30*(9-2*a^2)*y.^2 + 27 + 16*a^2 - 16*a^4);
R(:,:,12) = (3*u/(8*a^2*v*n))*(35*(1-a^2)^2*(18-36*a^2+67*a^4)*x.^4 + ...
    630*(1-2*a^2)*(1-2*a^2+2*a^4)*x.^2*y.^2 ...
    - 35*a^4*(49-98*a^2+67*a^4)*y.^4 - 30*(1-a^2)*(7-10*a^2-12*a^4+75*a^6-67*a^8)*x.^2 ...
    -30*a^2*(7-77*a^2+189*a^4-193*a^6+67*a^8)*y.^2 + ...
    a^2*(1-a^2)*(1-2*a^2)*(70-233*a^2+233*a^4));
R(:,:,13) = (sqrt(21)/(2*a*sqrt(1-3*a^2+4*a^4-2*a^6)))*(5*r.^2-3)* x .* y;
R(:,:,14) = 16*t*(735*(1-a^2)^4*x.^4 - 540*a^4*(1-a^2)^2*x.^2 .* y.^2 + 735*a^8*y.^4 ...
    - 90*a^2*(1-a^2)^3*(7-9*a^2)*x.^2 ...
    + 90 * a^6 *(1-a^2)*(2-9*a^2)*y.^2 + 3*a^4*(1-a^2)^2*(21-62*a^2+62*a^4));
R(:,:,15) = (sqrt(21)/(2*a^3*(1-a^2)*sqrt(1-3*a^2+4*a^4-2*a^6))) * ...
    (5*(1-a^2)^2*x.^2 - 5*a^4*y.^2 - a^2*(3 - 9*a^2 + 6*a^4)) .* x .* y;

z = zeros(resolution, resolution);
for n = 1:15
    z = z + handles.terms(n) * R(:,:,n);
end

plotZernike(x, y, z, handles);

guidata(hObject, handles);

%%
function drawHexZernikes(hObject, handles)

resolution = 2*handles.resolution;
r = linspace(1, 0, resolution)';
a = linspace(0, 2 * pi + 2*pi/resolution, resolution);

H = zeros(resolution, resolution,15);
H(:,:,1) = ones(resolution); 
H(:,:,2) = 2 * sqrt(6/5) * r * cos(a);
H(:,:,3) = 2 * sqrt(6/5) * r * sin(a);
H(:,:,4) = sqrt(5/43) * (-5 + 12 * r.^2) * ones(size(a));
H(:,:,5) = 2 * sqrt(15/7) * r.^2 * sin(2*a);
H(:,:,6) = 2 * sqrt(15/7) * r.^2 * cos(2*a);
H(:,:,7) = 4 * sqrt(42/3685) * (-14 * r + 25 * r.^3) * sin(a);
H(:,:,8) = 4 * sqrt(42/3685) * (-14 * r + 25 * r.^3) * cos(a);
H(:,:,9) = (4*sqrt(10)/3) * r.^3 * sin(3*a);
H(:,:,10) = 4 * sqrt(70/103) * r.^3 * cos(3*a);
H(:,:,11) = (3/sqrt(1072205)) * (737 - 5140 * r.^2 + 6020 * r.^4) * ones(size(a));
H(:,:,12) = (30 / sqrt(492583)) * (-249 * r.^2 + 392 * r.^4) * cos(2*a);
H(:,:,13) = (30 / sqrt(492583)) * (-249 * r.^2 + 392 * r.^4) * sin(2*a);
H(:,:,14) = (10/3) * sqrt(7/99258181) * (10 * (297 - 598 * r.^2) .* r.^2 * cos(2*a) + 5413 * r.^4 * cos(4*a));
H(:,:,15) = (10/3) * sqrt(7/99258181) * (-10 * (297 - 598 * r.^2) .* r.^2  * sin(2*a) + 5413 * r.^4 * sin(4*a));

result = zeros(resolution, resolution);
for n = 1:15
    result = result + handles.terms(n) * squeeze(H(:,:,n));
end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1);

xHex = [1/2; 1; 1/2; -1/2; -1; -1/2; 1/2];
yHex = [sqrt(3)/2; 0 ; -sqrt(3)/2; -sqrt(3)/2; 0; sqrt(3)/2; sqrt(3)/2];
in = inpolygon(x, y, xHex, yHex);
x(~in) = NaN;
y(~in) = NaN;
z(~in) = NaN;

plotZernike(x, y, z, handles);

%%
function drawEllipseZernikes(hObject, handles)

resolution = 2*handles.resolution;
r = linspace(1, 0, resolution)';
a = linspace(0, 2 * pi + 2*pi/resolution, resolution);
b = handles.ellipseAspectRatio;
alpha = sqrt(45 - 60*b^2 + 94*b^4 - 60*b^6 + 45*b^8);
beta = sqrt(1575 - 4800*b^2 + 12020*b^4 - 17280*b^6 - 21066*b^8 - 17280*b^10 ...
    + 12020*b^12 - 4800*b^14 + 1575*b^16);
gamma = sqrt(35 - 60*b^2 + 114*b^4 - 60*b^6 + 35*b^8);
delta = sqrt(5 - 6*b^2 + 5*b^4);

E = zeros(resolution, resolution,15);
E(:,:,1) = ones(resolution); 
E(:,:,2) = 2 * r * cos(a);
E(:,:,3) = (2 * r * sin(a))/b;
E(:,:,4) = (sqrt(3) / sqrt(3 - 2*b^2 + 3*b^4)) * (-1 - b^2 + 4*r.^2) * ones(size(r))';
E(:,:,5) = (sqrt(6)/b) * r.^2 * sin(2*a);
E(:,:,6) = (1/(2*b)) * r.^2 * sin(2*a);
E(:,:,7) = (4/(b*sqrt(5-6*b^2+9*b^4)))*(-1*(1+3*b^2)*r + 6*r.^3) * sin(a);
E(:,:,8) = (4/sqrt(9-6*b^2+5*b^4))*(-1*(3+b^2)*r + 6*r.^3)*cos(a);
E(:,:,9) = (1/(b^3*sqrt(5-6*b^2+9*b^4)))* ... 
    (3*(4*b^2*(1-b^2)*r -(5-2*b^2-3*b^4)*r.^3)*sin(a) + ...
    (5-6*b^2+9*b^4)*r.^3 * sin(3*a));
E(:,:,10) = (1/(b^2*sqrt(9-6*b^2+5*b^4))*(3*(4*b^2*(1-b^2)*r-(3+2*b^2-5*b^4)*r.^3)*...
    cos(a) + (9-6*b^2+5*b^4)*r.^3*cos(3*a)));
E(:,:,11) = sqrt(5)*(3 + 2*b^2 + 3*b^4 - 24*(1 + b^2)*r.^2 *ones(size(a)) + ...
    48*r.^4*ones(size(a)) - 12*(1-b^2)*r.^2 * cos(2*a)) * alpha;
E(:,:,12) = (sqrt(10)*alpha/(gamma*b^2))*(-3*r.^2+4*r.^4)*cos(2*a) + (sqrt(5/2)/(2*b^2*beta)) * ...
    (-12*b^2*(5-2*b^2+2*b^6-5*b^8) + 6*(15+125*b^2-194*b^4 + 194*b^6 - 125*b^8 - 15*b^10)*r.^2*ones(size(a)) + ...
    240*(-3 + 2*b^2 - 2*b^6 + 3*b^8)*r.^4*ones(size(a)) ...
    + 6*(75-155*b^2 + 174*b^4 - 134*b^6 + 55*b^8 - 15*b^10)*r.^2 * cos(2*a));
E(:,:,13) = (sqrt(10)/b)*(1/delta)*(-3*(1+b^2)*r.^2 + 8*r.^4)*sin(2*a);
E(:,:,14) = (sqrt(10)/(8*b^4*gamma))*(3*(1-b^2)^2 * (8*b^4-40*b^2*(1+b^2)*r.^2+5*(7+10*b^2+7*b^4)*r.^4)*ones(size(a)) ...
    +4*(6*b^2*(5-7*b^2+7*b^4-5*b^6)-5*(7-6*b^2+6*b^6-7*b^8)*r.^2).*r.^2*cos(2*a) + (35 - 60*b^2 + 114*b^4 ...
    -60*b^6 + 35*b^8)*r.^4*cos(4*a));
E(:,:,15) = (sqrt(10)/b^3)*(1/delta)*((6*b^2*(1-b^2)-5*(1-b^4)*r.^2).*r.^2*sin(2*a) + ...
    ((5 - 6*b^2+5*b^4)/2)*r.^4*sin(4*a));

result = zeros(resolution, resolution);
for n = 1:15
    result = result + handles.terms(n) * squeeze(E(:,:,n));
end

[x, y, z] = polar3d(result, 0, 2*pi, 0, 1, 1);
theta = linspace(0, 2*pi);
xEllipse = cos(theta);
yEllipse = b * sin(theta);

in = inpolygon(x, y, xEllipse, yEllipse);
x(~in) = NaN;
y(~in) = NaN;
z(~in) = NaN;

plotZernike(x, y, z, handles);


% --- Executes when entered data in editable cell(s) in uitable.
function uitable_CellEditCallback(hObject, eventdata, handles)
% hObject    handle to uitable (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) edited
%	PreviousData: previous data for the cell(s) edited
%	EditData: string(s) entered by the user
%	NewData: EditData or its converted form set on the Data property. Empty if Data was not changed
%	Error: error string when failed to convert EditData to appropriate value for Data
% handles    structure with handles and user data (see GUIDATA)
handles.table(eventdata.Indices(1), eventdata.Indices(2)) = eventdata.NewData;
handles.terms = [handles.table(:, 2)' handles.table(:,4)' handles.table(:,6)'];

guidata(hObject, handles);


% --- Executes on button press in clearButton.
function clearButton_Callback(hObject, eventdata, handles)
% hObject    handle to clearButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Coefficients of all the Zernike terms start out as zero
handles.table = zeros(5,6);
handles.table(:,1) = [1:5];
handles.table(:,3) = [6:10];
handles.table(:,5) = [11:15];
set(handles.uitable, 'data', handles.table);

handles.terms = zeros(1,15);

axes(handles.profile1)
cla reset
axes(handles.profile2)
cla reset
axes(handles.profile3)
cla reset

axes(handles.display)
cla reset

set(handles.display, 'Visible', 'on'); 
set(handles.profile1, 'Visible', 'off');
set(handles.profile2, 'Visible', 'off');
set(handles.profile3, 'Visible', 'off');

guidata(hObject, handles);

function plotModePanel_SelectionChangeFcn(hObject, eventdata)
 
%retrieve GUI data, i.e. the handles structure
handles = guidata(hObject); 
 
switch get(eventdata.NewValue,'Tag')   % Get Tag of selected object
    case 'surfButton'
        handles.plotStyle = 'SURFACE';
    case 'meshButton'
        handles.plotStyle = 'MESH';
    case 'phaseButton'
        handles.plotStyle = 'FRINGES';
    case 'profileButton'
        handles.plotStyle = 'PROFILE';
end

%updates the handles structure
guidata(hObject, handles);

%%
function plotZernike(x, y, z, handles)

axes(handles.profile1)
cla reset
axes(handles.profile2)
cla reset
axes(handles.profile3)
cla reset
axes(handles.display)
cla reset

switch handles.plotStyle
    case 'SURFACE';
        set(handles.display, 'Visible', 'on'); 
        set(handles.profile1, 'Visible', 'off');
        set(handles.profile2, 'Visible', 'off');
        set(handles.profile3, 'Visible', 'off');
        
        axes(handles.display)
        rotate3d on;
        axis(gca,'normal');
        colormap('default')
        surf(x, y, z)
        lighting phong
        shading interp
        light  
        axis equal
        xlim([-1 1])
        ylim([-1 1])
        view(2)
    case 'MESH'
        set(handles.display, 'Visible', 'on'); 
        set(handles.profile1, 'Visible', 'off');
        set(handles.profile2, 'Visible', 'off');
        set(handles.profile3, 'Visible', 'off');
        
        axes(handles.display)
        rotate3d on;
        colormap('default')
        meshc(x, y, z) 
        axis equal
        xlim([-1 1])
        ylim([-1 1])
        view(2)
    case 'FRINGES'
        set(handles.display, 'Visible', 'on'); 
        set(handles.profile1, 'Visible', 'off');
        set(handles.profile2, 'Visible', 'off');
        set(handles.profile3, 'Visible', 'off');
        
        axes(handles.display)
        rotate3d off
        a = [1 1 1];
        bFringe = sin(linspace(0, pi, 50));
        b = [];
    
        for t = 1:sum(handles.terms)
            b = [b bFringe];
        end
    
        if ~isempty(b)
            map = ((a'*b)');
            colormap(map)
            surf(x, y, zeros(size(x)), z, 'LineStyle', 'none')
            view([0 90])
            axis equal
            xlim([-1 1])
            ylim([-1 1])
        end
    case 'PROFILE'
        set(handles.display, 'Visible', 'off'); 
        set(handles.profile1, 'Visible', 'on');
        set(handles.profile2, 'Visible', 'on');
        set(handles.profile3, 'Visible', 'on');
        
        axes(handles.profile1)
        cla reset
        rotate3d off;
        colormap('default')
        surf(x, y, z)
        lighting phong
        shading interp
        light  
        hold on
        zmax = max(max(z));
        plot3([-1 1], [0 0], [zmax zmax], 'r');
        plot3([0 0], [-1 1], [zmax zmax], 'g');
        view([0 90])
        
        dims = size(z);
        
        axes(handles.profile2);
        
        x(isnan(x)) = [];
        y(isnan(y)) = [];
        z(isnan(z)) = [];
                
        xprofile = linspace(-1, 1, 40)';
        yprofile = linspace(0, 0, 40)';
        zprofile = griddata(x(1:10:end,1:10:end), y(1:10:end,1:10:end), z(1:10:end,1:10:end), xprofile, yprofile)';

        plot(xprofile, zprofile, 'r');
        ylim([-zmax zmax])
         
        axes(handles.profile3)
        
        xprofile = linspace(0, 0, 40)';
        yprofile = linspace(-1, 1, 40)';
        zprofile = griddata(x(1:10:end,1:10:end), y(1:10:end,1:10:end), z(1:10:end,1:10:end), xprofile, yprofile)';

        plot(yprofile, zprofile, 'g');
        ylim([-zmax zmax])
                
end

%%
%Following code written by:
%J M De Freitas. POLAR3D: A 3-Dimensional Polar Plot Function 
%               in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, 
%               Dorchester DT2 8XJ. UK. 2 June 2005. 

function [Xout, Yout, Zout] = polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,varargin)
%
% POLAR3D   Plots a 3D polar surface.
% 
%           POLAR3D(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale) plots 
%           the profiles as a mesh plot for Zin between radii Rho_min and 
%           Rho_max for polar angles equally spaced between theta_min and theta_max.
% 
%           POLAR3D(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,plotspec)
%           plots the profiles Zin between radii Rho_min and Rho_max for 
%           polar angles between theta_min and theta_max with a plot type 
%           specification. If plotspec = 'surf' a standard Matlab surface 
%           is plotted,whereas 'mesh', 'surfc' or 'meshc' will plot mesh, 
%           surface with countour, or mesh with contour, respectively.
%           The size of the squares on the mesh or surf plots is determined 
%           by meshscale. The default plot is a mesh plot.  
%            
% 
%           [Xout,Yout,Zout] = POLAR3D(Zin,theta_min,theta_max,Rho_min,
%           Rho_max, meshscale)returns Zout values corresponding to the
%           Cartesian positions (Xout,Yout).         
% 
% SYNTAX    polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale)
%           polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,plotspec)
%           polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale,interpspec)
%           polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,meshscale, plotspec,interpspec)
%           [Xout,Yout,Zout] = polar3d(Zin,theta_min,theta_max,Rho_min,Rho_max,...)
%           
% INPUT     Zin         input magnitude profiles where each column in Zin is 
%                       assumed to represent the radial information in the
%                       plot i.e. each column represents the information 
%                       along a radial line defined by theta, where theta 
%                       varies between theta_min and theta_max.
%                       
%                       Zin is a (M x N) matrix, where M and N are not 
%                       necessarily equal. If M is not equal to N then the
%                       data are interpolated to make them equal. The final
%                       size is determined by the larger value of (M,N).
% 
%                       The N columns of Zin are assumed to be equally
%                       spaced measurements of the radial components, where
%                       Zin(1,1) corresponds to (theta_min,Rho_max) and 
%                       Zin(M,1) corresponds to (theta_min,Rho_min), and so on. 
%                       Zin(1,N) corresponds to (theta_max,Rho_max) and 
%                       Zin(M,N) corresponds to (theta_max,Rho_min). Theta 
%                       increases in the anticlockwise direction.
% 
%           theta_min   the lower value in radians of the angular range 
%                       over which the data is defined. Theta_min is a 
%                       scalar quantity.
% 
%           theta_max   the upper value in radians of the angular range 
%                       over which the data is defined. Theta_max is a 
%                       scalar quantity.
% 
%           Rho_min     the lower value of the radial component for which 
%                       the data is defined. Rho_min is a scalar quantity.
% 
%           Rho_max     the upper value of the radial component for which 
%                       the data is defined. Rho_max is a scalar quantity.
% 
%           meshscale   a scalar that determines the size of the squares 
%                       on the mesh or surf plots. If meshscale is 1, the 
%                       mesh remains unchanged relative to the input grid; 
%                       if meshscale is 2, the size of the squares is doubled,
%                       if 3.2, say, it is scaled accordingly. Moreover, if  
%                       meshscale is less than 1, e.g. 0.2, the size of the  
%                       squares is reduced into a finer mesh. The dimensions
%                       of Xout, Yout and Zout are reduced or enlarged by
%                       meshscale.
% 
%           plotspec    = 'surf'    produces a surface plot.
%                       = 'surfc'   produces a surface plot with contour.
%                       = 'mesh'    produces a mesh plot.
%                       = 'meshc'   produces a mesh plot with countour.
%                       = 'contour' produces a 2D contour plot.
%                       = 'off' disengages plot function.
%                       = 'meshl'   produces a mesh plot with labelled axes.
%
%                                   It may be necessary for you to scale
%                                   Zin when this option is used because
%                                   the DataAspectRatio is 1:1:1. This is
%                                   done by replacing 'Zin' with 'n*Zin', 
%                                   e.g. polar3d(7.12*Zin,...) where the 
%                                   scaling factor in this case is 7.12.
%
%                                   Note also that when this option is
%                                   invoked, the data does not undergo any
%                                   interpolation; the same size matrix is
%                                   returned and plotted. 'meshscale' is
%                                   ignored.
% 
%           interpspec  = 'linear'  bilinear interpolation on Zin.
%                       = 'spline'  spline interpolation on Zin.
%                       = 'nearest'  nearest neighbour interpolation on Zin.
%                       = 'cubic'  bicubic interpolation on Zin.
% 
% OUTPUT    Zout        output magnitude profiles defined by Zin at
%                       positions (Xout,Yout). 
%                       
%                       Zout is square with dimensions determined by the 
%                       maximum dimension of the input matrix Zin. The 
%                       dimensions of Zout are reduced or enlarged by meshscale.   
%
%                       If 'meshl' option is used Zout is same size as Zin. 
% 
%           Xout        output X-positions corresponding to polar positions
%                       (rho,theta). Xout is square with dimensions  
%                       determined by the maximum dimension of the input 
%                       matrix Zin. The dimensions of Xout are reduced or 
%                       enlarged by meshscale.
%
%                       If 'meshl' option is used Xout is same size as Zin
% 
%           Yout        output Y-positions corresponding to polar positions
%                       (rho,theta). Yout is square with dimensions  
%                       determined by the maximum dimension of the input 
%                       matrix Zin. The dimensions of Yout are reduced or 
%                       enlarged by meshscale.
%
%                       If 'meshl' option is used Yout is same size as Zin
%
% See also POLAR, SPHERE3D, CYL3D, POL2CART and INTERP2

%           Written by JM DeFreitas, QinetiQ Ltd, Winfrith Technology
%           Centre, Dorchester DT2 8XJ, UK. jdefreitas@qinetiq.com.
% 
%           Revision 1.0 4 April 2005.
%           Released 5 April 2005. (Beta Release).
%
%           Revision 1.1 17 April 2005 
%           1. Inroduced new check on Zin to cover case of dimensions (M x 2) 
%              or (2 x N) matrix. These are not allowed.
%           2. Changed L2 so that when meshscale > 1 and theta ranges over 
%              360 degrees the data wraps around without spaces.
%           3. Removed Xout1, Yout1 and Zout1.
%           4. Changed 'theta(j,:) = ones([(L2/n) 1])*angl(j);' to
%              'theta(j,:) = ones([1 (L2/n)])*angl(j)'; so that it is
%              compatible with Matlab6 R12.1.
%           5. Reorganised meshgrid so that interpolation now works with
%              meshscale < 1.
%           6. Changed error traps from '((p ~= 1)&(q ~= 1))' to 
%              '((p ~= 1)|(q ~= 1))' where used.
%
%           Revision 1.2 2 March 2006 
%           1. Changed the 'contourf' option to 'meshl' option. This new
%              option allows the user to plot a standard mesh plot with
%              labelled polar axes. 
%           2. The 'meshl' option plots with a DataAspectRatio of 1:1:1 and
%              as such does not automatically scale Zin. It is important for
%              the user to suitably scale Zin in this eventuality.
%           3. The angles are labelled in degrees.
%              
%
%           Full Release 26 May 2005. All Rights Reserved. 
%     
%           Terms and Conditions of Use
% 
%           1.	This function is made available to Matlab users under the 
%               terms and conditions set out in the Matlab Exchange by 
%               The Mathworks, Inc.
%           2.	Where the use of POLAR3D is to be cited, the following is recommended:
%               J M De Freitas. POLAR3D: A 3-Dimensional Polar Plot Function 
%               in Matlab. Version 1.2. QinetiQ Ltd, Winfrith Technology Centre, Winfrith, 
%               Dorchester DT2 8XJ. UK. 2 June 2005. 
%           3.	No offer of warranty is made or implied by the author and 
%               use of this work implies that the user has agreed to take full 
%               responsibility for its use.
%

if (nargin < 6)
    disp('Polar3d Error: Too few input arguments.');
    return
elseif (nargin > 8)
    disp('Polar3d Error: Too many input arguments.');
    return
end

[p,q] = size(theta_min);
if (((p ~= 1)|(q ~= 1))|~isreal(theta_min))|ischar(theta_min)
    disp('Polar3d Error: theta_min must be scalar and real.');
    return
end
[p,q] = size(theta_max);
if (((p ~= 1)|(q ~= 1))|~isreal(theta_max))|ischar(theta_max)
    disp('Polar3d Error: theta_max must be scalar and real.');
    return
end
if theta_max <= theta_min
    disp('Polar3d Error: theta_max less than or equal theta_min.');
    return
end
if abs(theta_max - theta_min) > 2*pi
    disp('Polar3d Error: range of theta greater than 2pi.');
    return
end

[p,q] = size(Rho_max);
if (((p ~= 1)|(q ~= 1))|~isreal(Rho_max))|(ischar(Rho_max)|Rho_max < 0)
    disp('Polar3d Error: Rho_max must be scalar, positive and real.');
    return
end
[p,q] = size(Rho_min);
if (((p ~= 1)|(q ~= 1))|~isreal(Rho_min))|(ischar(Rho_min)|Rho_min < 0)
    disp('Polar3d Error: Rho_min must be scalar, positive and real.');
    return
end
if Rho_max <= Rho_min
    disp('Polar3d Error: Rho_max less than or equal Rho_min.');
    return
end

[p,q] = size(meshscale);
if (((p ~= 1)|(q ~= 1))|~isreal(meshscale))|ischar(meshscale)
    disp('Polar3d Warning: mesh scale must be scalar and real.');
    meshscale = 1;
end
if (meshscale <= 0)
    disp('Polar3d Warning: mesh scale must be scalar and positive.');
    meshscale = 1;
end

% Set up default plot and interpolation specifications.
str1 = 'mesh';
str2 = 'linear';
if length(varargin) == 2
% Sort out plot and interpolation specification if both strings given.   
    str1 = [varargin{1}(:)]';
    str2 = [varargin{2}(:)]';
    g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off');
    g2 = (~isequal(str1,'meshc')&~isequal(str1,'surfc'));
    g5 = (~isequal(str1,'contour')&~isequal(str1,'meshl'));
    g3 = (~isequal(str2,'cubic')&~isequal(str2,'linear'));
    g4 = (~isequal(str2,'spline')&~isequal(str2,'nearest'));
    if (g1&g2&g5)
        disp('Polar3d Warning: Incorrect plot specification. Default to mesh plot.');
        str1 = 'mesh';
    end
    if (g3&g4)
        disp('Polar3d Warning: Incorrect interpolation specification.'); 
        disp('Default to linear interpolation.');
        str2 = 'linear';
    end  
elseif length(varargin) == 1
% Sort out plot or interpolation specification from single string input.    
    str1 = [varargin{1}(:)]';
    g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off');
    g2 = (~isequal(str1,'meshc')&~isequal(str1,'surfc'));
    g5 = (~isequal(str1,'contour')&~isequal(str1,'meshl'));
    g3 = (~isequal(str1,'cubic')&~isequal(str1,'linear'));
    g4 = (~isequal(str1,'spline')&~isequal(str1,'nearest'));    
    if (g1&g2)&(g3&g4&g5)
        disp('Polar3d Error: Incorrect plot or interpolation specification.');
        return
    elseif isequal(str1,'cubic')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'linear')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'spline')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'nearest')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'off')
        str2 = 'linear';   
    end 
end;

% Check if dimensions of input data are acceptable.
[r,c] = size(Zin');
if (r < 5)&(c < 5)
    disp('Polar3d Error: Input matrix dimensions must be greater than (4 x 4).');
    return    
end
% Check if input data has two rows or columns or less.
if (r < 3)|(c < 3)
    disp('Polar3d Error: One or more input matrix dimensions too small.');
    return    
end

% Transpose and setup input magnitude matrix.
temp = Zin';
for j = 1:c
    P(:,j) = temp(:,c-j+1); % swap columns over
end
Zin = P; 
[r,c] = size(Zin);

% Check if meshscale is compatible with dimensions of input data.
scalefactor = round(max(r,c)/meshscale);
if scalefactor < 3
    disp('Polar3d Error: mesh scale incompatible with dimensions of input data.');
    return    
end

% Set up meshgrid corresponding to larger matrix dimension of Zin
% for interpolation if required.
if ~isequal(str1,'meshl')
    n = meshscale;
    if r > c
        L = r;
        L2 = fix(L/n)*n;
        step = r/(c-1);
        [X1,Y1] = meshgrid(0:step:r,1:r);
        if n < 1
            [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
        else
            [X,Y] = meshgrid(1:n:L2,1:n:L2); 
        end
        T = interp2(X1,Y1,Zin,X,Y,str2);
    elseif c > r
        L = c;
        L2 = fix(L/n)*n;
        step = c/(r-1);
        [X1,Y1] = meshgrid(1:c,0:step:c);
        if n < 1
            [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
        else
            [X,Y] = meshgrid(1:n:L2,1:n:L2); 
        end
        T = interp2(X1,Y1,Zin,X,Y,str2);
    else
        L = r;
        L2 = fix(L/n)*n;
        [X1,Y1] = meshgrid(1:r,1:r);
        if n < 1
            [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
        else
            [X,Y] = meshgrid(1:n:L2,1:n:L2); 
        end
        T = interp2(X1,Y1,Zin,X,Y,str2);
    end

    [p,q] = size(T);
    L2 = max(p,q);

    % Set up angles
    angl = theta_min:abs(theta_max-theta_min)/(L2-1):theta_max;
    for j = 1:L2  
        theta(j,:) = ones([1 L2])*angl(j);
    end

    % Set up radial components
    Rho = Rho_min:abs(Rho_max-Rho_min)/(L2-1):Rho_max;

    % Convert to Cartesian coordinates
    for j = 1:L2  
        [Xout(j,:) Yout(j,:) Zout(j,:)] = pol2cart(theta(j,:),Rho,T(j,:));
    end
end %if str1 ~= 'meshl'

% Plot Cartesian surface
switch str1;
case 'mesh'
    colormap([0 0 0]);
    mesh(Xout,Yout,Zout);
    axis off;
    grid off;
case 'meshc'
    colormap([0 0 0]);
    meshc(Xout,Yout,Zout);
    axis off;
    grid off;
case 'surf'
    surf(Xout,Yout,Zout);
    axis off;
    grid off;
case 'surfc'
    surfc(Xout,Yout,Zout);
    axis off;
    grid off;
    hold off
case 'contour' 
    axis equal;
    h = polar([theta_min theta_max], [Rho_min Rho_max]);
    delete(h)
    hold on
    contour(Xout,Yout,Zout,20);
    hold off
    colorbar;
case 'meshl' 
    % Set up mesh plot with polar axes and labels 
    angl = theta_min:abs(theta_max-theta_min)/(r-1):theta_max;
    Rho = ones([1 c])*Rho_max*1.005;
    X = Rho'*cos(angl);
    Y = Rho'*sin(angl);
    Z = zeros(size(X));
    % set up output data - this is exactly the same as the input
    Rho = Rho_min:abs(Rho_max-Rho_min)/(c-1):Rho_max;
    Xout = Rho'*cos(angl);
    Yout = Rho'*sin(angl);
    Zout = Zin';
    % plot the data
    axis equal;
    mesh(Xout,Yout,Zout);
    hold on
    % plot the axis
    mesh(X,Y,Z,'edgecolor',[0 0 0]);
    hold on;
    % set up tic mark and labels
    ticangle = round(theta_min*18/pi)*pi/18:pi/18:round(theta_max*18/pi)*pi/18;
    ticlength = [Rho_max*1.005 Rho_max*1.03];
    Xtic = ticlength'*cos(ticangle);
    Ytic = ticlength'*sin(ticangle);
    Ztic = zeros(size(Xtic));
    Xlbl = Rho_max*1.1*cos(ticangle);
    Ylbl = Rho_max*1.1*sin(ticangle);
    Zlbl = zeros(size(Xlbl));
    line(Xtic,Ytic,Ztic,'Color',[0 0 0]);
    if abs(theta_min-theta_max)==2*pi 
        Ntext = round(length(ticangle)/2)-1;
    else
        Ntext = round(length(ticangle)/2);
    end
    for i = 1:Ntext
        text(Xlbl(2*i-1),Ylbl(2*i-1),Zlbl(2*i-1),...
            num2str(ticangle(2*i-1)*180/pi),...
            'horizontalalignment','center')
    end
    set(gca,'DataAspectRatio',[1 1 1])  
    axis off;
    grid off; 
end
return


Contact us