%
%   OILTANK.M   Modeling of Flow from a Damaged Oil Tank 
%
%   This demo models and evaluates the flow of oil from a horizontal tank.  This
%   simulation requires several numerical procedures including:
%     1. numerical integration of a fairly complicated function (QUADL)
%     2. spline fitting, root finding, and piecewise polynomial interpolation (SPLINE,
%        FZERO, and PPVAL)
%     3. numerical solution of a first order differential equation (ODE23)
%   
%   The cross sectional area of the horizontal cylindrical tank is given by
%             A(h) = 2*int[x(y)dy]  for 0 < y < h
%   where     x(y) = sqrt(Dy - y^2)
%   This integral is performed numerically for several values of h using QUADL.
%   This gives a relationship for A(h) and V(h) = L*A(h).
%
%   The volume of oil in the tank, V(h) = L*A(h), is then fit to a set of piecewise
%   cubic polynomials using SPLINE.  Then, knowing the initial volume of oil, one can 
%   determine the initial height of oil in the tank, ho.  This is done using FZERO
%   to find the root of f(h) = V(h) - Vo = 0.  PPVAL is also used here to perform the 
%   piecewise polynomial interpolation - evaluating V(h) for arbitrary h.
%
%   With the initial height known, ODE23 is used to numerically integrate the 
%   differential equation that represents a mass balance on the oil in the tank 
%   [written in terms of h(t)].  The equation of interest here is
%             dh/dt = -(Aout/2L)*sqrt(2g/(D-h))    with h(0) = ho   
%
%   Finally, now that h(t) is known, PPVAL is used again to determine V(t) = V(h(t))
%   and to evaluate the amount of oil that is spilled from the tank, Vs(t) = Vo-V(t).
%
%   This main m-file also refers to the following function files:
%      oiltanka.m  -  function file for QUADL routine
%      oiltankb.m  -  function file for FZERO routine
%      oiltankc.m  -  function file for ODE23 routine
%
%   File prepared by J. R. White, UMass-Lowell (July 2003)
%

%
%   getting started
      clear all, close all, nfig = 0;
      global PP Vo
%
%   specify problem data and convert to consistent units
      DIAM = 2.1;                  % tank diameter (m)
      LENGTH = 5.6;                % tank length (m)
      Vo = 5000;                   % initial volume of oil (gallons)
      cg2m3 = 3.785*0.001;         % conversion factor (gallons to m^3)
      Vo = Vo*cg2m3;               % initial volume of oil (gallons converted to m^3)
      Dout = 1.0;                  % diameter of opening (inches)
      cin2m = 2.54/100;            % conversion factor (inches to m)
      Dout = Dout*cin2m;           % diameter of opening (inches converted to m)
      AOUT = pi*Dout*Dout/4;       % area of opening (m^2)
      GACCEL = 9.8;                % gravitational accel (m^2/sec)
      Tresp = 50;                  % response time (minutes)
      Tresp = Tresp*60;            % response time (minutes converted to seconds)
      hmin = 0;   hmax = DIAM;     % min & max value for height of oil (m)
%
%   find volume of oil versus height of oil 
      Nh1 = 11;   hd1 = linspace(hmin,hmax,Nh1);   Vh1 = zeros(size(hd1));
      for i = 2:Nh1
        Vh1(i) = 2*LENGTH*quadl('oiltanka',hmin,hd1(i),[],[],DIAM);
      end 
%
%   find piecewise cubic polynomial fit to discrete data
      PP = spline(hd1,Vh1);
%
%   plot oil volume versus oil height (over fine grid for smooth curve)
      Nh2 = 51;   hd2 = linspace(hmin,hmax,Nh2);   Vh2 = ppval(PP,hd2);
      Vhg = Vh2/cg2m3;      % oil volume in gallons
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(hd2,Vhg,'b-','LineWidth',2),grid
      range = axis;   range(2)= hmax;   axis(range);
      title('OilTank:  Volume of Oil Versus Oil Height')                               
      xlabel('Oil Height (m)'),ylabel('Oil Volume (gallons)')
%
%   knowing the initial oil volume, Vo, determine the initial height, ho
      ho = fzero('oiltankb',hmax);
      disp('Initial oil height (m) = '),   ho
%
%   now, with ho known, solve the differential balance equation for h(t)
      to = 0;   tf = 120*60;  tol = 0.00001;
      options = odeset('RelTol',tol);
      [t,h] = ode23('oiltankc',[to tf],ho,options,AOUT,LENGTH,GACCEL,DIAM);
%
%   plot h(t)
      subplot(2,2,3),plot(t/60,h,'b-','LineWidth',2),grid
      range = axis;   range(2)= 120;   axis(range);
      title('Height of Oil Versus Time')                               
      xlabel('Time (min)'),ylabel('Oil Height (m)')
%
%   now find and plot amount of oil spilled
      Vs = Vo-ppval(PP,h);
      Vsg = Vs/cg2m3;      % spilled oil volume in gallons
      subplot(2,2,4),plot(t/60,Vsg,'b-','LineWidth',2),grid
      range = axis;   range(2)= 120;   axis(range);
      title('Volume of Oil Spilled Versus Time')                               
      xlabel('Time (min)'),ylabel('Oil Spilled (gallons)')
%
      range = axis;
      xt = range(1) + 0.40*(range(2)-range(1));
      yt1 = range(3) + 0.20*(range(4)-range(3)); 
      yt2 = range(3) + 0.08*(range(4)-range(3));
      text(xt,yt1,['V_o = ',num2str(Vo/cg2m3),' gallons'])
      text(xt,yt2,['D_{out} = ',num2str(Dout/cin2m),' inches'])
%
%  end simulation