root/adcp/trunk/adcp/nortek/nortek_waveplot.m

function [E_WI,I_WI,nortek_WI,fig1,fig2]=nortek_waveplot(EMEP,IMLM,sysinfo,nortekspec) 

%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 frequency and dir resolution for those generated in DIWASP
freqres=0.01;
freqs=[0.01:0.01:.4];
dirres=2;
dirs=[-180:2:180];

%what is the frequency and dir resolution for those created by nortek
%freqres is the same
nortekfreqs=[0.01:0.01:.4];
nortekdirres=4;
nortekdirs=[270:-4:-86];
nortek.S=nortekspec;
nortek.dirs=nortekdirs;
nortek.freqs=nortekfreqs;
nortek.xaxisdir=90;

% plot the spectrum generated through DIWASP 

scrsz = get(0,'ScreenSize');
fig1=figure('Position',[scrsz]);
subplot(2,2,1);
subplotspec(EMEP,4);
title('EMEP AST beam and radial velocities');

subplot(2,2,2);
subplotspec(IMLM,4);
title('IMLM AST beam and radial velocities');

subplot(2,2,3);
subplotspec(nortek,4);
title('nortek quick wave spectrum');


%calculate just the dir energy spectrum on a single graph

EMEPdir=sum(EMEP.S)*freqres;
IMLMdir=sum(real(IMLM.S))*freqres;
nortekdir=sum(nortek.S)*freqres;


%calculate just the frequency energy spectrum
EMEPfreq=sum(EMEP.S')*dirres;
IMLMfreq=sum(real(IMLM.S)')*dirres;
nortekfreq=sum(nortek.S')*nortekdirres;

%Find the maximum for the directional spectrum so we can set up the proper
%x-axis
[maxvalue,maxindex] = max(EMEPdir);
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;
    %for nortek
    Aindex=find(nortek.dirs < 0);               
    Bindex=find((-1 < nortek.dirs) & (nortek.dirs < 361));
    Bindex2=find(nortek.dirs > 360);
    nortek.dirs(Bindex2)=nortek.dirs(Bindex2)-360;
    nortek.dirs(Aindex)=nortek.dirs(Aindex)+360;

    %plot the directional energy spectrum
    fig2=figure('Position',[scrsz]);
    subplot(1,2,1);
    h1 = plot(dirs(index2),EMEPdir(index2),'b');
    hold on
    h1a= plot(dirs(index1),EMEPdir(index1),'b');
    h2 = plot(dirs(index2),IMLMdir(index2),'r');
    h2a= plot(dirs(index1),IMLMdir(index1),'r');
    %plot the nortek data
    h3 = plot(nortek.dirs(Bindex2),nortekdir(Bindex2),'k');
    h3a = plot(nortek.dirs(Bindex),nortekdir(Bindex),'k');
    h3b= plot(nortek.dirs(Aindex),nortekdir(Aindex),'k');
else
    %for diwasp spectra do nothing
    
    %for nortek
    Aindex=find(nortek.dirs > 180);
    Bindex=find(nortek.dirs < 181);
    nortek.dirs(Aindex)=nortek.dirs(Aindex)-360;
    %plot the directional energy spectrum
    fig2=figure('Position',[scrsz]);
    subplot(1,2,1);
    h1 = plot(dirs,EMEPdir,'b');
    hold on
    h2 = plot(dirs,IMLMdir,'r');
    %plot the nortek data
    h3a = plot(nortek.dirs(Bindex),nortekdir(Bindex),'k');
    h3b = plot(nortek.dirs(Aindex),nortekdir(Aindex),'k');
    
end

legend([h1, h2, h3a],'EMEP','IMLM','nortek quickwave','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,EMEPfreq,'b');
hold on
plot(freqs,IMLMfreq,'r');
plot(nortekfreqs,nortekfreq,'k');
legend('EMEP','IMLM','nortek quickwave','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 EMEP AST and radial velocities

%calculate the 0,1,2 moments
m0=sum(EMEPfreq*freqres);                              
m1=sum(freqs.*EMEPfreq*freqres);
m2=sum((freqs.^2).*EMEPfreq*freqres);
% Calculate the Sig wave height
EMEP_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(EMEPfreq);
EMEP_Tp=1/(freqs(I));

%Calculate the Direction of the peak period DTp
[P,I]=max(real(EMEP.S(I,:)));
EMEP_DTp=dirs(I);

%Calculate the Dominant Direction Dp
[P,I]=max(EMEPdir);
EMEP_Dp=dirs(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['EMEP']);
disp(['SigH (meters): ' num2str(EMEP_Hsig)]);
disp(['peak period (seconds): ' num2str(EMEP_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(EMEP_DTp, EMEP.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(EMEP_Dp, EMEP.xaxisdir))]);
disp(['  ']);

E_WI.hsig=EMEP_Hsig;
E_WI.tp=EMEP_Tp;
E_WI.dtp=compangle(EMEP_DTp, EMEP.xaxisdir);
E_WI.dp=compangle(EMEP_Dp, EMEP.xaxisdir);


% For IMLM AST and radial velocities

%calculate the 0,1,2 moments
m0=sum(IMLMfreq*freqres);                              
m1=sum(freqs.*IMLMfreq*freqres);
m2=sum((freqs.^2).*IMLMfreq*freqres);
% Calculate the Sig wave height
IMLM_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(IMLMfreq);
IMLM_Tp=1/(freqs(I));

%Calculate the Direction of the peak period DTp
[P,I]=max(real(IMLM.S(I,:)));
IMLM_DTp=dirs(I);

%Calculate the Dominant Direction Dp
[P,I]=max(IMLMdir);
IMLM_Dp=dirs(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['IMLM']);
disp(['SigH (meters): ' num2str(IMLM_Hsig)]);
disp(['peak period (seconds): ' num2str(IMLM_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(IMLM_DTp, IMLM.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(IMLM_Dp, IMLM.xaxisdir))]);
disp(['  ']);

I_WI.hsig=IMLM_Hsig;
I_WI.tp=IMLM_Tp;
I_WI.dtp=compangle(IMLM_DTp, IMLM.xaxisdir);
I_WI.dp=compangle(IMLM_Dp, IMLM.xaxisdir);



%for nortek generated spectrum

%calculate the 0,1,2 moments
m0=sum(nortekfreq*freqres);                              
m1=sum(nortekfreqs.*nortekfreq*freqres);
m2=sum((nortekfreqs.^2).*nortekfreq*freqres);
% Calculate the Sig wave height
nortek_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(nortekfreq);
nortek_Tp=1/(nortekfreqs(I));

%Calculate the Direction of the peak period DTp
[P,I]=max(real(nortek.S(I,:)));
nortek_DTp=nortekdirs(I);

%Calculate the Dominant Direction Dp
[P,I]=max(nortekdir);
nortek_Dp=nortekdirs(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['nortek']);
disp(['SigH (meters): ' num2str(nortek_Hsig)]);
disp(['peak period (seconds): ' num2str(nortek_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(nortek_DTp, nortek.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(nortek_Dp, nortek.xaxisdir))]);
disp(['  ']);

nortek_WI.hsig=nortek_Hsig;
nortek_WI.tp=nortek_Tp;
nortek_WI.dtp=compangle(nortek_DTp, nortek.xaxisdir);
nortek_WI.dp=compangle(nortek_Dp, nortek.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);