function NewtonMechanics

  % NewtonMechanics visualizes the trajectory of the solution to the
  % problem m*(d^2x/dt^2) = -f(dx/dt) where f(x)=-alpha-beta*x denotes the
  % friction force and m the mass of the considered object. The solution
  % consists of the position x and its velocity dx/dt of the object at
  % time t. The trajectory itself is determined by the ode solver ode45;
  % the direction field (velocity vector field) corresponding to the
  % system is generated by the quiver function. The parameters v0 (initial
  % velocity), m (mass) and alpha,beta (friction) can be varied freely.

  % ----- Initialize setting -------------------------------------------- %
  set(0,'defaulttextinterpreter','latex');
  set(0,'DefaultLineLineWidth',2);
  clear vars; close all; clc;
  % ----- Initialize data ----------------------------------------------- %
  alpha = 1;
  beta = 1;
  mass = 10;
  velo = 10;
  % ----- Create a figure and axes -------------------------------------- %
  figure('units','pixels','position',[100 100 600 450])
  axes('units','pixels','position',[50 120 500 300],'xlim',[0 25],      ...
    'ylim',[0 12]);
  xlabel('position $x(t)$'); ylabel('velocity $\dot y(t)$');
  y = SolveOde(alpha,beta,mass,velo);
  DrawSolution(alpha,beta,mass,y);
  % ----- Create SliderAlpha -------------------------------------------- %
  SliderAlpha = uicontrol('Style','slider','Min',0,'Max',10,'Value',    ...
    alpha,'Position',[440 20 120 20],'Callback',@ValueAlpha);
  % ----- Create EditAlpha ---------------------------------------------- %
  EditAlpha = uicontrol('Style','edit','String',sprintf('%1.4f',alpha), ...
    'Position',[370 20 60 20]);
  % ----- Create TextAlpha ---------------------------------------------- %
  TextAlpha = uicontrol('Style','text','String','alpha','Position',     ...
    [320 20 40 20],'Backgroundcolor',[.8 .8 .8]);
  % ----- Create SliderBeta --------------------------------------------- %
  SliderBeta = uicontrol('Style','slider','Min',0,'Max',10,'Value',     ...
    beta,'Position',[440 50 120 20],'Callback',@ValueBeta);
  % ----- Create EditBeta ----------------------------------------------- %
  EditBeta = uicontrol('Style','edit','String',sprintf('%1.4f',beta),   ...
    'Position',[370 50 60 20]);
  % ----- Create TextBeta ----------------------------------------------- %
  TextBeta = uicontrol('Style','text','String','beta','Position',     ...
    [320 50 40 20],'Backgroundcolor',[.8 .8 .8]);
  % ----- Create SliderMass --------------------------------------------- %
  SliderMass = uicontrol('Style','slider','Min',0,'Max',10,'Value',     ...
    mass,'Position',[140 50 120 20],'Callback',@ValueMass);
  % ----- Create EditMass ----------------------------------------------- %
  EditMass = uicontrol('Style','edit','String',sprintf('%1.4f',mass),   ...
    'Position',[70 50 60 20]);
  % ----- Create TextMass ----------------------------------------------- %
  TextMass = uicontrol('Style','text','String','mass','Position',       ...
    [20 50 40 20],'Backgroundcolor',[.8 .8 .8]);
  % ----- Create SliderVelo --------------------------------------------- %
  SliderVelo = uicontrol('Style','slider','Min',0,'Max',10,'Value',     ...
    velo,'Position',[140 20 120 20],'Callback',@ValueVelo);
  % ----- Create EditVelo ----------------------------------------------- %
  EditVelo = uicontrol('Style','edit','String',sprintf('%1.4f',velo),   ...
    'Position',[70 20 60 20]);
  % ----- Create TextVelo ----------------------------------------------- %
  TextVelo = uicontrol('Style','text','String','velo','Position',       ...
    [20 20 40 20],'Backgroundcolor',[.8 .8 .8]);
% ----- Create ValueAlpha --------------------------------------------- %
  function ValueAlpha(~,~)
    alpha = get(SliderAlpha,'Value');
    beta = get(SliderBeta,'Value');
    mass = get(SliderMass,'Value');
    velo = get(SliderVelo,'Value');
    set(EditAlpha,'String',sprintf('%1.4f',alpha));
    y = SolveOde(alpha,beta,mass,velo);
    DrawSolution(alpha,beta,mass,y);
  end
  % ----- Create ValueBeta ---------------------------------------------- %
  function ValueBeta(~,~)
    alpha = get(SliderAlpha,'Value');
    beta = get(SliderBeta,'Value');
    mass = get(SliderMass,'Value');
    velo = get(SliderVelo,'Value');
    set(EditBeta,'String',sprintf('%1.4f',beta));
    y = SolveOde(alpha,beta,mass,velo);
    DrawSolution(alpha,beta,mass,y);
  end
  % ----- Create ValueMass ---------------------------------------------- %
  function ValueMass(~,~)
    alpha = get(SliderAlpha,'Value');
    beta = get(SliderBeta,'Value');
    mass = get(SliderMass,'Value');
    velo = get(SliderVelo,'Value');
    set(EditMass,'String',sprintf('%1.4f',mass));
    y = SolveOde(alpha,beta,mass,velo);
    DrawSolution(alpha,beta,mass,y);
  end
  % ----- Create ValueVelo ---------------------------------------------- %
  function ValueVelo(~,~)
    alpha = get(SliderAlpha,'Value');
    beta = get(SliderBeta,'Value');
    mass = get(SliderMass,'Value');
    velo = get(SliderVelo,'Value');
    set(EditVelo,'String',sprintf('%1.4f',velo));
    y = SolveOde(alpha,beta,mass,velo);
    DrawSolution(alpha,beta,mass,y);
  end
  % ----- Create ode data ----------------------------------------------- %
  function f = EvalForce(~,y,para)
    f = [y(2);-(1/para(3))*(para(1)+para(2)*y(2)^2)];
  end
  % ----- Solve ode ----------------------------------------------------- %
  function y = SolveOde(alpha,beta,mass,velo)
    [~,y] = ode45(@(t,y) EvalForce(t,y,[alpha,beta,mass]),[0,25],[0,velo]);
  end
  % ----- Visualize results --------------------------------------------- %
  function DrawSolution(alpha,beta,mass,y)
    hold off;
    [t1 t2] = meshgrid(0:1.5:25,0:0.5:12);
    X = zeros(size(t1)); Y = zeros(size(t1)); 
    for i = 1:size(t1,1)
      for j = 1:size(t1,2)
        Z = EvalForce(0,[t1(i,j);t2(i,j)],[alpha,beta,mass]);
        Z = Z/(norm(Z));
        X(i,j) = Z(1);
        Y(i,j) = Z(2);
      end
    end
    quiver(t1,t2,X,Y,0.5);
    hold on;
    plot(y(:,1),y(:,2),'r');
    xlim([0,25]); ylim([0,12]);
    xlabel('position $x(t)$'); ylabel('velocity $\dot x(t)$');
    title('direction field of $m\ddot x=-\alpha-\beta\dot x^2$');
  end
end

% ----------------------------------------------------------------------- %
% Martin Gubisch
% 23.11.2014 --- 19:35