root/DPWavesProc/trunk/DPWavesProc/test_code/specmultiplot_compare.m

function [Winfo]=specmultiplot_compare(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','IMLM UVW','EMEP range','IMLM range','EMEP radial','IMLM radial','wavesmon'}';
Winfo.hsig=[];
Winfo.hconf=[];
Winfo.peakP=[];
Winfo.dirP=[];
Winfo.Ddir=[];
Winfo.Spectrum.EMEPuvw=[];
Winfo.Spectrum.IMLMuvw=[];
Winfo.Spectrum.EMEPrange=[];
Winfo.Spectrum.IMLMrange=[];
Winfo.Spectrum.EMEPradial=[];
Winfo.Spectrum.IMLMradial=[];
Winfo.Spectrum.wavesmon=[];

%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]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now for IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_uvw_F01_D2, EPout]=dirspec(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]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now with IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_range_F01_D2, EPout]=dirspec(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]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0});

    %now with IMLM
    EP.method='IMLM';
    EP.iter=3;
    [I_radial_F01_D2, EPout]=dirspec(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_compare(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_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;
    
    %Make sure wavesmon spectrum is formatted right before outputting, change to m^2/Hz/deg
    wmon.S=spec/(360*1000*1000);
    
    %Fill up the Winfo data structure

    hsig=[vertcat(E_UVW_WI.hsig,I_UVW_WI.hsig,E_range_WI.hsig,I_range_WI.hsig,E_radial_WI.hsig,I_radial_WI.hsig,wmon_WI.hsig)];
    hconf=[vertcat(E_UVW_WI.hconf,I_UVW_WI.hconf,E_range_WI.hconf,I_range_WI.hconf,E_radial_WI.hconf,I_radial_WI.hconf,[NaN NaN])];
    peakP=[vertcat(E_UVW_WI.tp,I_UVW_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,I_UVW_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,I_UVW_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.hconf=[horzcat(Winfo.hconf,hconf)];
    Winfo.peakP=[horzcat(Winfo.peakP,peakP)];
    Winfo.dirP=[horzcat(Winfo.dirP,dirP)];
    Winfo.Ddir=[horzcat(Winfo.Ddir,Ddir)];
    Winfo.Spectrum.EMEPuvw=[cat(3,Winfo.Spectrum.EMEPuvw,E_uvw_F01_D2.S)];
    Winfo.Spectrum.IMLMuvw=[cat(3,Winfo.Spectrum.IMLMuvw,I_uvw_F01_D2.S)];
    Winfo.Spectrum.EMEPrange=[cat(3,Winfo.Spectrum.EMEPrange,E_range_F01_D2.S)];
    Winfo.Spectrum.IMLMrange=[cat(3,Winfo.Spectrum.IMLMrange,I_range_F01_D2.S)];
    Winfo.Spectrum.EMEPradial=[cat(3,Winfo.Spectrum.EMEPradial,E_radial_F01_D2.S)];
    Winfo.Spectrum.IMLMradial=[cat(3,Winfo.Spectrum.IMLMradial,I_radial_F01_D2.S)];
    Winfo.Spectrum.wavesmon=[cat(3,Winfo.Spectrum.wavesmon,wmon.S)];
    
end

%now create the time series plots of the wave info
[fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot_compare(Winfo);

saveas(fig_sigh,['sigh_' dateinfo '.fig']);
saveas(fig_pp,['peakp_' dateinfo '.fig']);
saveas(fig_dp,['dirpeak_' dateinfo '.fig']);
saveas(fig_dd,['domdir_' dateinfo '.fig']);

close all;