function forcedvibsim
%%
%% Simulator for 2nd order forced vibrations in MATLAB.
%% Written by Abhijit Dasgupta for MTH-3720 (UDM).
%% 
%% This simulation tool in MATLAB displays a second order
%% forced vibration system's response to sinusoidal input,
%% and has the following features:
%% 
%% - Fully simulate the sinusoidal response of any spring-mass-damper or
%%   any series RLC circuit (when time and frequency units are normalized)
%%
%% - View input and response together with gain and phase-lag on Bode plots
%% 
%% - Real-time slider control of parameters (input frequency, damping, etc)
%% 
%% - Option to display response split into transient and steady-state parts
%% 
%% - Simulate undamped response (beats, pure resonance) when delta set to 0
%% 
%% - Simulate free response by setting input amplitude to 0
%%
%% - For convenience, the natural undamped frequency omega_0 is fixed at 1,
%%   so that the input frequency remains normalized
%%
global hf hs hax hvals ptrTypes ptrID cc1b numFrames hvals DRAGstr ...
	t paramvals idx FRQ AMP IPH DMP QVL GNV PHL TLIM YLIM GLIM PLIM ...
	ctrlprm rspnsdisptyp tlimopt tscalar BodeGain BodePhase isSmooth;

isSmooth = 'off';

ptrTypes    = {'hand', 'lrdrag', 'arrow'};  % possible cursor pointers
ptrID       = 3;                            % default pointer

%% Slider colors
cbg = [1, 1, 1];
cfg = [0.1, 0.5, 0.05];

% Define colormap for the figure
cm = [
.3, .3, .3;
.9, .9, .9;
.2, .6, .2;
.8,  1, .8;
.2, .2, .6;
.8, .8,  1];

cc1b = [3 4 4 3]; % patch color for PLAYING

numFrames = 601;

tscalar = [ 1, 2, 5, 10, 20, 50, 100, 200, 500 ];

%% Control parameters strings

CPstrs = {...
	'Inp Frqncy \omega',...
	'Inp Ampltd A',...
	'Inp Phase \theta',...
      	'delta = Damping/(2m)',...
 	'Init Pos y(0)',...
     	'Init Vel y''(0)',...
      	'Q Factor',...
     	'Rspns Gain',...
       	'Rspns Phs-Lag'...
	};

DRAGstr = {...
	'Drag to select Inp Frqncy \omega',...
	'Drag to select Inp Ampltd A',...
	'Drag to select Inp Phase \theta',...
      	'Drag to select delta',...
 	'Drag to select Init Pos y(0)',...
     	'Drag to select Init Vel y''(0)',...
      	'Drag to select Q Factor',...
     	'Drag to select Rspns Gain',...
       	'Drag to select Rspns Phase-Lag'...
	};


FRQ = 1;
AMP = 2;
IPH = 3;
DMP = 4;
POS = 5;
VEL = 6;
QVL = 7;
GNV = 8;
PHL = 9;

YLIM = 5;
GLIM = 8.5;
PLIM = 3.5;
TLIM = tscalar(6);

t = (linspace(0,TLIM,801))';

%% col vecs
paramranges = [...
		[0,2.5]',...
		[0,YLIM]',...
	       	[0,10*pi]',...
	       	[0,2]',...
	       	[-YLIM,YLIM]',...
	       	[-YLIM,YLIM]'...
];
%% Fill in paramvals column by column (first col for FRQ, etc)
for k=1:6
    paramvals(:,k) = (linspace(paramranges(1,k), paramranges(2,k), numFrames))';
end
%% Refill damping values:  Make damping values log/exp scaled, and
%% fill in corresponding Q values
if 1 == 1
  base = 0.02 ;
  d0 = paramranges(1, DMP);
  d1 = paramranges(2, DMP);
  paramvals(:,DMP) = (d0 + (d1-d0)*log(linspace(1,base,numFrames))/log(base))';
end
%% Fill in Q values
paramvals(1,QVL) = inf;
for k=2:numFrames
    paramvals(k,QVL) = 0.5/paramvals(k,DMP);
end

