root/DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot.m

function [Winfo,spec2d,fd,thetad,td]=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
            
%Make a directory for the data
mkdir('tempdir');

%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 radial'};
Winfo.hsig=[];
Winfo.hconf=[];
Winfo.peakP=[];
Winfo.dirP=[];
Winfo.Ddir=[];
Winfo.Spectrum.EMEPradial=[];

%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');
 

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

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

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

    [E_radial_WI,fig1,fig2]=waveplot(radialE,sysinfo); 
    
    saveas(fig1,['polar1_' dateinfo '.fig']);
    saveas(fig2,['dirfreq1_' dateinfo '.fig']);
    
    movefile(['polar1_' dateinfo '.fig'],'tempdir');
    movefile(['dirfreq1_' dateinfo '.fig'],'tempdir');
    close all;
    
    
    %Fill up the Winfo data structure

    hsig=[vertcat(E_radial_WI.hsig)];
    hconf=[vertcat(E_radial_WI.hconf)];
    peakP=[vertcat(E_radial_WI.tp)];
    dirP=[vertcat(E_radial_WI.dtp)];
    Ddir=[vertcat(E_radial_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.EMEPradial=[cat(3,Winfo.Spectrum.EMEPradial,radialE.S)];
end

Winfo.Spectrum.EMEPradial=permute(Winfo.Spectrum.EMEPradial,[3,1,2]);

spec2d=Winfo.Spectrum.EMEPradial;
%account for if there are 2 (just one sample) or 3 (time series) dimensions
numdimensions=ndims(spec2d);
if numdimensions == 3
    spec2d=spec2d(:,:,1:180);
else
    spec2d=spec2d(:,1:180);
end

%this changes the coordinates from axis to compass (and from direction TO,
%to direction FROM)
% NEED TO ADJUST THIS TO AUTOMATICALLY ACOUNT FOR USER INPUTTED SM.dirs!!!
newspec2d=zeros(size(spec2d));
for i= 1:size(spec2d,1)
    if numdimensions == 3
        oldspec=squeeze(spec2d(i,:,:));
    else
        oldspec=spec2d;
    end

    newspec=zeros(size(oldspec));

    %this is the NEW way to organize the code for SM.dirs=[0:2:358] (edited
    %on 3/27/09
    for j = 1:136
        newspec(:,j)=oldspec(:,137-j);
    end
    for j = 137:180
        newspec(:,j)=oldspec(:,181-(j-136));
    end
    
    if numdimensions == 3
        newspec2d(i,:,:)=newspec; 
    else
        newspec2d=newspec;
    end   
end

Winfo.Spectrum.EMEPradial=newspec2d;

%Use this to organize spec IF SM.dirs=[-180:2:178]
    %for j = 1:46
    %    newspec(:,j)=oldspec(:,47-j);
    %end
    %for j = 47:180
    %    newspec(:,j)=oldspec(:,181-(j-46));
    %end

%*************************************************************************
%these following variables can be output if you wish to continue to process the data through
%the XWAVES software.

spec2d=newspec2d*(360/(2*pi));
fd=[0.01:0.01:0.4];
thetad=[0:2:358];
date1=datevec(Winfo.time);
td=date1(:,1:5);

save(['specdata_' dateinfo],'Winfo','spec2d','fd','thetad','td');
movefile(['specdata_' dateinfo '.mat'],'tempdir');

%**************************************************************************

%now create the time series plots of the wave info
[fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot(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']);

movefile(['sigh_' dateinfo '.fig'],'tempdir');
movefile(['peakp_' dateinfo '.fig'],'tempdir');
movefile(['dirpeak_' dateinfo '.fig'],'tempdir');
movefile(['domdir_' dateinfo '.fig'],'tempdir');

movefile('tempdir',strcat('procdata_',dateinfo));

close all;