%% Drag 
function drag(~,~)
close all

global Input1 Input2 s_1 s_2 s_3 s_4 s_5 fontsize linewidth markersize displayscale; %Define this to be global so subfunction can input and access figure

% displayscale scales pixels and fonts, displayscale ~ 1 for a projector,
% displayscale ~ 2 for an HD monitor
displayscale = 1.2;

% main figure that hold the GUI
figure('position', displayscale*[0 0 700 350], 'name', 'Drag', 'NumberTitle', 'off');

% better display values for plotting
fontsize=displayscale*10;
linewidth=displayscale;
markersize=displayscale*5;

% subplots to hold figures
s_1 = subplot('position',[0.08 0.3 0.33 0.65]);
xlabel('x [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
ylabel('y [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

s_2 = subplot('position',[0.48 0.57 0.2 0.38]);
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('x [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

s_3 = subplot('position',[0.75 0.57 0.2 0.38]);
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('y [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

s_4 = subplot('position',[0.48 0.1 0.2 0.38]);
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('v$_{\rm x}$ [m/s] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

s_5 = subplot('position',[0.75 0.1 0.2 0.38]);
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('v$_{\rm y}$ [m/s] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

% descriptive text for parameters
uibutton('Style', 'text', 'String', 'v$_{\rm xo}$/v$_{\rm t}$',...
'Position', displayscale*[40 45 40 20],'FontSize',fontsize,'Interpreter','latex');

uibutton('Style', 'text', 'String', 'v$_{\rm yo}$/v$_{\rm t}$',...
'Position', displayscale*[40 20 40 20],'FontSize',fontsize,'Interpreter','latex');

% text input for varying initial parameters.
Input1 = uicontrol('Style', 'edit', 'Position', displayscale*[90 45 40 20],'FontSize',fontsize);

Input2 = uicontrol('Style', 'edit', 'Position', displayscale*[90 20 40 20],'FontSize',fontsize);

% button to run the simulation
uicontrol('Style', 'pushbutton', 'String', 'Run', 'FontSize',fontsize,...
'Position', displayscale*[140 30 60 30],'Callback', @MakePlots);

% button to clear past simulation
uicontrol('Style', 'pushbutton', 'String', 'Clear', 'FontSize',fontsize,...
'Position', displayscale*[210 30 60 30],'Callback', @drag);

%% Called by drag to do the plotting
function MakePlots(~,~)
global Input1 Input2 s_1 s_2 s_3 s_4 s_5 fontsize linewidth markersize;

% gets the values of the initial parameters from the sliders
vxo = get(Input1,'String');
vyo = get(Input2,'String');

vxo = str2double(vxo);
vyo = str2double(vyo);

refine = true;
N_fine = 100;
N = 100;

g = 9.8;
vt = 1;     %all velocities normalized to a terminal velocity of 1

R = 2*vxo*vyo/g;    %vacuum range
to = 0;             %starting time for all solutions
tf = 2.5*R/vxo;     %ending time all solutions, n times the vacuum range time
tau = vt/g;         %e-folding time
xinf = vxo*tau;     %x asymptote for linear drag

time = 0:tf/N:tf;
x = vxo*tau*(1-exp(-time/tau)); %linear drag trajectory
y = (vyo+vt)*tau*(1-exp(-time/tau)) - vt*time;
xvac = vxo*time;    %vacuum (no drag) trajectory
yvac = vyo*time-0.5*g*time.^2;

options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6]);
%options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6]);   %Matlab default values
[t,v] = ode45(@accelerations,[to tf],[vxo vyo],options);
%the syntax for any ODE integrator is:
%[t, f(t)] = odesolver(@(function that returns [t, df/dt]), [(start time) (end time)], [f(t=0) df/dt(t=0)])

vxq = v(:,1);
vyq = v(:,2);

%refining velocities from ODE integration to improve position accuracy
if refine
    t_fine = linspace(t(1),t(end),N_fine*length(t))';
    vxq = spline(t,vxq,t_fine);
    vyq = spline(t,vyq,t_fine);
    t = t_fine;
end

dt = diff(t);
dt(end+1) = dt(end);
xfine = cumsum(dt.*vxq)';
yfine = cumsum(dt.*vyq)';
t = t';

% end up with positions and velocities for quadratic drag at same times as
% linear and vacuum samples
xq = interp1(t,xfine,time,'spline','extrap');
yq = interp1(t,yfine,time,'spline','extrap');
vxq = interp1(t,vxq,time,'spline','extrap');
vyq = interp1(t,vyq,time,'spline','extrap');

vx = diff(x)./diff(time);
vy = diff(y)./diff(time);
vxvac = diff(xvac)./diff(time);
vyvac = diff(yvac)./diff(time);
vx(end+1) = vx(end);
vy(end+1) = vy(end);
vxvac(end+1) = vxvac(end);
vyvac(end+1) = vyvac(end);

% plot of trajectory, y vs. x
subplot(s_1)
plot(xvac,yvac,'b--','LineWidth',linewidth)
hold all
h_1 = plot(x(1),y(1),'r','LineWidth',linewidth);
h_3 = plot(xq(1),yq(1),'m','LineWidth',linewidth);
%legend({'vacuum','linear','quadratic'},'FontSize',fontsize,'Interpreter','latex');
h_2 = plot(x(1),y(1),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
h_4 = plot(xq(1),yq(1),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
xlim([0 R])
ylim([0 max(yvac)])
%ylim([min(yvac) max(yvac)])
xlabel('x [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
ylabel('y [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

% the other plots show position and momentum components vs. time
subplot(s_2)
plot(time,xvac,'b--','LineWidth',linewidth)
hold all
plot(time,xinf.*ones(1,length(time)),'k:','LineWidth',linewidth)
h_5 = plot(time(1),x(1),'r','LineWidth',linewidth);
h_6 = plot(time(1),x(1),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
h_7 = plot(time(1),xq(1),'m','LineWidth',linewidth);
h_8 = plot(time(1),xq(1),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
xlim([time(1) time(end)])
ylim([0 max(xvac)])
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('x [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(s_3)
plot(time,yvac,'b--','LineWidth',linewidth)
hold all
h_9 = plot(time(1),y(1),'r','LineWidth',linewidth);
h_10 = plot(time(1),y(1),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
h_11 = plot(time(1),yq(1),'m','LineWidth',linewidth);
h_12 = plot(time(1),yq(1),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
xlim([time(1) time(end)])
ylim([min(yvac) max(yvac)])
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('y [m] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(s_4)
plot(time,vxvac,'b--','LineWidth',linewidth)
hold all
plot(time,vt*ones(1,length(time)),'k:','LineWidth',linewidth)
h_13 = plot(time(1),vx(1),'r','LineWidth',linewidth);
h_14 = plot(time(1),vx(1),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
h_15 = plot(time(1),vxq(1),'m','LineWidth',linewidth);
h_16 = plot(time(1),vxq(1),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
xlim([time(1) time(end)])
ylim([0 1.2*max(vxvac)])
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('v$_{\rm x}$ [m/s] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(s_5)
plot(time,vyvac,'b--','LineWidth',linewidth)
hold all
plot(time,-vt*ones(1,length(time)),'k:','LineWidth',linewidth)
h_17 = plot(time(1),vy(1),'r','LineWidth',linewidth);
h_18 = plot(time(1),vy(1),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
h_19 = plot(time(1),vyq(1),'m','LineWidth',linewidth);
h_20 = plot(time(1),vyq(1),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
xlim([time(1) time(end)])
ylim([min(vyvac) max(vyvac)])
xlabel('t [s]','FontSize',fontsize,'Interpreter','latex')
ylabel('v$_{\rm y}$ [m/s] ($\times$ v$_{\rm t}$ [m/s])','FontSize',fontsize,'Interpreter','latex')
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

% looping through time to make an on-screen animation
for i=2:N
    pause(0.1)
    subplot(s_1)
    delete(h_1)
    delete(h_2)
    delete(h_3)
    delete(h_4)
    h_1 = plot(x(1:i),y(1:i),'r','LineWidth',linewidth);
    h_2 = plot(x(i),y(i),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
    h_3 = plot(xq(1:i),yq(1:i),'m','LineWidth',linewidth);
    h_4 = plot(xq(i),yq(i),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
    
    subplot(s_2)
    delete(h_5)
    delete(h_6)
    delete(h_7)
    delete(h_8)
    h_5 = plot(time(1:i),x(1:i),'r','LineWidth',linewidth);
    h_6 = plot(time(i),x(i),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
    h_7 = plot(time(1:i),xq(1:i),'m','LineWidth',linewidth);
    h_8 = plot(time(i),xq(i),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');

    subplot(s_3)
    delete(h_9)
    delete(h_10)
    delete(h_11)
    delete(h_12)
    h_9 = plot(time(1:i),y(1:i),'r','LineWidth',linewidth);
    h_10 = plot(time(i),y(i),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
    h_11 = plot(time(1:i),yq(1:i),'m','LineWidth',linewidth);
    h_12 = plot(time(i),yq(i),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');

    subplot(s_4)
    delete(h_13)
    delete(h_14)
    delete(h_15)
    delete(h_16)
    h_13 = plot(time(1:i),vx(1:i),'r','LineWidth',linewidth);
    h_14 = plot(time(i),vx(i),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
    h_15 = plot(time(1:i),vxq(1:i),'m','LineWidth',linewidth);
    h_16 = plot(time(i),vxq(i),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');

    subplot(s_5)
    delete(h_17)
    delete(h_18)
    delete(h_19)
    delete(h_20)
    h_17 = plot(time(1:i),vy(1:i),'r','LineWidth',linewidth);
    h_18 = plot(time(i),vy(i),'ro','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','r','MarkerFaceColor','r');
    h_19 = plot(time(1:i),vyq(1:i),'m','LineWidth',linewidth);
    h_20 = plot(time(i),vyq(i),'mo','LineWidth',linewidth,'MarkerSize',markersize,'MarkerEdgeColor','m','MarkerFaceColor','m');
end

%% x and y accelerations under the influence of quadratic drag
function dvdt = accelerations(~,v)
vt = 1;
g = 9.8;
c = g/vt^2; % quadratic drag constant, m=1 is implicit

dvdt = zeros(2,1);    % a column vector for [(dv/dt)_x, (dv/dt)_y]
dvdt(1) = -c*v(1)*sqrt(v(1)^2+v(2)^2);             % this is (dv/dt)_x
dvdt(2) = -g-c*v(2)*sqrt(v(1)^2+v(2)^2);    % this is (dv/dt)_y, m=1 is implicit

%% The following was taken from http://www.mathworks.com/matlabcentral/fileexchange/10743-uibutton--gui-pushbuttons-with-better-labels 
function [hout,ax_out] = uibutton(varargin)
%uibutton: Create pushbutton with more flexible labeling than uicontrol.
% Usage:
%   uibutton accepts all the same arguments as uicontrol except for the
%   following property changes:
%
%     Property      Values
%     -----------   ------------------------------------------------------
%     Style         'pushbutton', 'togglebutton' or 'text', default =
%                   'pushbutton'.
%     String        Same as for text() including cell array of strings and
%                   TeX or LaTeX interpretation.
%     Interpreter   'tex', 'latex' or 'none', default = default for text()
%     Rotation      text rotation angle, default = 0
%
% Syntax:
%   handle = uibutton('PropertyName',PropertyValue,...)
%   handle = uibutton(parent,'PropertyName',PropertyValue,...)
%   [text_obj,axes_handle] = uibutton('Style','text',...
%       'PropertyName',PropertyValue,...)
%
% uibutton creates a temporary axes and text object containing the text to
% be displayed, captures the axes as an image, deletes the axes and then
% displays the image on the uicontrol.  The handle to the uicontrol is
% returned.  If you pass in a handle to an existing uicontol as the first
% argument then uibutton will use that uicontrol and not create a new one.
%
% If the Style is set to 'text' then the axes object is not deleted and the
% text object handle is returned (as well as the handle to the axes in a
% second output argument).
%
% See also UICONTROL.

% Version: 1.9, 5 November 2010
% Author:  Douglas M. Schwarz
% Email:   dmschwarz=ieee*org, dmschwarz=urgrad*rochester*edu
% Real_email = regexprep(Email,{'=','*'},{'@','.'})


% Detect if first argument is a uicontrol handle.
keep_handle = false;
if nargin > 0
	h = varargin{1};
	if isscalar(h) && ishandle(h) && strcmp(get(h,'Type'),'uicontrol')
		keep_handle = true;
		varargin(1) = [];
	end
end

% Parse arguments looking for 'Interpreter' property.  If found, note its
% value and then remove it from where it was found.
interp_value = get(0,'DefaultTextInterpreter');
rotation_value = get(0,'DefaultTextRotation');
arg = 1;
remove = [];
while arg <= length(varargin)
	v = varargin{arg};
	if isstruct(v)
		fn = fieldnames(v);
		for i = 1:length(fn)
			if strncmpi(fn{i},'interpreter',length(fn{i}))
				interp_value = v.(fn{i});
				v = rmfield(v,fn{i});
			elseif strncmpi(fn{i},'rotation',length(fn{i}))
				rotation_value = v.(fn{i});
				v = rmfield(v,fn{i});
			end
		end
		varargin{arg} = v;
		arg = arg + 1;
	elseif ischar(v)
		if strncmpi(v,'interpreter',length(v))
			interp_value = varargin{arg+1};
			remove = [remove,arg,arg+1]; %#ok<AGROW>
		elseif strncmpi(v,'rotation',length(v))
			rotation_value = varargin{arg+1};
			remove = [remove,arg,arg+1]; %#ok<AGROW>
		end
		arg = arg + 2;
	elseif arg == 1 && isscalar(v) && ishandle(v) && ...
			any(strcmp(get(h,'Type'),{'figure','uipanel'}))
		arg = arg + 1;
	else
		error('Invalid property or uicontrol parent.')
	end
end
varargin(remove) = [];

% Create uicontrol, get its properties then hide it.
if keep_handle
	set(h,varargin{:})
else
	h = uicontrol(varargin{:});
end
s = get(h);
if ~any(strcmp(s.Style,{'pushbutton','togglebutton','text'}))
	delete(h)
	error('''Style'' must be pushbutton, togglebutton or text.')
end
set(h,'Visible','off')

% Create axes.
parent = get(h,'Parent');
ax = axes('Parent',parent,...
	'Units',s.Units,...
	'Position',s.Position,...
	'XTick',[],'YTick',[],...
	'XColor',s.BackgroundColor,...
	'YColor',s.BackgroundColor,...
	'Box','on',...
	'Color',s.BackgroundColor);
% Adjust size of axes for best appearance.
set(ax,'Units','pixels')
pos = round(get(ax,'Position'));
if strcmp(s.Style,'text')
	set(ax,'Position',pos + [0 1 -1 -1])
else
	set(ax,'Position',pos + [4 4 -8 -8])
end
switch s.HorizontalAlignment
	case 'left'
		x = 0.0;
	case 'center'
		x = 0.5;
	case 'right'
		x = 1;
end
% Create text object.
text_obj = text('Parent',ax,...
	'Position',[x,0.5],...
	'String',s.String,...
	'Interpreter',interp_value,...
	'Rotation',rotation_value,...
	'HorizontalAlignment',s.HorizontalAlignment,...
	'VerticalAlignment','middle',...
	'FontName',s.FontName,...
	'FontSize',s.FontSize,...
	'FontAngle',s.FontAngle,...
	'FontWeight',s.FontWeight,...
	'Color',s.ForegroundColor);

% If we are creating something that looks like a text uicontrol then we're
% all done and we return the text object and axes handles rather than a
% uicontrol handle.
if strcmp(s.Style,'text')
	delete(h)
	if nargout
		hout = text_obj;
		ax_out = ax;
	end
	return
end

% Determine parent figure and move it to main screen, if necessary.
% (Necessary because of bug in getframe.)
fig = parent;
while ~strcmp(get(fig,'Type'),'figure')
	fig = get(fig,'Parent');
end
fig_pos = get(fig,'Position');
movegui(fig)

% Capture image of axes and then delete the axes.
frame = getframe(ax);
delete(ax)

% Move figure back to original position.
set(fig,'Position',fig_pos)

% Build RGB image, set background pixels to NaN and put it in 'CData' for
% the uicontrol.
if isempty(frame.colormap)
	rgb = frame.cdata;
else
	rgb = reshape(frame.colormap(frame.cdata,:),[pos([4,3]),3]);
end
size_rgb = size(rgb);
rgb = double(rgb)/255;
back = repmat(permute(s.BackgroundColor,[1 3 2]),size_rgb(1:2));
isback = all(rgb == back,3);
rgb(repmat(isback,[1 1 3])) = NaN;
set(h,'CData',rgb,'String','','Visible',s.Visible)

% Assign output argument if necessary.
if nargout
	hout = h;
end