function waveplot3(EMEP,IMLM,rangeE,rangeI,spec) 

%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
sangle=281+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 

subplot(2,3,1);
subplotspec(EMEP.SMout,4);
title('EMEP Range');

subplot(2,3,2);
subplotspec(IMLM.SMout,4);
title('IMLM Range');

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

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

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

%calculate just the dir energy spectrum on a single graph

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

%calculate just the frequency energy spectrum
EMEPfreq=sum(EMEP.SMout.S')*dirres;
IMLMfreq=sum(real(IMLM.SMout.S)')*dirres;
EMEPrangefreq=sum(rangeE.SMout.S')*dirres;
IMLMrangefreq=sum(real(rangeI.SMout.S)')*dirres;
wmonfreq=sum(wmon.S')*wmondirres;

%plot the directional energy spectrum
figure;
subplot(1,2,1);
plot(dirs,EMEPdir,'b');
hold on
plot(dirs,EMEPrangedir,'r');
plot(dirs,IMLMrangedir,'g');
%need to make sure the wmon plot has the same x-axis as the others
Aindex=find(wmon.dirs > 180);
Bindex=find(wmon.dirs <181);
wmon.dirs(Aindex)=wmon.dirs(Aindex)-360;

plot(wmon.dirs(Aindex),wmondir(Aindex),'k');
plot(wmon.dirs(Bindex),wmondir(Bindex),'k');
axis(axis);
plot(dirs,IMLMdir,'c');
legend('EMEP range','EMEP range cpsd','IMLM range cpsd','wavesmon','','IMLM range');
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 range','EMEP range cpsd','IMLM range cpsd','wavesmon','IMLM range');
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 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.SMout.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 range']);
disp(['SigH (meters): ' num2str(Hsig)]);
disp(['peak period (seconds): ' num2str(Tp)]);
disp(['Dir of peak period: ' num2str(compangle(DTp, EMEP.SMout.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(Dp, EMEP.SMout.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.SMout.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 range']);
disp(['SigH (meters): ' num2str(Hsig)]);
disp(['peak period (seconds): ' num2str(Tp)]);
disp(['Dir of peak period: ' num2str(compangle(DTp, IMLM.SMout.xaxisdir))]);
disp(['Dominant Direction: ' num2str(compangle(Dp, IMLM.SMout.xaxisdir))]);
disp(['  ']);

% for Range of EMEP

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

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

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

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

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

%for Range of IMLM

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

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

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

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

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

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

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

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

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

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


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