function [Winfo]=specmultiplot(dateinfo1,dateinfo2,timestep); 
%where dateinfo1 and 2 are the strings of the first and last date extension
%of the pressure,range,etc files and timestep is the time in hours between
%each burst

%This function takes the pressure, range, orbit and sysinfo files created
%in python and outputs a Winfo data structure with the sig wave height,
%peak period, direction of peak period and dominant dir.  It also generates
%polar plots and 2d freq and directional plots of a variety of different
%estimation methods, resolutions and the wavesmon output. 


%first set up the time vector going from dateinfo1 to dateinfo2
Winfo.time=[];

date1=datevec(dateinfo1, 'yymmddHHMM');
date2=datevec(dateinfo2,'yymmddHHMM');

if date1 == date2
    Winfo.time=[datenum(date1)];
else
    Winfo.time=[datenum(date1)];
    newdate=date1;
    while datenum(newdate) ~= datenum(date2)
        newdate=newdate+[0 0 0 timestep 0 0];
        Winfo.time=[horzcat(Winfo.time,datenum(newdate))];
    end
end
            
      
%set up time, which is Winfo.time as date str in the yymmddHHMM format
time=datestr(Winfo.time, 'yymmddHHMM');

%set up data structure for wave info
Winfo.setup={'EMEP UVW','EMEP UVW F005','EMEP UVW D01','IMLM UVW','IMLM UVW F005','IMLM UVW D01','EMEP range','IMLM range','EMEP radial','IMLM radial','wavesmon'}';
Winfo.hsig=[];
Winfo.peakP=[];
Winfo.dirP=[];
Winfo.Ddir=[];

%Load the data and run the script
for i=1:length(time(:,1))
    dateinfo=time(i,:);

    pressure=strcat('pressure_',dateinfo,'.txt');
    range=strcat('range_',dateinfo,'.txt');
    orbit=strcat('orbit_',dateinfo,'.txt');
    sysinfo=strcat('sysinfo_',dateinfo,'.txt');
    spec=strcat('DSpec',dateinfo,'.txt');

    pressure=load(pressure);
    range=load(range);
    orbit=load(orbit);
    sysinfo=load(sysinfo);
    spec=load(spec);

    %set up data with uvw and pressure, freq at 0.01 and dir at 2
    [ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,1);

    %run diwasp to generate this spectrum
    [E_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now for IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now for IMLM with increased Freq resolution
    SM.freqs=[0.005:0.005:0.4];
    [I_uvw_F005_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now for IMLM with increased Dir resolution
    SM.freqs=[0.01:0.01:0.4];
    SM.dirs=[-180:1:180];
    [I_uvw_F01_D1, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now for EMEP with increased Freq res and then Dir res
    EP.method='EMEP';
    EP.iter=100;
    SM.freqs=[0.005:0.005:0.4];
    SM.dirs=[-180:2:180];
    [E_uvw_F005_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    SM.freqs=[0.01:0.01:0.4];
    SM.dirs=[-180:1:180];
    [E_uvw_F01_D1, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});


    %set up data with ranges, freq and dir at default
    [ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,2);

    % run diwasp to generate this spectrum with EMEP
    [E_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now with IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %set up data with radial velocities, freq and dir at default
    [ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,3);

    % run diwasp to generate this spectrum with EMEP
    [E_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now with IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %Use waveplot to plot the polar, freq and dir plots of the spectra

    [E_UVW_WI,I_UVW_WI,E_range_WI,I_range_WI,wmon_WI,fig1,fig2]=waveplot(E_uvw_F01_D2,I_uvw_F01_D2,E_range_F01_D2,I_range_F01_D2,spec,sysinfo);
    
    saveas(fig1,['polar1_' dateinfo '.fig']);
    saveas(fig2,['dirfreq1_' dateinfo '.fig']);
    
    [E_UVW_F_WI,E_UVW_D_WI,I_UVW_F_WI,I_UVW_D_WI,fig3,fig4]=waveplot2(E_uvw_F01_D2,I_uvw_F01_D2,E_uvw_F005_D2,I_uvw_F005_D2,E_uvw_F01_D1,I_uvw_F01_D1,sysinfo);

    saveas(fig3,['polar2_' dateinfo '.fig']);
    saveas(fig4,['dirfreq2_' dateinfo '.fig']);
    
    [E_radial_WI,I_radial_WI,fig5,fig6]=radialwaveplot(E_radial_F01_D2,I_radial_F01_D2,spec,sysinfo);

    saveas(fig5,['polar3_' dateinfo '.fig']);
    saveas(fig6,['dirfreq3_' dateinfo '.fig']);
    
    close all;
    
    %Fill up the Winfo data structure

    hsig=[vertcat(E_UVW_WI.hsig,E_UVW_F_WI.hsig,E_UVW_D_WI.hsig,I_UVW_WI.hsig,I_UVW_F_WI.hsig,I_UVW_D_WI.hsig,E_range_WI.hsig,I_range_WI.hsig,E_radial_WI.hsig,I_radial_WI.hsig,wmon_WI.hsig)];
    peakP=[vertcat(E_UVW_WI.tp,E_UVW_F_WI.tp,E_UVW_D_WI.tp,I_UVW_WI.tp,I_UVW_F_WI.tp,I_UVW_D_WI.tp,E_range_WI.tp,I_range_WI.tp,E_radial_WI.tp,I_radial_WI.tp,wmon_WI.tp)];
    dirP=[vertcat(E_UVW_WI.dtp,E_UVW_F_WI.dtp,E_UVW_D_WI.dtp,I_UVW_WI.dtp,I_UVW_F_WI.dtp,I_UVW_D_WI.dtp,E_range_WI.dtp,I_range_WI.dtp,E_radial_WI.dtp,I_radial_WI.dtp,wmon_WI.dtp)];
    Ddir=[vertcat(E_UVW_WI.dp,E_UVW_F_WI.dp,E_UVW_D_WI.dp,I_UVW_WI.dp,I_UVW_F_WI.dp,I_UVW_D_WI.dp,E_range_WI.dp,I_range_WI.dp,E_radial_WI.dp,I_radial_WI.dp,wmon_WI.dp)];
    
    Winfo.hsig=[horzcat(Winfo.hsig,hsig)];
    Winfo.peakP=[horzcat(Winfo.peakP,peakP)];
    Winfo.dirP=[horzcat(Winfo.dirP,dirP)];
    Winfo.Ddir=[horzcat(Winfo.Ddir,Ddir)];
    
end

% change Winfo.time back to a matlab datenum

%Winfo.time=datenum(Winfo.time);