root/DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplot.m

function [E_UVW_WI,I_UVW_WI,E_range_WI,I_range_WI,wmon_WI,fig1,fig2]=waveplot(EMEP,IMLM,rangeE,rangeI,spec,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;

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

% plot the spectrum generated through DIWASP 

scrsz = get(0,'ScreenSize');
fig1=figure('Position',[scrsz]);
subplot(2,3,1);
subplotspec(EMEP,4);
title('EMEP uvw');

subplot(2,3,2);
subplotspec(IMLM,4);
title('IMLM uvw');

subplot(2,3,3);
subplotspec(wmon,4);
title('Wavesmon output');

subplot(2,3,4);
subplotspec(rangeE,4);
title('EMEP Range');

subplot(2,3,5);
subplotspec(rangeI,4);
title('IMLM Range');

%calculate just the dir energy spectrum on a single graph

EMEPdir=sum(EMEP.S)*freqres;
IMLMdir=sum(real(IMLM.S))*freqres;
EMEPrangedir=sum(rangeE.S)*freqres;
IMLMrangedir=sum(real(rangeI.S))*freqres;
wmondir=sum(wmon.S)*wmonfreqres;

%calculate just the frequency energy spectrum
EMEPfreq=sum(EMEP.S')*dirres;
IMLMfreq=sum(real(IMLM.S)')*dirres;
EMEPrangefreq=sum(rangeE.S')*dirres;
IMLMrangefreq=sum(real(rangeI.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=EMEP.degF;
chiUp=chi2inv(.975,degF);
chiLow=chi2inv(.025,degF);
coeffUp=degF/chiLow;
coeffLow=degF/chiUp;

%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 wavesmon
    Aindex=find(wmon.dirs < 0);               
    Bindex=find((-1 < wmon.dirs) & (wmon.dirs < 361));
    Bindex2=find(wmon.dirs > 360);
    wmon.dirs(Bindex2)=wmon.dirs(Bindex2)-360;
    wmon.dirs(Aindex)=wmon.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),EMEPrangedir(index2),'r');
    h2a= plot(dirs(index1),EMEPrangedir(index1),'r');
    h3 = plot(dirs(index2),IMLMrangedir(index2),'g');
    h3a= plot(dirs(index1),IMLMrangedir(index1),'g');
    %plot the wavesmon data
    h4 = plot(wmon.dirs(Bindex2),wmondir(Bindex2),'k');
    h4a = plot(wmon.dirs(Bindex),wmondir(Bindex),'k');
    h4b= plot(wmon.dirs(Aindex),wmondir(Aindex),'k');
    axis(axis);
    h5 = plot(dirs(index2),IMLMdir(index2),'c');
    h5a= plot(dirs(index1),IMLMdir(index1),'c');

else
    %for diwasp spectra do nothing
    
    %for wavesmon
    Aindex=find(wmon.dirs > 180);
    Bindex=find(wmon.dirs < 181);
    wmon.dirs(Aindex)=wmon.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,EMEPrangedir,'r');
    h3 = plot(dirs,IMLMrangedir,'g');
    %plot the wavesmon data
    h4 = plot(wmon.dirs(Aindex),wmondir(Aindex),'k');
    h4a = plot(wmon.dirs(Bindex),wmondir(Bindex),'k');
    axis(axis);
    h5 = plot(dirs,IMLMdir,'c');
end

legend([h1, h2, h3, h4, h5],'EMEP uvw','EMEP range','IMLM range','wavesmon','IMLM uvw','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,EMEPrangefreq,'r');
plot(freqs,IMLMrangefreq,'g');
plot(wmon.freqs,wmonfreq,'k');
axis(axis);
plot(freqs,IMLMfreq,'c');
legend('EMEP uvw','EMEP range','IMLM range','wavesmon','IMLM uvw','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 uvw
% Calculate the Sig wave height
EMEP_UVW_Hsig=Hsig(EMEP);
%Use the function HsigConf.m to calculate the sigH confidence limits
EMEP_HsConf=HsigConf(EMEP);

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

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['EMEP uvw']);
disp(['SigH (meters): ' num2str(EMEP_UVW_Hsig)]);
disp(['SigH 95% confidence limits: ' num2str(EMEP_HsConf)]);
disp(['peak period (seconds): ' num2str(EMEP_UVW_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(EMEP_UVW_DTp, EMEP.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(EMEP_UVW_Dp, EMEP.xaxisdir))]);
disp(['  ']);

E_UVW_WI.hsig=EMEP_UVW_Hsig;
E_UVW_WI.hconf=EMEP_HsConf;
E_UVW_WI.tp=EMEP_UVW_Tp;
E_UVW_WI.dtp=compangle(EMEP_UVW_DTp, EMEP.xaxisdir);
E_UVW_WI.dp=compangle(EMEP_UVW_Dp, EMEP.xaxisdir);


% For IMLM uvw

% Calculate the Sig wave height
IMLM_UVW_Hsig=Hsig(IMLM);
%Use the function HsigConf.m to calculate the sigH confidence limits
IMLM_HsConf=HsigConf(IMLM);

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

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['IMLM uvw']);
disp(['SigH (meters): ' num2str(IMLM_UVW_Hsig)]);
disp(['SigH 95% confidence limits: ' num2str(IMLM_HsConf)]);
disp(['peak period (seconds): ' num2str(IMLM_UVW_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(IMLM_UVW_DTp, IMLM.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(IMLM_UVW_Dp, IMLM.xaxisdir))]);
disp(['  ']);

I_UVW_WI.hsig=IMLM_UVW_Hsig;
I_UVW_WI.hconf=IMLM_HsConf;
I_UVW_WI.tp=IMLM_UVW_Tp;
I_UVW_WI.dtp=compangle(IMLM_UVW_DTp, IMLM.xaxisdir);
I_UVW_WI.dp=compangle(IMLM_UVW_Dp, IMLM.xaxisdir);

% for Range of EMEP

% Calculate the Sig wave height
EMEP_range_Hsig=Hsig(rangeE);
%Use the function HsigConf.m to calculate the sigH confidence limits
EMEP_range_HsConf=HsigConf(rangeE);

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

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['EMEP range']);
disp(['SigH (meters): ' num2str(EMEP_range_Hsig)]);
disp(['SigH 95% confidence limits: ' num2str(EMEP_range_HsConf)]);
disp(['peak period (seconds): ' num2str(EMEP_range_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(EMEP_range_DTp, rangeE.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(EMEP_range_Dp, rangeE.xaxisdir))]);
disp(['  ']);

E_range_WI.hsig=EMEP_range_Hsig;
E_range_WI.hconf=EMEP_range_HsConf;
E_range_WI.tp=EMEP_range_Tp;
E_range_WI.dtp=compangle(EMEP_range_DTp, rangeE.xaxisdir);
E_range_WI.dp=compangle(EMEP_range_Dp, rangeE.xaxisdir);


%for Range of IMLM

% Calculate the Sig wave height
IMLM_range_Hsig=Hsig(rangeI);
%Use the function HsigConf.m to calculate the sigH confidence limits
IMLM_range_HsConf=HsigConf(rangeI);

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

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['IMLM range']);
disp(['SigH (meters): ' num2str(IMLM_range_Hsig)]);
disp(['SigH 95% confidence limits: ' num2str(IMLM_range_HsConf)]);
disp(['peak period (seconds): ' num2str(IMLM_range_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(IMLM_range_DTp, rangeI.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(IMLM_range_Dp, rangeI.xaxisdir))]);
disp(['  ']);

I_range_WI.hsig=IMLM_range_Hsig;
I_range_WI.hconf=IMLM_range_HsConf;
I_range_WI.tp=IMLM_range_Tp;
I_range_WI.dtp=compangle(IMLM_range_DTp, rangeI.xaxisdir);
I_range_WI.dp=compangle(IMLM_range_Dp, rangeI.xaxisdir);

% for wavesmon

%calculate the 0,1,2 moments
m0=sum(wmonfreq*wmonfreqres);                              
m1=sum(wmon.freqs.*wmonfreq*wmonfreqres);
m2=sum((wmon.freqs.^2).*wmonfreq*wmonfreqres);
% Calculate the Sig wave height
wmon_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(wmonfreq);
wmon_Tp=1/(wmon.freqs(I));

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

%Calculate the Dominant Direction Dp
[P,I]=max(wmondir);
wmon_Dp=wmon.dirs(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['Wavesmon output']);
disp(['SigH (meters): ' num2str(wmon_Hsig)]);
disp(['peak period (seconds): ' num2str(wmon_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(wmon_DTp, wmon.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(wmon_Dp, wmon.xaxisdir))]);
disp(['  ']);

wmon_WI.hsig=wmon_Hsig;
wmon_WI.tp=wmon_Tp;
wmon_WI.dtp=compangle(wmon_DTp, wmon.xaxisdir);
wmon_WI.dp=compangle(wmon_Dp, wmon.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);