% 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