%% The idx vector provides a snapshot of all the 6 parameters
idx(1) = ceil(numFrames/10);	%% freq = 0.25
idx(2) = ceil(2*numFrames/5);	%% amp = 2.0
idx(3) = 1;			%% inp phase = 0
idx(4) = 199;			%% delta = 0.2
idx(5) = ceil(6*numFrames/10);	%% y(0) = 1
idx(6) = round(numFrames/2);	%% y'(0) = 0


    hf = figure(...
      'Renderer'                        , 'painters', ...
      'Numbertitle'                     , 'off', ...
      'DockControls'                    , 'on', ...
      'Toolbar'                         , 'none', ...
      'Menubar'                         , 'none', ...
      'Color'                           , [.8 .8 .8], ...
      'Colormap'                        , cm, ...
      'Visible'				, 'off', ...
      'Units'                           , 'pixels', ...
      'Position'			, [100,100,900,730], ...
      'Color'                           , [.8 .8 .8], ...
      'DoubleBuffer'                    , 'on', ...
      'WindowButtonMotionFcn'           , @motionFcn, ...
      'Tag'                             , 'MainFigure', ...
      'defaultUIControlBackgroundColor' , [.8 .8 .8], ...
      'defaultAxesUnits'                , 'pixels', ...
      'defaultPatchEdgeColor'           , 'none');
      %% 'Colormap'                        , cm, ...

    % Time line bar
    hax(1) = axes(...
      'Parent'          , hf        , ...
      'Position'        , [30,35,560,25],...
      'Drawmode','fast',...
      'Visible'         , 'off'         , ...
      'Tag'             , 'TimelineAxes', ...
      'XLim'            , [0 1]         , ...
      'YLim'            , [-1 1]        );
     title('');
     xlabel('');
     ylabel('');
     set(hax(1), 'Visible'         , 'off' );
    patch([0, 0, 1, 1], [-.5, .5, .5, -.5], cbg, ...
      'Tag'             , 'timeLine', ...
      'ButtonDownFcn'   , @initiateDragFcn, ...
      'Parent'          , hax(1));
    patch([0, 0, 0, 0], [-.5, .5, .5, -.5], cc1b, ...
      'HitTest'         , 'off', ...
      'Tag'             , 'TimelinePatch', ...
      'CDataMapping'    , 'direct', ...
      'FaceColor'       , 'interp', ...
      'Parent'          , hax(1));
    % Message text boxes
    htxt(1) = uicontrol(...
      'Parent'              , hf, ...
      'Position'            , [280,55,120,20],...
      'style'               , 'text', ...
      'Tag'                 , 'FrameText', ...
      'HorizontalAlignment' , 'left', ...
      'FontWeight'          , 'Bold');
    htxt(2) = uicontrol(...
      'Parent'              , hf, ...
      'Position'            , [200,10,300,25],...
      'style'               , 'text', ...
      'Tag'                 , 'StatusText', ...
      'HorizontalAlignment' , 'center');
    % left range for control parameter
    htxt(3) = uicontrol(...
      'Parent'              , hf, ...
      'Position'            , [20,15,40,20],...
      'style'               , 'text', ...
      'Tag'                 , 'LeftRan', ...
      'HorizontalAlignment' , 'left', ...
      'FontWeight'          , 'Bold');
    htxt(4) = uicontrol(...
      'Parent'              , hf, ...
      'Position'            , [560,15,40,20],...
      'style'               , 'text', ...
      'Tag'                 , 'RightRan', ...
      'HorizontalAlignment' , 'right', ...
      'FontWeight'          , 'Bold');

    %%%%
    for k=0:2
        for l=1:3
	    m = 3*k + l;
	    hvals(m) = uicontrol(...
      		'Parent'              , hf, ...
      		'Position'            , [4+k*220,195-25*l,55,25],...
      		'style'               , 'text', ...
      		'HorizontalAlignment' , 'right', ...
      		'FontWeight'          , 'Bold');
    	    hlabs(m) = uicontrol(...
      		'Parent'              , hf, ...
      		'Position'            , [64+k*220,195-25*l,155-10*(k-1+abs(k-1)),25],...
      		'style'               , 'text', ...
      		'HorizontalAlignment' , 'left', ...
		'String'              , CPstrs(m),...
      		'FontWeight'          , 'Bold');
	end
    end

