%
%   ED96.M   Plots data for Hurricane Edouard (late Aug 1996) using Matlab
%
%   This is a demo that illustrates some aspects of programming within the 
%   Matlab environment.  The goal here is simply to plot some data that were 
%   downloaded from the Web related to hurricane Edouard (which started 
%   Aug 21 1996).  This hurricane came up just off the east coast and caused  
%   some minor damage (high winds, heavy rain, local flooding, etc.) in the 
%   local area.
%
%   The Matlab functions for reading mixed text and numerical data from ascii 
%   files are not particularly easy to use (especially for the beginner).  Therefore, 
%   instead of reading the raw data directly, some editing of the original file 
%   was performed.  Thus, there are two files with basically the same data (the 
%   second file will be used here).
%     1. raw data file:  edouard96.dat (http://weather.unisys.com/hurricane/)
%     2. m-file containing data in modified format:  ed96dat.m
%   The original data was modified with the Wordpad editor...
%
%   A brief study of the raw data shows a correlation of the pressure and wind
%   speed versus time.  Thus two data representations are presented:
%     1. pressure and wind speed versus time 
%          a. two plots on single page -> subplots
%          b. use of double y axes
%     2. wind speed versus pressure
%
%   File prepared by J. R. White, UMass-Lowell (last update: Sept. 2017)
%

      clear all,   close all,   nfig =0;
%
%   First let's execute script file ed96dat.m to get the raw data into Matlab
%   Note: This defines a matrix, A, that has 7 columns (see original and 
%         modified data files), of which only the last 3 columns are used here.
%          col 5 - time of day (in military time so we will need to work on this)
%          col 6 - wind speed in knots
%          col 7 - pressure in millibar
%
      ed96dat           % this defines matrix A (a script file can call another script)
      tt = A(:,5);      % relative time (need to work on this)
      ws = A(:,6)*1.15; % wind speed (mph)
      press = A(:,7);   % pressure (mB) 
%
%   Create a time vector that represents cumulative time
      NT = length(tt);   t = zeros(size(tt));
      for i = 2:NT
        if tt(i) > tt(i-1);   t(i) = t(i-1)+(tt(i)-tt(i-1));   end
        if tt(i) < tt(i-1);   t(i) = t(i-1)+(tt(i)+(24-tt(i-1)));   end
      end
%
%   Plot wind speed and pressure versus time (use subplots)
      nfig = nfig+1;   figure(nfig)
      subplot(2,1,1),plot(t,press,'go','LineWidth',2), grid
      ylabel('Pressure (millibar)')
      r = axis;  r(3) = 900;  r(4) = 1050;  axis(r)  % resets y-axis limits
      title('Ed96: Data for Hurricane Edouard (starting Aug 21, 1996)')
      subplot(2,1,2),plot(t,ws,'rs','LineWidth',2), grid, 
      xlabel('Cumulative Time (hours)'),  ylabel('Wind Speed (mph)')
%
%   Plot wind speed & pressure using two y axes - this is a little advanced
      nfig = nfig+1;   figure(nfig)
      [ax,h1,h2] = plotyy(t,press,t,ws); grid
      xlabel('Cumulative Time (hours)'),
      title('Ed96: Data for Hurricane Edouard (starting Aug 21, 1996)')
      axes(ax(1)), ylabel('Pressure (millibar)'),  set(ax(1),'YColor',[0 0 0])
      set(ax(1),'YLim',[900 1050],'YTick',[900 950 1000 1050], ...
                 'YTickLabel',[900 950 1000 1050]');
      axes(ax(2)), ylabel('Wind Speed (mph)'),     set(ax(2),'YColor',[0 0 0])
      set(ax(2),'YLim',[0 150],'YTick',[0 50 100 150],'YTickLabel',[0 50 100 150]');
      set(h1,'Marker','o','MarkerEdgeColor','g','LineStyle','none','LineWidth',2);
      set(h2,'Marker','s','MarkerEdgeColor','r','LineStyle','none','LineWidth',2);
      legend([h1 h2],'Pressure','Wind Speed','location','Southeast')
%
%   Plot wind speed versus pressure and fit a best straight line through the data
      [coeff] = polyfit(press,ws,1);                % does linear fit
      wsfit = polyval(coeff,press);                 % does polynomial evaluation
      nfig = nfig+1;   figure(nfig)
      plot(press,ws,'go',press,wsfit,'b-','LineWidth',2), grid
      r = axis;  r(3) = 25;  r(4) = 150;  axis(r)    % resets y-axis limits
      title('Ed96: Data for Hurricane Edouard (starting Aug 21, 1996)')
      xlabel('Pressure (millibar)'),ylabel('Wind Speed (mph)')
      legend('Data points','Linear fit')
%
%   end of program
