function varargout = FreqSystemIdent(varargin)
% FREQSYSTEMIDENT M-file for FreqSystemIdent.fig
%      FREQSYSTEMIDENT, by itself, creates a new FREQSYSTEMIDENT or raises the existing
%      singleton*.
%
%      H = FREQSYSTEMIDENT returns the handle to a new FREQSYSTEMIDENT or the handle to
%      the existing singleton*.
%
%      FREQSYSTEMIDENT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FREQSYSTEMIDENT.M with the given input arguments.
%
%      FREQSYSTEMIDENT('Property','Value',...) creates a new FREQSYSTEMIDENT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before FreqSystemIdent_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to FreqSystemIdent_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
% System Identifcation GUI - written by Dr Richard Mitchell, May 2011

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

% Last Modified by GUIDE v2.5 09-Apr-2014 14:52:03

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @FreqSystemIdent_OpeningFcn, ...
                   'gui_OutputFcn',  @FreqSystemIdent_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 FreqSystemIdent is made visible.
function FreqSystemIdent_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 FreqSystemIdent (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes FreqSystemIdent wait for user response (see UIRESUME)
% uiwait(handles.figure1);
setupsystem(handles);      % set up the system based on number in edit1

% --- Outputs from this function are returned to the command line.
function varargout = FreqSystemIdent_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 pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
setupsystem(handles);      % set up the system based on number in edit1


function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (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 edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
setupsystem(handles);       % set up the system based on number in edit1

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (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 button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Estimate gain
[m,p] = getmpdata(handles);                     % get m data from plot
if get(handles.checkbox1, 'value') > 0.5, mg = m(1); else mg = 1; end
hint = ['Low f gain is ', num2str(m(1))];       % define hint string
donewelement(handles, 1, mg, {'Gain'}, hint);   % include est of gain

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% take off lead
findleadorlag(handles, 0, 1, 45, {'TE'});   % look for lead where phase 45

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% take off lag
findleadorlag (handles, 0, 2, -45, {'TA'});  % look for lag where phase -45

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% take off integrator
donewelement(handles, 2, [], {'Integrator'}, '');       % include pure integrator

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% look for two pole
findquadpolezero(handles, 2, -90);      % look for quad pole, phase -90

% --- Executes on button press in pushbutton14.
function pushbutton14_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton14 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% find quadratic zero
findquadpolezero(handles, 1, 90);       % look for quad zero, phase +90

% --- Executes on button press in pushbutton7.
function pushbutton7_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% find lag using m
findleadorlag (handles, 1, 2, 1/sqrt(2), {'TA'});   
                            % look for lag where gain is dc/sqrt(2)

% --- Executes on button press in pushbutton8.
function pushbutton8_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% find lead using m
findleadorlag (handles, 1, 1, sqrt(2), {'TE'});
                            % look for lead where gain is dc * sqrt(2)

% --- Executes on button press in pushbutton9.
function pushbutton9_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton9 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% reveal system, ie show original TF
[numsys, densys] = getsystem(handles.text4, ' is ');
    % display string representation of system stored in text4.userdata
    
% --- Executes on button press in pushbutton10.
function pushbutton10_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton10 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% find a lead/lag or lag/lead
findleadlag (handles);      % fid Te and Ta by looking max phase lead/lag

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (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 edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double
tryestimate(handles);       % try estimate of element with para in edit2
                            % after user has adjusted parameter estimate
                            
% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (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 edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (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 edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.
tryestimate(handles);       % try estimate of element with para in edit3
                            % after user has adjusted parameter estimate

function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (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 button press in pushbutton12.
function pushbutton12_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton12 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% accept formally the latest estimate
set(handles.uipanel1, 'visible', 'off');    % hide parameter panel
set(handles.uipanel2, 'visible', 'on');     % show command panel


% --- Executes on button press in pushbutton13.
function pushbutton13_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton13 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Remove current guess
cancelestimate(handles);                    % take last estimate off
set(handles.uipanel1, 'visible', 'off');    % hide parametr panel
set(handles.uipanel2, 'visible', 'on');     % show command panel


function edit4_Callback(hObject, eventdata, handles)
% hObject    handle to edit4 (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 edit4 as text
%        str2double(get(hObject,'String')) returns contents of edit4 as a double

% --- Executes during object creation, after setting all properties.
function edit4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit4 (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 edit5_Callback(hObject, eventdata, handles)
% hObject    handle to edit5 (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 edit5 as text
%        str2double(get(hObject,'String')) returns contents of edit5 as a double


% --- Executes during object creation, after setting all properties.
function edit5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit5 (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 edit6_Callback(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit6 as text
%        str2double(get(hObject,'String')) returns contents of edit6 as a double


% --- Executes during object creation, after setting all properties.
function edit6_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit6 (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 edit7_Callback(hObject, eventdata, handles)
% hObject    handle to edit7 (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 edit7 as text
%        str2double(get(hObject,'String')) returns contents of edit7 as a double
setupsystem(handles);       % set up the system based on number in edit1

% --- Executes during object creation, after setting all properties.
function edit7_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit7 (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 button press in togglebutton1.
function togglebutton1_Callback(hObject, eventdata, handles)
% hObject    handle to togglebutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
print -dmeta 
% Copy to Clipboard


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

% Hint: get(hObject,'Value') returns toggle state of checkbox1


% ***********************************************

function findleadorlag (handles, mnotp, urow, uval, pstr)
% finds TC of lead or lag, by searching for uval in m or p graph
[m, p, w] = getmpdata(handles);     % get m,p,w from gain/phase graphs
if mnotp == 1, 
    TC = 1/interp1(m, w, uval);        % find time constant searching m
    hint = ['Gain is ', num2str(uval), ' at ', num2str(1/TC)];
else
    TC = 1/interp1(p, w, uval);        % find time constant searching p
    hint = ['Phase is ', num2str(uval), ' at ', num2str(1/TC)];
end
if get(handles.checkbox1, 'value') < 0.5, TC = 1; end
    % if not autocalculate parameters set TC to default value
donewelement (handles, 2+urow, TC, pstr, hint); % install element

function findleadlag (handles)
% find leadlag element
[m, p, w] = getmpdata(handles);     % get m,p,w from gain/phase graphs

if max(p) > p(1) + 5                % is Lead-lag
     aang = max(p) * pi/180;                % find angle at max phase
     llw = w(p==max(p));                    % ang freq where occurs
else                                % is Lag/Lead
     aang = min(p) * pi/180;                % angle at max phase lag
     llw = w(p==min(p));                    % and ang freq there
end;
afac = sqrt((1 - sin(aang)) / (1 + sin(aang))); % width of lead/lag
TE = 1/(llw * afac);                % find lead tc
TA = afac / llw;                    % and lag tc
if get(handles.checkbox1, 'value') < 0.5, 
    TE = 1;
    TA = 1;
end
hint = ['Phase is ', num2str(round(aang*180/pi)), ' at ', num2str(llw)];
donewelement(handles, 5, [TE, TA], {'TE', 'TA'}, hint); % install element

function findquadpolezero(handles, urow, pval)
% find quadratic pole/zero   
% urow = 1 if zero, so pval = +90 
% urow = 2 if poles, so pval = -90
% : find zeta, wn, but report as two poles/zero if zeta>=1
[m, p, w] = getmpdata(handles);

mg = 2*interp1(p, m, pval);     % gain where phase is +/- 90
wn = interp1(p, w, pval);       % wn is where phase is +/-90
zeta = m(1)/ mg;                % zeta is dc gain / 2 *gain at wn
hint = ['Phase is ', num2str(pval), ' at ', num2str(wn), ...
        ' when gain is ', num2str(mg/2), ' Low f gain is ', num2str(m(1))];

if zeta>=1,                     % if damping ratio > 1, model as two TCs
    pstr = {'T1', 'T2'}; 
    z2 = sqrt(zeta^2 - 1);
    pvals = (zeta + [z2,  - z2]) / wn;
else
    pstr = {'zeta', 'wn'}; 
    pvals = [zeta, wn];
end
if get(handles.checkbox1, 'value') < 0.5, 
    zeta = 0.5;
    wn = 1;
end
donewelement(handles, 5+urow, pvals, pstr, hint);

function [m, p, w] = getmpdata(handles)
% get [m,p,w] from remaining gain / phase graphs;
% limit them according to max w
ch = get(handles.axes4);                % from to be Id'ed phase graph
w = get(ch.Children, 'XData');          % get ang freqs
p = get(ch.Children, 'YData');          % get phase graphs
ch = get(handles.axes3);                % and from the gain graph
m = get(ch.Children, 'YData');          % get gains
maxw = str2num(get (handles.edit4, 'string'));
if maxw < max(w)                        % if specify only work up to maxw
    m = m(w<=maxw);                     % exclude values at higher w
    p = p(w<=maxw);
    w = w(w<=maxw);
end

function donewelement (handles, option, para, pstr, hintstr)
% user has selected an estimate, this is installed as follows
para = NsigFig(para, 2);
set(handles.text17, 'string', hintstr);
u = get(handles.text5, 'userdata');
set(handles.pushbutton13, 'userdata', u);    % remember what u was
set(handles.pushbutton12, 'userdata', option);    % remember what option was
if length(para) == 2,                               % if two parameters
    set(handles.edit3, 'string', num2str(para(2))); % set and show them
    set(handles.edit3, 'visible', 'on');
    set(handles.text8, 'string', pstr(2));
    set(handles.text8, 'visible', 'on');
else                                                % hide edit3/text8
    set(handles.edit3, 'visible', 'off');
    set(handles.text8, 'visible', 'off');
    set(handles.edit3, 'string', '0');
end
set(handles.text7, 'string', pstr(1));
if ~isempty(para),                                  % if one para
    set(handles.edit2, 'string', num2str(para(1))); % show it
    set(handles.edit2, 'visible', 'on');
else
    set(handles.edit2, 'string', '0');              % hide text8, edit2/3
    set(handles.edit2, 'visible', 'off');
    set(handles.edit3, 'visible', 'off');
    set(handles.text8, 'visible', 'off');
end
set(handles.uipanel1, 'visible', 'on');
set(handles.uipanel2, 'visible', 'off');
tryestimate (handles);                      % show effect of including estimate

function cancelestimate(handles)            
% undo estimate when user hits Cancel
u = get(handles.pushbutton13, 'userdata');
set(handles.text5, 'userdata', u);
goplot(handles);

function tryestimate (handles)
% show on graphs the effect of including estimate of next element
para = [str2num(get(handles.edit2, 'string')), ...
        str2num(get(handles.edit3, 'string'))];
u = get(handles.pushbutton13, 'userdata');
option = get(handles.pushbutton12, 'userdata');
switch option                                   % what type is it?
    case 1, u(1,[0:1]+findfreeu(u,1)) = [1, para(1)];   % just gain
    case 2, u(2,findfreeu(u,2)) = 2;                    % pure integrator
    case 3, u(1,[0:1]+findfreeu(u,1)) = [3, para(1)];   % zero
    case 4, u(2,[0:1]+findfreeu(u,2)) = [3, para(1)];   % pole
    case 5, u(1,[0:1]+findfreeu(u,1)) = [3, para(1)];   % lead lag
            u(2,[0:1]+findfreeu(u,2)) = [3, para(2)];
    case {6, 7}                                         % quadratic zero/pole
            urow = option - 5;
            if strcmp(get(handles.text7, 'string'), 'T1')
                u(urow,[0:1]+findfreeu(u,urow)) = [3, para(1)];
                u(urow,[0:1]+findfreeu(u,urow)) = [3, para(2)];
            else
                u(urow,[0:2]+findfreeu(u,urow)) = [4, para(1), para(2) ];
            end    
end
set(handles.text5, 'userdata', u);
goplot(handles);
    
function col = findfreeu (u, row)
% find next free column in give row in u
pos = 1;
while u(row,pos) > 0, pos = pos + 1; end    % keep looking for first 0
col = pos;

function setupsystem(handles)
set(handles.uipanel1, 'visible', 'off');    % hide parameter panel
set(handles.uipanel2, 'visible', 'on');     % show command panel
sysnum = get(handles.edit1, 'String');
if sysnum > '0'                             % if >0 hide user panel
    set(handles.uipanel3, 'visible', 'off');
else
    set(handles.uipanel3, 'visible', 'on');
end;    
if sysnum == '0', 
    u = getuser(handles.edit5, handles.edit6);  % get TF from user
else
    u = getusersys(handles, sysnum);            % get from system number
end    
set(handles.text4, 'userdata', u);              % string for what system is

set(handles.text5, 'UserData', zeros(2, 20));   % clear array used in Id
set(handles.text5, 'String', 'System unknown'); % estimate of system now

goplot(handles);                                % plot data
ch = get(handles.axes4);
w = get(ch.Children, 'XData');                  % find max w of graphs
set(handles.edit4, 'String', num2str(max(w)));  % string for maximum w

function goplot(handles)
% plot the two sets of bode plots
[numsys, densys] = getsystem(handles.text4, '');            % get n/d poly
[numest, denest] = getsystem(handles.text5, ' est as ');    % for est
[m, p, w] = bode(numsys, densys);                           % get freq resp
[me, pe] = bode(numest, denest, w);

myplot(handles.axes1, handles.axes2, m, p, w, me, pe);      
    % plot actual response and that of model so far
myplot(handles.axes3, handles.axes4, m./me, p-pe, w);
    % plot actual minus that so far ... ie what is left to show

function [numarr, denarr] = getsystem(texthandle, xtrstr)
% function to get polynom,ials and string representations of system
u = get(texthandle, 'userdata');
[numstr, numarr] = depoly(u(1, :));
[denstr, denarr] = depoly(u(2,:));
if ~isempty(xtrstr) 
    set(texthandle, 'String', ['System ', xtrstr, numstr, ' / ', denstr]);
else
    set(texthandle, 'String', '');
end

function myplot(axm, axp, m, p, w, me, pe)
% function to plot m and me vs w on axm, and p and pe on axp
% if me and pe not given, just plot m and p
% add appropriate labelling
axes(axm); cla;
if nargin>5, ch = loglog(w,m, w, me); else ch = loglog(w, m); end
set(ch, 'Linewidth', 2)
set(gca, 'xlim', [min(w) max(w)], 'xtick', [min(w) sqrt(min(w)*max(w)) max(w)]);
axes(axp); cla;
if nargin>5, ch = semilogx(w,p, w,pe); else ch = semilogx(w,p); end
ptop = max([0 45*ceil(max(p)/45)]);
pbot = 45*floor(min(p)/45);
set(ch, 'Linewidth', 2)
set(gca, 'ylim', [pbot-10, ptop+10], 'ytick', [pbot:45:ptop]);
set(gca, 'xlim', [min(w) max(w)], 'xtick', [min(w) sqrt(min(w)*max(w)) max(w)]);

function [pstr, parr] = depoly(pinfo)
% pinfo is array defining a polynomial
% function turns this into its string and polynomial representation
pstr = '';
parr = 1;
ct = 1;
while ct<=length(pinfo)
    if pinfo(ct) == 1             % gain
        pstr = [pstr, '*', num2str(pinfo(ct+1))];
        parr = parr * pinfo(ct+1);
        ct = ct + 2;
    elseif pinfo(ct) == 2           % pure s term
        pstr = [pstr, '*s'];
        parr = conv(parr, [1 0]);
        ct = ct + 1;
    elseif pinfo(ct) == 3           % time const
        ptc = pinfo(ct+1);
        pstr = [pstr, '*(s*', num2str(ptc), ' + 1)'];
        parr = conv(parr, [ptc 1]);
        ct = ct + 2;
    elseif pinfo(ct) == 4           % s/wn&^2 + 2zs/wn + 1
        zeta = pinfo(ct+1);
        wn = pinfo(ct+2);
        pstr = [pstr, '*(s^2*', num2str(1/wn^2), ' + s* ', num2str(2*zeta/wn), ' + 1)'];
        parr = conv(parr, [1/wn^2 2*zeta/wn 1]);
        ct = ct + 3;
    else
        ct = length(pinfo)+1;    % go to end
    end    
end
if ~isempty(pstr) && pstr(1)=='*', pstr = pstr(2:end); end

function u = getusersys (handles, num)
% Call FreqSystems function to get system defined by num
resp = FreqSystems (10+round(str2num(get(handles.edit1, 'String'))), get(handles.edit7, 'String'));
u(1, :) = doupoly(resp.pnum);
u(2, :) = doupoly(resp.pden);


function u = getuser (numedit, denedit)
% get user defined system ... using strings in numedit and denedit
u(1, :) = getuserpoly(numedit);
u(2, :) = getuserpoly(denedit);

function urow = getuserpoly (polyedit)
% convert string of numbers into polynomial array
urow = doupoly(str2num(get(polyedit, 'string')));

function urow = doupoly(actpoly)
% convert polynomial into data defining each corner freq
urow = zeros(1, 20);
urow(1) = 1;                                        % first value
urow(2) = actpoly(find(actpoly>0, 1, 'last' ));     % last one
pct = 3;
polyroots = -roots(actpoly);                        % find roots of rest
rct = 1;
while rct<=length(polyroots)                        % incorporate each root
    if polyroots(rct) == 0                          % r = 0
        urow(pct) = 2;
        pct = pct + 1;
    elseif isreal(polyroots(rct))                   % real root
        urow(pct) = 3;
        urow(pct+1) = 1/polyroots(rct);
        pct = pct + 2;
    else                                            % complex root
        urow(pct) = 4;
        urow(pct+1) = real(polyroots(rct));
        urow(pct+2) = 1/abs(polyroots(rct));
        pct = pct + 3;
        rct = rct + 1;
    end
    rct = rct + 1;    
end


function ans = NsigFig(X, SF)
% return number X to SF significant figures
N= num2str(X,SF);
ans= str2num(N);