% 
   %  Construct the components.
   %% htext = uicontrol('Style','text','String','Select Q (damping)',...
   	  %% 'Position',[480,180,120,25]);

   hpop(1) = uicontrol('Style','popupmenu',...
	'FontWeight', 'Bold',...
	'String',...
	{'Input Freq omega',...
	'Input Amplitude A',...
	'Input Phase theta', ...
	'delta = Damping/(2m)',...
	'y(0) : Initial Position',...
	'y''(0) : Initial Velocity'},...
	'Position',[270,80,250,25],...
	'Tag', 'ControlParam',...
	'Callback',{@set_CtrlParam_Callback});
   htxt(5) = uicontrol('Style','text',...
	'FontWeight', 'Bold',...
   	'String','Control Parameter: ',...
	'Position',[80,80,190,20]);

   hpop(2) = uicontrol('Style','popupmenu',...
	'FontWeight', 'Bold',...
	'String',{'Show Full Response y(t) = y_c(t) + Y(t)',...
	'Show Only Steady State part Y(t)',...
	'Split Transient y_c(t) and Steady State Y(t)'}, ...
	'Position',[275,275,310,25],...
	'Tag', 'DisplayOption',...
	'Callback',{@set_Display_Option_Callback});
	%% 'Position',[100,200,60,20],...
   htxt(6) = uicontrol('Style','text',...
	'FontWeight', 'Bold',...
   	'String','Display Style for Response: ',...
	'Position',[35,275,240,20]);

   hpop(3) = uicontrol('Style','popupmenu',...
	'FontWeight', 'Bold',...
	'String',{'Time End: 1 ',...
	'Time End: 2 ',...
	'Time End: 5 ',...
	'Time End: 10 ',...
	'Time End: 20 ',...
	'Time End: 50 ',...
	'Time End: 100 ',...
	'Time End: 200 ',...
	'Time End: 500 '},...
	'Position',[480 320 125 15],...
	'Tag', 'TlimOption',...
	'Callback',{@set_Tlim_Option_Callback});

   % 
   %% Input and Response
   hax(2) = axes(...
   	'Units','Pixels',...
	'Position',[60,360,540,340],...
        'Tag', 'MainPlot',...
       	'XLim',[0, TLIM],...
	'YLim',[-YLIM, YLIM],...
   	'Drawmode','fast'); 
   title('Input and Response');
   xlabel('Normalize Time \omega_0 t');
   ylabel('');
   %% Bode Gain
   hax(3) = axes(...
   	'Units','Pixels',...
	'Position',[680,360,200,340],...
   	'Tag', 'GainPlot'); 
   %% Bode Phase Lage
   hax(4) = axes('Units','Pixels','Position',[680,60,200,160],...
   'Tag', 'PhasePlot'); 
   %% Dynamic label for values
   %ha4 = axes('Units','Pixels','Position',[60,40,540,75]); 
   %% ODE equation IVP
   hax(5) = axes('Units','Pixels','Position',[240,215,320,35]); 
   htxt(7) = uicontrol('Style','text',...
	'FontWeight', 'Bold',...
   	'String','Forced Vibrations ODE: ',...
	'Position',[70,215,160,20]);
   %% Input freq display
   %ha6 = axes('Units','Pixels','Position',[680,285,200,25]); 
   %% Input amplitude display
   %ha7 = axes('Units','Pixels','Position',[60,285,180,25]); 
   %% Initial conditions
   %ha8 = axes('Units','Pixels','Position',[265,185,200,25]);
   %% Damping and Q
   %ha9 = axes('Units','Pixels','Position',[480,186,120,24]);
