Index: DPWavesProc/trunk/DPWavesProc/8marrayconf.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/8marrayconf.m (revision 216)
@@ -1,0 +1,3 @@
+%calculate just the frequency energy spectrum
+arrayrawfreq=sum(arrayraw.S')*dirres;
+array8mfreq=sum(array8m.S')*dirres;
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/Hsig.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/Hsig.m (revision 216)
@@ -1,0 +1,20 @@
+function Hs=Hsig(varargin)
+%DIWASP V1.1 function to calculate significant wave height
+%
+%Hs=Hsig(SM)
+%
+%Hs is 4 times the zeroth moment wave height of spectral matrix SM
+%
+%"help data_structures" for information on the DIWASP data structures
+
+%Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth
+
+
+if nargin==1
+    SM=varargin{1};
+    df=SM.freqs(2)-SM.freqs(1);ddir=abs(SM.dirs(2)-SM.dirs(1));S=SM.S;
+elseif nargin==3
+    df=varargin{2};ddir=varargin{3};S=varargin{1};
+end
+
+Hs=sqrt(16*sum(sum(real(S)*df*ddir)));
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/HsigConf.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/HsigConf.m (revision 216)
@@ -1,0 +1,29 @@
+function HsConf=HsigConf(SM)
+
+%This is a modification of the DIWASP Hsig function used to calculate the
+%95% confidence bounds for the significant wave height  
+%
+%Added on 9/17/08 by G. Dusek
+
+df=SM.freqs(2)-SM.freqs(1);
+ddir=abs(SM.dirs(2)-SM.dirs(1));
+S=SM.S;
+numFreq=size(SM.freqs,1);
+
+
+%the sig wave height degF will be different since we are summing over all of the
+%frequency bands.
+hsigDegF=SM.degF*numFreq;
+
+%Calculate the lower and upper confidence coefficients
+chiUp=chi2inv(.975,hsigDegF);
+chiLow=chi2inv(.025,hsigDegF);
+
+coeffUp=hsigDegF/chiLow;
+coeffLow=hsigDegF/chiUp;
+
+%Calculate the upper and low conf bounds on the sig wave height
+confUp=sqrt(coeffUp*16*sum(sum(real(S)*df*ddir)));
+confLow=sqrt(coeffLow*16*sum(sum(real(S)*df*ddir)));
+HsConf=[confLow,confUp];
+
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/Winfo_plot1.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/Winfo_plot1.m (revision 216)
@@ -1,0 +1,48 @@
+function [fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot1(Winfo);
+
+%This function takes the wave info (sigH, peak period, direction of peak and dominant direction
+%for a series of bursts and plots it as a function of time
+
+
+%plot sig wave height
+scrsz = get(0,'ScreenSize');
+fig_sigh=figure('Position',[scrsz]);
+h1=plot(Winfo.time,Winfo.hsig(1,:),'b.-');
+
+datetick('x',15,'keepticks');
+legend([h1],'EMEP radial','location','best');
+title(['Significant Wave Height starting on ',datestr(Winfo.time(:,1),0)]);
+xlabel('time');
+ylabel('Wave Height (m)');
+
+%plot peak period
+fig_pp=figure('Position',[scrsz]);
+h1=plot(Winfo.time,Winfo.peakP(1,:),'b.-');
+
+datetick('x',15,'keepticks');
+legend([h1],'EMEP radial','location','best');
+title(['Peak Period starting on ',datestr(Winfo.time(:,1),0)]);
+xlabel('time');
+ylabel('Period (seconds)');
+
+%plot dir of peak period
+fig_dp=figure('Position',[scrsz]);
+h1=plot(Winfo.time,Winfo.dirP(1,:),'b.-');
+
+datetick('x',15,'keepticks');
+legend([h1],'EMEP radial','location','best');
+title(['Direction of Peak Period starting on ',datestr(Winfo.time(:,1),0)]);
+xlabel('time');
+ylabel('Direction (Degrees True)');
+
+%plot the dominant direction
+
+fig_dd=figure('Position',[scrsz]);
+h1=plot(Winfo.time,Winfo.Ddir(1,:),'b.-');
+
+datetick('x',15,'keepticks');
+legend([h1],'EMEP radial','location','best');
+title(['Dominant Direction starting on ',datestr(Winfo.time(:,1),0)] );
+xlabel('time');
+ylabel('Direction (Degrees True)');
+
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/confplot.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/confplot.m (revision 216)
@@ -1,0 +1,99 @@
+function [fig1]=confplot(EMEP,rangeE,radialE,spec,sysinfo)
+
+%set up the wavesmon data in a structure
+
+%what is the magnetic variation to nearest degree
+magvar=10;
+
+%first change to m^2/Hz/deg
+wmon.S=spec/(360*1000*1000);
+
+%what is the start angle
+heading=sysinfo(18,:);
+heading=heading/100;
+sangle=heading+magvar;
+
+%what is the frequency and dir resolution for those generated in DIWASP
+freqres=0.01;
+freqs=[0.01:0.01:.4];
+dirres=2;
+dirs=[-180:2:180];
+
+%set up the directions
+adj_angle=90-sangle+360+180;
+wmon.dirs=[adj_angle:-4:-(356-adj_angle)]; 
+wmondirres=4;
+wmon.xaxisdir=90;
+wmon.freqs=[0.00781250:0.00781250:1];
+wmonfreqres=0.00781250;
+
+%calculate just the frequency energy spectrum
+EMEPfreq=sum(EMEP.S')*dirres;
+EMEPrangefreq=sum(rangeE.S')*dirres;
+EMEPradialfreq=sum(real(radialE.S)')*dirres;
+wmonfreq=sum(wmon.S')*wmondirres;
+
+%Compute the coefficient for the upper and lower error bounds for the power
+%spectrum assuming 95% confidence.  Added on 9/17/08
+degF=radialE.degF;
+chiUp=chi2inv(.975,degF);
+chiLow=chi2inv(.025,degF);
+coeffUp=degF/chiLow;
+coeffLow=degF/chiUp;
+logCoeffUp=log10(coeffUp);
+logCoeffLow=log10(coeffLow);
+
+%calculate the conf limits throughout the frequency spectrum
+EMEPradialfreqUP=EMEPradialfreq*coeffUp;
+EMEPradialfreqLOW=EMEPradialfreq*coeffLow;
+
+%same thing for the Puvw and range
+degF2=25;
+chiUp2=chi2inv(.975,degF2);
+chiLow2=chi2inv(.025,degF2);
+coeffUp2=degF2/chiLow2;
+coeffLow2=degF2/chiUp2;
+logCoeffUp2=log10(coeffUp2)
+logCoeffLow2=log10(coeffLow2)
+
+%calculate the conf limits throughout the frequency spectrum
+EMEPfreqUP=EMEPfreq*coeffUp2;
+EMEPfreqLOW=EMEPfreq*coeffLow2;
+EMEPrangefreqUP=EMEPrangefreq*coeffUp2;
+EMEPrangefreqLOW=EMEPrangefreq*coeffLow2;
+
+scrsz = get(0,'ScreenSize');
+fig1=figure('Position',[scrsz]);
+
+%plot the frequency energy spectrum
+h1=semilogy(freqs,EMEPfreq,'b');
+hold on
+plot([0.1,0.1],[10^logCoeffLow2,10^logCoeffUp2],'b.-');
+%plot(freqs,EMEPfreqUP,'b--');
+%plot(freqs,EMEPfreqLOW,'b--');
+
+h2=semilogy(freqs,EMEPrangefreq,'r');
+plot([0.11,0.11],[10^logCoeffLow2,10^logCoeffUp2],'r.-');
+%plot(freqs,EMEPrangefreqUP,'r--');
+%plot(freqs,EMEPrangefreqLOW,'r--');
+
+h3=semilogy(freqs,EMEPradialfreq,'g');
+plot([0.12,0.12],[10^logCoeffLow,10^logCoeffUp],'g.-');
+%plot(freqs,EMEPradialfreqUP,'g--');
+%plot(freqs,EMEPradialfreqLOW,'g--');
+ylim([10^-2,10^1]);
+%set(gca,'ytick',[-4 -3 -2 -1 0 1]);
+axis(axis);
+h4=semilogy(wmon.freqs,wmonfreq,'k');
+legend([h1,h2,h3,h4],'EMEP uvw','EMEP range','EMEP radial','wavesmon','location','best');
+title('directional wave spectrum integrated over direction');
+xlabel('frequency in Hz');
+ylabel('m^2 / hz');
+
+
+
+
+
+
+
+
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot1.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot1.m (revision 216)
@@ -1,0 +1,159 @@
+function [Winfo,spec2d,fd,thetad,td]=specmultiplot1(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]=waveplottest1(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]);
+
+%this is added in to format the data for Xwaves (using radial velocities)
+
+spec2d=Winfo.Spectrum.EMEPradial;
+spec2d=spec2d(:,:,1:180);
+%firstspec=spec2d(:,:,91:180);
+%nextspec=spec2d(:,:,1:90);
+%spec2d=cat(3,firstspec,nextspec);
+spec2d=spec2d*(360/(2*pi));
+
+%this changes the coordinates from axis to compass
+newspec2d=zeros(size(spec2d));
+for i= 1:size(spec2d,1)
+    oldspec=squeeze(spec2d(i,:,:));
+
+    newspec=zeros(size(oldspec));
+
+    for j = 1:46
+        newspec(:,j)=oldspec(:,47-j);
+    end
+    for j = 47:180
+        newspec(:,j)=oldspec(:,181-(j-46));
+    end
+    
+    %firstspec=newspec(:,91:180);
+    %secondspec=newspec(:,1:90);
+    %newspec=cat(2,firstspec,secondspec);
+    
+    newspec2d(i,:,:)=newspec;    
+end
+spec2d=newspec2d;
+
+
+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_plot1(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;
+
+
+
+
Index: DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplottest1.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplottest1.m (revision 216)
@@ -1,0 +1,130 @@
+function [E_radial_WI,fig1,fig2]=waveplottest1(radialE,sysinfo) 
+
+%make sure to load the EMEP, IMLM and wavesmon samples
+%EMEP=load('emep'), IMLM=load('imlm'), spec=load('wavesmon')
+%rangeE=load('EMEP_range'), rangeI=load('IMLM_range')
+
+%set up the wavesmon data in a structure
+
+%what is the magnetic variation to nearest degree
+magvar=10;
+
+%what is the start angle
+heading=sysinfo(18,:);
+heading=heading/100;
+sangle=heading+magvar;
+
+%what is the frequency and dir resolution for those generated in DIWASP
+freqres=0.01;
+freqs=[0.01:0.01:.4];
+dirres=2;
+dirs=[-180:2:180];
+
+%set up the directions
+adj_angle=90-sangle+360+180;
+
+% plot the spectrum generated through DIWASP 
+
+scrsz = get(0,'ScreenSize');
+fig1=figure('Position',[scrsz]);
+subplot(1,1,1);
+subplotspec(radialE,4);
+title('radial velocity data');
+
+
+%calculate just the dir energy spectrum on a single graph
+
+
+EMEPradialdir=sum(radialE.S)*freqres;
+
+%calculate just the frequency energy spectrum
+EMEPradialfreq=sum(radialE.S')*dirres;
+
+%Find the maximum for the directional spectrum so we can set up the proper
+%x-axis
+[maxvalue,maxindex] = max(EMEPradialdir);
+maxdir=dirs(maxindex);
+% set up the x-axis for all of the spectra depending on the max
+if ((100 < maxdir) | (maxdir < -100));
+    %for diwasp spectra
+    index1=find(dirs < 0);
+    index2=find(dirs > -1);
+    dirs(index1)=dirs(index1) +360;
+    %plot the directional energy spectrum
+    fig2=figure('Position',[scrsz]);
+    subplot(1,2,1);
+    h1 = plot(dirs(index2),EMEPradialdir(index2),'b');
+    hold on
+    h1a= plot(dirs(index1),EMEPradialdir(index1),'b');
+else
+    %for diwasp spectra do nothing
+    
+    %plot the directional energy spectrum
+    fig2=figure('Position',[scrsz]);
+    subplot(1,2,1);
+    h1 = plot(dirs,EMEPradialdir,'b');
+end
+
+legend([h1],'EMEP radial','location','best');
+title('directional wave spectrum integrated over frequency');
+xlabel('axis angle (degrees true)');
+ylabel('m^2 / deg');
+
+%plot the frequency energy spectrum
+subplot(1,2,2);
+plot(freqs,EMEPradialfreq,'b');
+legend('EMEP radial','location','best');
+title('directional wave spectrum integrated over direction');
+xlabel('frequency in Hz');
+ylabel('m^2 / hz');
+
+% ______Calculate and display the wave parameters SigH, Tp, Dp,DTp_______
+
+
+% for Radial of EMEP
+
+% Calculate the Sig wave height
+EMEP_radial_Hsig=Hsig(radialE);
+%Use the function HsigConf.m to calculate the sigH confidence limits
+EMEP_radial_HsConf=HsigConf(radialE);
+
+% Calculate the peak period Tp
+[P,I]=max(EMEPradialfreq);
+EMEP_radial_Tp=1/(freqs(I));
+
+%Calculate the Direction of the peak period DTp
+[P,I]=max(real(radialE.S(I,:)));
+EMEP_radial_DTp=dirs(I);
+
+%Calculate the Dominant Direction Dp
+[P,I]=max(EMEPradialdir);
+EMEP_radial_Dp=dirs(I);
+
+%Display on the screen the SigH,Tp,Dp,DTp
+disp(['EMEP radial']);
+disp(['SigH (meters): ' num2str(EMEP_radial_Hsig)]);
+disp(['SigH 95% confidence limits: ' num2str(EMEP_radial_HsConf)]);
+disp(['peak period (seconds): ' num2str(EMEP_radial_Tp)]);
+disp(['Dir of peak period: ' num2str(compangle(EMEP_radial_DTp, radialE.xaxisdir))]);
+disp(['Dominant Direction: ' num2str(compangle(EMEP_radial_Dp, radialE.xaxisdir))]);
+disp(['  ']);
+
+E_radial_WI.hsig=EMEP_radial_Hsig;
+E_radial_WI.hconf=EMEP_radial_HsConf;
+E_radial_WI.tp=EMEP_radial_Tp;
+E_radial_WI.dtp=compangle(EMEP_radial_DTp, radialE.xaxisdir);
+E_radial_WI.dp=compangle(EMEP_radial_Dp, radialE.xaxisdir);
+
+
+
+%function to change from axis angles to compass bearings
+
+function angle=compangle(angle,xaxisdir)
+angle=xaxisdir*ones(size(angle))-angle;
+angle=angle+360*(angle<0);
+angle=angle-360*(angle>360);
+
+
+
+
+
Index: DPWavesProc/trunk/DPWavesProc/compangle.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/compangle.m (revision 216)
@@ -1,0 +1,6 @@
+%function to change from axis angles to compass bearings
+
+function angle=compangle(angle,xaxisdir)
+angle=xaxisdir*ones(size(angle))-angle;
+angle=angle+360*(angle<0);
+angle=angle-360*(angle>360);
Index: DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/private/HsigConf.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/diwasp_1_1GD/private/HsigConf.m (revision 216)
@@ -1,0 +1,29 @@
+function HsConf=HsigConf(SM)
+
+%This is a modification of the DIWASP Hsig function used to calculate the
+%95% confidence bounds for the significant wave height  
+%
+%Added on 9/17/08 by G. Dusek
+
+df=SM.freqs(2)-SM.freqs(1);
+ddir=abs(SM.dirs(2)-SM.dirs(1));
+S=SM.S;
+numFreq=size(SM.freqs,1);
+
+
+%the sig wave height degF will be different since we are summing over all of the
+%frequency bands.
+hsigDegF=SM.degF*numFreq;
+
+%Calculate the lower and upper confidence coefficients
+chiUp=chi2inv(.975,hsigDegF);
+chiLow=chi2inv(.025,hsigDegF);
+
+coeffUp=hsigDegF/chiLow;
+coeffLow=hsigDegF/chiUp;
+
+%Calculate the upper and low conf bounds on the sig wave height
+confUp=sqrt(coeffUp*16*sum(sum(real(S)*df*ddir)));
+confLow=sqrt(coeffLow*16*sum(sum(real(S)*df*ddir)));
+HsConf=[confLow,confUp];
+
Index: DPWavesProc/trunk/DPWavesProc/nortek/nortek_HsigConf.m
===================================================================
---  (revision )
+++ DPWavesProc/trunk/DPWavesProc/nortek/nortek_HsigConf.m (revision 216)
@@ -1,0 +1,29 @@
+function HsConf=nortek_HsigConf(SM)
+
+%This is a modification of the DIWASP Hsig function used to calculate the
+%95% confidence bounds for the significant wave height  
+%
+%Added on 9/17/08 by G. Dusek
+
+df=SM.freqs(2)-SM.freqs(1);
+ddir=abs(SM.dirs(2)-SM.dirs(1));
+S=SM.S;
+numFreq=size(SM.freqs,1);
+
+
+%the sig wave height degF will be different since we are summing over all of the
+%frequency bands.
+hsigDegF=2*SM.degF*numFreq;
+
+%Calculate the lower and upper confidence coefficients
+chiUp=chi2inv(.975,hsigDegF);
+chiLow=chi2inv(.025,hsigDegF);
+
+coeffUp=hsigDegF/chiLow;
+coeffLow=hsigDegF/chiUp;
+
+%Calculate the upper and low conf bounds on the sig wave height
+confUp=sqrt(coeffUp*16*sum(sum(real(S)*df*ddir)));
+confLow=sqrt(coeffLow*16*sum(sum(real(S)*df*ddir)));
+HsConf=[confLow,confUp];
+