function [E_UVW_F_WI,E_UVW_D_WI,I_UVW_F_WI,I_UVW_D_WI,fig3,fig4]=waveplot2(EMEP,IMLM,freqE,freqI,dirE,dirI,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')

%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;
freqresH=0.005;
freqs=[0.01:0.01:.4];
freqsH=[0.005:0.005:0.4];
dirres=2;
dirresH=1;
dirs=[-180:2:180];
dirsH=[-180:1:180];

%set up the directions
adj_angle=90-sangle+360+180;

% plot the spectrum generated through DIWASP 
scrsz = get(0,'ScreenSize');
fig3=figure('Position',[scrsz]);
subplot(2,3,1);
subplotspec(EMEP,4);
title('EMEP UVW');

subplot(2,3,2);
subplotspec(freqE,4);
title('EMEP UVW, freq = 0.005');

subplot(2,3,3);
subplotspec(dirE,4);
title('EMEP UVW, dir res = 1');

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

subplot(2,3,5);
subplotspec(freqI,4);
title('IMLM UVW, freq = 0.005');

subplot(2,3,6);
subplotspec(dirI,4);
title('IMLM UVW, dir res = 1');

%calculate just the dir energy spectrum on a single graph

EMEPdir=sum(EMEP.S)*freqres;
IMLMdir=sum(real(IMLM.S))*freqres;
freqEdir=sum(freqE.S)*freqresH;
dirEdir=sum(dirE.S)*freqres;
freqIdir=sum(real(freqI.S))*freqresH;
dirIdir=sum(real(dirI.S))*freqres;

%calculate just the frequency energy spectrum
EMEPfreq=sum(EMEP.S')*dirres;
IMLMfreq=sum(real(IMLM.S)')*dirres;
freqEfreq=sum(freqE.S')*dirres;
dirEfreq=sum(dirE.S')*dirresH;
freqIfreq=sum(real(freqI.S)')*dirres;
dirIfreq=sum(real(dirI.S)')*dirresH;

%plot the directional energy spectrum
fig4=figure('Position',[scrsz]);
subplot(1,2,1);
plot(dirs,EMEPdir,'b');
hold on
plot(dirs,freqEdir,'r');
plot(dirsH,dirEdir,'g');
axis(axis);
plot(dirs,IMLMdir,'c');
plot(dirs,freqIdir,'m');
plot(dirsH,dirIdir,'k');


legend('EMEP uvw','EMEP uvw freq=0.005','EMEP uvw dir=1','IMLM uvw','IMLM uvw freq=0.005','IMLM uvw dir=1','location','best');
title('directional wave spectrum integrated over frequency');
xlabel('axis angle (degrees true)');
ylabel('m^2 / hz');

%plot the frequency energy spectrum
subplot(1,2,2);
plot(freqs,EMEPfreq,'b');
hold on
plot(freqsH,freqEfreq,'r');
plot(freqs,dirEfreq,'g');
plot(freqs,IMLMfreq,'c');
plot(freqsH,freqIfreq,'m');
plot(freqs,dirIfreq,'k');
legend('EMEP uvw','EMEP uvw freq=0.005','EMEP uvw dir=1','IMLM uvw','IMLM uvw freq=0.005','IMLM uvw dir=1','location','best');
title('directional wave spectrum integrated over direction');
xlabel('frequency in Hz');
ylabel('m^2 / deg');

% ______Calculate and display the wave parameters SigH, Tp, Dp,DTp_______

%For EMEP uvw

%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
Hsig=4*sqrt(m0);

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

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

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

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

% For IMLM uvw

%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
Hsig=4*sqrt(m0);

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

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

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

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

% for EMEP freq=0.005

%calculate the 0,1,2 moments
m0=sum(freqEfreq*freqresH);                              
m1=sum(freqsH.*freqEfreq*freqresH);
m2=sum((freqsH.^2).*freqEfreq*freqresH);
% Calculate the Sig wave height
E_UVW_F_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(freqEfreq);
E_UVW_F_Tp=1/(freqsH(I));

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['EMEP UVW freq=0.005']);
disp(['SigH (meters): ' num2str(E_UVW_F_Hsig)]);
disp(['peak period (seconds): ' num2str(E_UVW_F_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(E_UVW_F_DTp, freqE.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(E_UVW_F_Dp, freqE.xaxisdir))]);
disp(['  ']);

E_UVW_F_WI.hsig=E_UVW_F_Hsig;
E_UVW_F_WI.tp=E_UVW_F_Tp;
E_UVW_F_WI.dtp=compangle(E_UVW_F_DTp, freqE.xaxisdir);
E_UVW_F_WI.dp=compangle(E_UVW_F_Dp, freqE.xaxisdir);

%for EMEP uvw dir=1

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

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

%Calculate the Direction of the peak period DTp
[P,I]=max(real(dirE.S(I,:)));
E_UVW_D_DTp=dirsH(I);

%Calculate the Dominant Direction Dp
[P,I]=max(dirEdir);
E_UVW_D_Dp=dirsH(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['EMEP UVW dir=1']);
disp(['SigH (meters): ' num2str(E_UVW_D_Hsig)]);
disp(['peak period (seconds): ' num2str(E_UVW_D_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(E_UVW_D_DTp, dirE.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(E_UVW_D_Dp, dirE.xaxisdir))]);
disp(['  ']);

E_UVW_D_WI.hsig=E_UVW_D_Hsig;
E_UVW_D_WI.tp=E_UVW_D_Tp;
E_UVW_D_WI.dtp=compangle(E_UVW_D_DTp, dirE.xaxisdir);
E_UVW_D_WI.dp=compangle(E_UVW_D_Dp, dirE.xaxisdir);

% for IMLM freq=0.005

%calculate the 0,1,2 moments
m0=sum(freqIfreq*freqresH);                              
m1=sum(freqsH.*freqIfreq*freqresH);
m2=sum((freqsH.^2).*freqIfreq*freqresH);
% Calculate the Sig wave height
I_UVW_F_Hsig=4*sqrt(m0);

% Calculate the peak period Tp
[P,I]=max(freqIfreq);
I_UVW_F_Tp=1/(freqsH(I));

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

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

%Display on the screen the SigH,Tp,Dp,DTp
disp(['IMLM UVW freq=0.005']);
disp(['SigH (meters): ' num2str(I_UVW_F_Hsig)]);
disp(['peak period (seconds): ' num2str(I_UVW_F_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(I_UVW_F_DTp, freqI.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(I_UVW_F_Dp, freqI.xaxisdir))]);
disp(['  ']);

I_UVW_F_WI.hsig=I_UVW_F_Hsig;
I_UVW_F_WI.tp=I_UVW_F_Tp;
I_UVW_F_WI.dtp=compangle(I_UVW_F_DTp, freqI.xaxisdir);
I_UVW_F_WI.dp=compangle(I_UVW_F_Dp, freqI.xaxisdir);

%for IMLM uvw dir=1

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

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

%Calculate the Direction of the peak period DTp
[P,I]=max(real(dirI.S(I,:)));
I_UVW_D_DTp=dirsH(I);

%Calculate the Dominant Direction Dp
[P,I]=max(dirIdir);
I_UVW_D_Dp=dirsH(I);

%Display on the screen the SigH,Tp,Dp,DTp
disp(['IMLM UVW dir=1']);
disp(['SigH (meters): ' num2str(I_UVW_D_Hsig)]);
disp(['peak period (seconds): ' num2str(I_UVW_D_Tp)]);
disp(['Dir of peak period: ' num2str(compangle(I_UVW_D_DTp, dirI.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(I_UVW_D_Dp, dirI.xaxisdir))]);
disp(['  ']);

I_UVW_D_WI.hsig=I_UVW_D_Hsig;
I_UVW_D_WI.tp=I_UVW_D_Tp;
I_UVW_D_WI.dtp=compangle(I_UVW_D_DTp, freqI.xaxisdir);
I_UVW_D_WI.dp=compangle(I_UVW_D_Dp, freqI.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);