%   
%
   for k=1:length(hax)
	   ManAxMode(hax(k));
	   set(hax(k), 'Visible'         , 'off' );
   end
   %
   hsld(1) =  uicontrol('Style','slider',...
        'Max',1,'Min',0,'Value',0.12,...
        'SliderStep',[0.01 0.1],...
        'Position',[15 585 15 120],...
        'Callback',@ylim_slider_callback);
   hsld(2) =  uicontrol('Style','slider',...
        'Max',1,'Min',0,'Value',0.25,...
        'SliderStep',[0.01 0.1],...
        'Position',[635 585 15 120],...
        'Callback',@glim_slider_callback);

   % Initialize the GUI.
   % Change units to normalized so components resize automatically.
   % set([f,ha1,ha2,ha3,ha4,ha5,ha6,ha7,ha8,ha9,hsa,hsb,hsc,hsd,...
   set([hf,hax,htxt,hvals,hlabs,hpop,hsld], 'Units','normalized');
   %
   % Assign the GUI a name to appear in the window title.
   set(hf,'Name','Sinusoidal Forced Vibrations')

   %% Get all handles for tag reference
   hs = guihandles(hf);
%


   % Initialize popoup menu hpopb and associated global var rspnsdisptyp
   ctrlprm = FRQ;
   set(hs.ControlParam, 'Value', ctrlprm);
   rspnsdisptyp = 1;
   set(hs.DisplayOption, 'Value', rspnsdisptyp);
   tlimopt = 6;
   set(hs.TlimOption, 'Value', tlimopt);
   TLIM = tscalar(tlimopt);
   t = (linspace(0,TLIM,801))';
   % Show ode
   dispode(hax(5), '');
   % Show Damping and Q
   % dispdq(ha9, 'Damping and Q');
   stackbode(1);
   stackfrms;
   drawset;

   sldval = log(YLIM*9.01/100+0.99)/log(10);
   set(hsld(1), 'Value', sldval);
   sldval = log(GLIM*9.01/100+0.99)/log(10);
   set(hsld(2), 'Value', sldval);

   update;

   for k=1:6
       set(hvals(k), 'String', sprintf('%5.3f', paramvals(idx(k), k)));
   end
   set(hvals(QVL), 'String', sprintf('%5.3f', paramvals(idx(DMP), QVL)));
   set(hvals(GNV), 'String', sprintf('%5.3f', BodeGain(idx(FRQ), idx(DMP))));
   set(hvals(PHL), 'String', sprintf('%5.3f', BodePhase(idx(FRQ), idx(DMP))));

   update;

   % Move the GUI to the center of the screen.
   movegui(hf,'center')
   % Make the GUI visible.
    set(hf, 'Visible' , 'on');


%  Callbacks for simple_gui. These callbacks automatically
%  have access to component handles and initialized data 
%  because they are nested at a lower level.


%%   %  Pop-up menu callbacks. Read the pop-up menu Value property
%%   %  to determine which item is currently displayed and make it
%%   %  the current Q value.
%%

function ylim_slider_callback(hObject,eventdata)
	global hax YLIM;
	sldval = get(hObject,'Value');
	YLIM = 100*(10^sldval-0.99)/9.01;
	%% sldval = log(YLIM*9.01/100+0.99)/log(10);
	%% set(hs.MainPlot, 'YLim', [-YLIM,YLIM]);
	set(hax(2), 'YLim', [-YLIM,YLIM]);

function glim_slider_callback(hObject,eventdata)
	global hax GLIM;
	sldval = get(hObject,'Value');
	GLIM = 100*(10^sldval-0.99)/9.01;
	%% set(hs.GainPlot, 'YLim', [0,GLIM]);
	set(hax(3), 'YLim', [0,GLIM]);


function set_CtrlParam_Callback(source,eventdata) 
	global ctrlprm;
	%% str = get(source, 'String');
	ctrlprm = get(source,'Value');
	stackfrms;
	drawset;
	update;

function set_Display_Option_Callback(source,eventdata) 
	global rspnsdisptyp;
	rspnsdisptyp = get(source,'Value');
	drawset;
	update;

function set_Tlim_Option_Callback(source,eventdata) 
	global tlimopt;
	%% str = get(source, 'String');
	tlimopt = get(source,'Value');
	stackfrms;
	drawset;
	update;


function initiateDragFcn(varargin)
    global hs;
    set(hs.MainFigure, 'WindowButtonMotionFcn', @timelineFcn);
    timelineFcn;
    set(hs.MainFigure, 'WindowButtonUpFcn', @stopDragFcn);


function timelineFcn(varargin)
    global hs numFrames idx ctrlprm;
    pt = get(hs.TimelineAxes, 'CurrentPoint');
    xVal = max([0, min([1, pt(1)])]);
    idx(ctrlprm) = max([1, round(numFrames * xVal)]);
    update;


% Called when the click-n-drag has completed
function stopDragFcn(varargin)
    global hs;
    set(hs.MainFigure, 'WindowButtonMotionFcn', @motionFcn);
    set(hs.MainFigure, 'WindowButtonUpFcn'    , ''        );


function motionFcn(varargin)
    global hs ptrTypes ptrID DRAGstr ctrlprm;
    %
    loc = get(hs.TimelineAxes, 'CurrentPoint');
    if isInRect(loc(1,1:2), [0, 1, -0.8, 0.8]) % mouse over time line bar
      if ptrID ~=2
        ptrID = 2;
        setptr(hs.MainFigure, ptrTypes{ptrID});
        set(hs.StatusText, 'String', DRAGstr(ctrlprm));
      end
    else                                          % everywhere else
      if ptrID ~= 3
        ptrID = 3;
        setptr(hs.MainFigure, ptrTypes{ptrID});
      end
      set(hs.StatusText, 'String', '');
    end


% Checks to see if pointer is in a rectangular region
function out = isInRect(xy, rect)
    out = (xy(1) > rect(1)) & (xy(1) < rect(2)) & ...
      (xy(2) > rect(3)) & (xy(2) < rect(4));


% Updates lines based on the signal type
function update
    global hs hvals numFrames idx ctrlprm BodeGain BodePhase YtStack ...
    	paramvals rspnsdisptyp FRQ DMP QVL GNV PHL;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    frm = idx(ctrlprm);
    val = paramvals(frm, ctrlprm);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Update time line slider bar
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    set(hs.TimelinePatch , 'XData', [0, 0, frm/numFrames, frm/numFrames]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Update frame number text
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    set(hs.FrameText, 'String', sprintf('Value: %5.3f', val));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Update the main axes plot
    %% t = input ; yc = transient ; Yp = steady state;
    %% Colors: Input = grn, Yp = red, yc = blue, y = magenta
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if rspnsdisptyp == 1	% plot combined y = yc + Yp
	%% set(hs.Lines, {'YData'} , YCELL_14(:, :, frm)'                 );
	set(hs.Lines(2), 'YData', YtStack(:,1,frm));
	set(hs.Lines(1), 'YData', YtStack(:,4,frm));
    elseif rspnsdisptyp == 2	% plot only Yp
	%% set(hs.Lines, {'YData'} , YCELL_13(:, :, frm)'                 );
	set(hs.Lines(2), 'YData', YtStack(:,1,frm));
	set(hs.Lines(1), 'YData', YtStack(:,3,frm));
    else			% split yc and Yp
	%% set(hs.Lines, {'YData'} , YCELL_123(:, :, frm)'                 );
	set(hs.Lines(3), 'YData', YtStack(:,1,frm));
	set(hs.Lines(1), 'YData', YtStack(:,2,frm));
	set(hs.Lines(2), 'YData', YtStack(:,3,frm));
    end
    %%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Update Control Param numerical value text fields
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Updating parameter values numeric text fields:
	% hval is a vector of length 6 + 3 = 9, with hval(1) being
	% a handle to the text element (axes or direct) for showing
	% the numerical value of the input freq, etc, etc,
	% and hval(7) = hval(QFAC) is a handle to a text element
	% for showing the numerical value of Q, hval(8) = hval(GAIN)
	% for showing gain, and hval(9) = hval(PHASE) for phase-lag
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    set(hvals(ctrlprm), 'String', sprintf('%5.3f', val));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For special parameters DMP and FRQ, update gain curves and markers,
    % and also the additional numerical value text fields
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (ctrlprm == DMP)
	set(hvals(QVL), 'String', sprintf('%5.3f', paramvals(frm, QVL)));
	set(hs.GLin, 'YData', BodeGain(:,frm));
	set(hs.PLin, 'YData', BodePhase(:,frm));
    end
    if (ctrlprm == DMP | ctrlprm == FRQ)
	frq = paramvals(idx(FRQ),FRQ);
	gain = BodeGain(idx(FRQ), idx(DMP));
	phsl = BodePhase(idx(FRQ), idx(DMP));
	set(hvals(GNV), 'String', sprintf('%5.3f', gain));
	set(hvals(PHL), 'String', sprintf('%5.3f', phsl));
	set(hs.GMrk, 'XData', frq, 'YData', gain);
	set(hs.PMrk, 'XData', frq, 'YData', phsl);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% The ODE IVP
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dispode(ax, ttl)
axes(ax);
plot([0], [0]);
set(ax, 'XTick', []);
set(ax, 'YTick', []);
set(ax, 'Xlim', [0,10], 'Ylim', [0, 1]);
text(0.5, 0.25,...
    'y'''' + 2{\delta}y'' + {\omega_0^2}y = {\omega_0^2}(A cos(\omega{t} - \theta));  with  \omega_0 = 1. ');
title(ttl);
% 
% function dispdq(ax, ttl)
% axes(ax);
% plot([0], [0]);
% set(ax, 'XTick', []);
% set(ax, 'YTick', []);
% set(ax, 'Xlim', [0,10], 'Ylim', [0, 1]);
% text(2, 0.3, 'Q := \omega_0/(2\delta)');
% title(ttl);
%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Solve IVP  y'' + 2*dmp*y' + (wn^2)*y = wn^2 * (A * cos(w*t - inphs))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yc, Yp] = solvivp(w, amp, inphs, dmp, ic1, ic2, wn)
	global t;	% series of time points column vector (so yc, Yp also col vecs)
	inp = amp*cos(w*t - inphs);
	if (dmp > 0) % DAMPED SYSTEM
		denom = (wn^2 - w^2)^2 + 4*(dmp^2)*(w^2);
		gain = wn^2/sqrt(denom);
		phslag = acos((wn^2 - w^2)/sqrt(denom));
		Yp = gain*amp*cos(w*t - inphs - phslag);
		%% Yz = (wn^2)*(wn^2 - w^2)/denom;
		%% Ypz = 2*dmp*(wn^2)*(w^2)/denom;
		Yz = gain*amp*cos(inphs + phslag);
		Ypz = w*gain*amp*sin(inphs + phslag);
		if (dmp > wn)
			if (wn == 0)
				wn = 0.5e-5;
			end
			dsr = sqrt(dmp^2 - wn^2);  % offset of each root from their midpt -dmp
			rt1 = -dmp + dsr;
			rt2 = -dmp - dsr;
			c1 = (ic2 - ic1*rt2 - Ypz + rt2*Yz)/(2*dsr);
			c2 = (ic1*rt1 - ic2 + Ypz - rt1*Yz)/(2*dsr);
			yc = c1*exp(rt1*t) + c2*exp(rt2*t);
		elseif (dmp == wn)
			c1 = ic1 - Yz;
			c2 = ic1*dmp + ic2 - dmp*Yz - Ypz;
			yc = (c1 + c2*t).*exp(-dmp*t);
		else	% dmp < wn
			nu = sqrt(wn^2 - dmp^2);  % quasi-freq
			c1 = ic1 - Yz;
			c2 = (ic1*dmp + ic2 -  dmp*Yz - Ypz)/nu;
			yc = (c1*cos(nu*t) + c2*sin(nu*t)).*exp(-dmp*t);
		end
	elseif (dmp == 0)	% UNDAMPED SYSTEM
		if (w ~= wn)	% beats;
			gain = 1/abs(wn^2 - w^2);
			if (w < wn)
				phslag = 0;
			else
				phslag = pi;
			end
			Yp = (amp/(wn^2-w^2))*cos(w*t - inphs);
			Yz = (amp/(wn^2-w^2))*cos(inphs);
			Ypz = (amp*w/(wn^2-w^2))*sin(inphs);
		else	% pure resonance
			gain = inf;
			phslag = pi/2;
			Yp = (amp/(2*wn))*t.*sin(wn*t - inphs);
			Yz = 0;
			Ypz = (-amp/(2*wn))*sin(inphs);
		end
		c1 = ic1 - Yz;
		c2 = (ic2 - Ypz)/wn;
		yc = c1*cos(wn*t) + c2*sin(wn*t);
	else % (dmp < 0)
		gain = nan;
		phslag = nan;
		yc = nan*t;
		Yp = nan*t;
	end


%% Construct all Bode values (need to be done just once at start of program)
function stackbode(wn)
	global paramvals BodeGain BodePhase numFrames idx FRQ DMP GLIM;
	for k = 1:numFrames
		for l = 1:numFrames
			w = paramvals(l, FRQ);
			dmp = paramvals(k, DMP);
			denom = (wn^2 - w^2)^2 + 4*(dmp^2)*(w^2);
			if (dmp == 0)
				if (w < wn)
					gain = 1/abs(wn^2 - w^2);
					phslag = 0;
				elseif (w > wn)
					gain = 1/abs(wn^2 - w^2);
					phslag = pi;
				else % pure resonance
					gain = inf;
					phslag = pi/2;
				end
			else
				gain = wn^2/sqrt(denom);
				phslag = acos((wn^2 - w^2)/sqrt(denom));
			end
			if (gain < GLIM)
				BodeGain(l,k) = gain;
			else
				%% BodeGain(l,k) = inf;
				BodeGain(l,k) = gain;
			end
			BodePhase(l,k) = phslag;
		end
	end


%% Call this at init and whenever ctrlprm changes
function stackfrms
	global t tscalar tlimopt paramvals idx ctrlprm YtStack FRQ AMP IPH YLIM TLIM numFrames ;
	TLIM = tscalar(tlimopt);
	t = (linspace(0,TLIM,801))';
	%% Parameter ranges defined (may be global)
	%% prmran = [0, 2.5; 0, 2; 0, 2*pi; 0, 3; -3, 3; -3, 3];
	%% p = linspace(prmran(ctrlprm,1), prmran(ctrlprm,2), numFrames);
	%% u = [w, amp, inphs, dmp, ic1, ic2];  %% GLOBAL VALS
	u = [paramvals(idx(1),1), paramvals(idx(2),2), paramvals(idx(3),3),...
			paramvals(idx(4),4),paramvals(idx(5),5),paramvals(idx(6),6)];
	wn = 1;
	for k = 1:numFrames
		u(ctrlprm) = paramvals(k,ctrlprm);
		w = u(FRQ);
		amp = u(AMP);
		inphs = u(IPH);
		[yc, Yp] = solvivp(u(1), u(2), u(3), u(4), u(5), u(6), wn);
		YtStack(:,1,k) = amp*cos(w*t - inphs);
		YtStack(:,2,k) = yc ;
		YtStack(:,3,k) = Yp ;
		YtStack(:,4,k) = yc + Yp ;
		%% YCELL_14 = num2cell( YtStack(:,[1,4],:), 1);
		%% YCELL_13 = num2cell( YtStack(:,[1,3],:), 1);
		%% YCELL_123 = num2cell( YtStack(:,[1,2,3],:), 1);
		if 0 == 1
		for j = 1:length(t)
			for l = 1:4
				if (abs(YtStack(j,l,k)) >= YLIM)
					YtStack(j,l,k) = nan;
				end
			end
		end
		end
	end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Function to plot at initialization or reset (change of control param)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function drawset
	global hf hs hax rspnsdisptyp t YtStack BodeGain BodePhase ...
	FRQ DMP QVL GNV PHL paramvals idx ctrlprm ...
	tscalar tlimopt TLIM YLIM GLIM PLIM isSmooth numFrames ;
	TLIM = tscalar(tlimopt);
	t = (linspace(0,TLIM,801))';
	frm = idx(ctrlprm);
	% ax = hs.MainPlot;
	ax = hax(2);
	axes(ax);
	%% t = input ; yc = transient ; Yp = steady state;
	%% Colors: Input = grn, Yp = red, yc = blue, y = magenta
	if rspnsdisptyp == 1		% plot combined y = yc + Yp
    	hplt = plot(t, YtStack(:,1,frm), 'g-',...
		t, YtStack(:,4,frm), 'm-', 'LineWidth', 2, 'Tag', 'Lines');
    	%% legend(ax, 'Input A cos(\omega{t}-\theta)', 'Response y_c(t) + Y(t)');
	elseif rspnsdisptyp == 2	% plot only Yp
    	hplt = plot(t, YtStack(:,1,frm), 'g-',...
   		t, YtStack(:,3,frm), 'r-', 'LineWidth', 2, 'Tag', 'Lines');
    	%% legend(ax, 'Input A cos(\omega{t}-\theta)', 'Steady State Y(t)');
	else					% split yc and Yp
    	hplt = plot(t, YtStack(:,1,frm), 'g-',...
		t, YtStack(:,3,frm), 'r-',...
		t, YtStack(:,2,frm), 'b-', 'LineWidth', 2, 'Tag', 'Lines');
    	%% legend(ax, 'Input A cos(\omega{t}-\theta)', 'Steady State Y(t)', 'Transient y_c(t)');
	end
	if strcmp(isSmooth, 'on')
		try
			set(hL, 'LineSmoothing', isSmooth);
		catch
			disp('Anti-aliasing not supported.');
		end
	end
	set(ax, 'FontWeight', 'Bold');
	xlabel('Normalize Time \omega_0 t');
	ylabel('');
	set(ax, 'Xlim', [0, TLIM], 'Ylim', [-YLIM, YLIM]);
	grid(ax, 'on');
	title('Input: Green, Full-Rspns: Magenta, SS: Red, Tramsient: Blue');
	%%%
	%% BodeGain
	%%%
	% ax = hs.GainPlot;
	ax = hax(3);
	axes(ax);
	hold off;
	plot(paramvals(:,FRQ), BodeGain(:,idx(DMP)),...
		'k-', 'LineWidth', 2, 'Tag', 'GLin');
	hold on;
	plot(paramvals(idx(FRQ),FRQ), BodeGain(idx(FRQ),idx(DMP)),...
		'rd', 'LineWidth', 3, 'MarkerSize', 8, 'Tag', 'GMrk');
	set(ax, 'FontWeight', 'Bold');
	set(ax, 'Xlim', [paramvals(1,FRQ),paramvals(numFrames,FRQ)],...
		'Ylim', [0, GLIM]);
	xlabel('Normalized Input Freq \omega/\omega_0');
	ylabel('Gain');
	grid(ax, 'on');
	hold off;
	title('SS Gain (Bode Plot)');
	%%%
	%% BodePhase
	%%%
	% ax = hs.PhasePlot;
	ax = hax(4);
	axes(ax);
	hold off;
	plot(paramvals(:,FRQ), BodePhase(:,idx(DMP)),...
		'k-', 'LineWidth', 2, 'Tag', 'PLin');
	hold on;
	plot(paramvals(idx(FRQ),FRQ), BodePhase(idx(FRQ),idx(DMP)),...
		'md', 'LineWidth', 3, 'MarkerSize', 8, 'Tag', 'PMrk');
	set(ax, 'FontWeight', 'Bold');
	set(ax, 'Xlim', [paramvals(1,FRQ),paramvals(numFrames,FRQ)],...
		'Ylim', [0, PLIM]);
	xlabel('Normalized Input Freq \omega/\omega_0');
	ylabel('Phase Lag');
	grid(ax, 'on');
	hold off;
	title('SS Phase-Lag (Bode Plot)');
	set(hs.LeftRan, 'String',...
		sprintf('%4.2f', paramvals(1,ctrlprm)));
	set(hs.RightRan, 'String',...
		sprintf('%4.2f', paramvals(numFrames,ctrlprm)));
	hs = guihandles(hf);



function ManAxMode(h)
    % Do not set CameraViewAngleMode, DataAspectRatioMode,
    % and PlotBoxAspectRatioMode to aviod exposing a bug
    pn = {'ALimMode',...
        'CameraPositionMode','CameraTargetMode',...
        'CameraUpVectorMode','CLimMode',...
        'TickDirMode','XLimMode',...
        'YLimMode','ZLimMode',...
        'XTickMode','YTickMode',...
        'ZTickMode','XTickLabelMode',...
        'YTickLabelMode','ZTickLabelMode'};
    for k = 1:15
        pv(k) = {'manual'};
    end
    set(h,pn,pv);