function [E_radial_WI,fig1,fig2]=waveplot(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')


%what is the frequency and dir resolution for those generated in DIWASP
freqres=0.01;
freqs=[0.01:0.01:.4];
dirres=2;
dirs=[0:2:360];


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

%__________________________________________________________________________
%Use this if you want to calculate confidence bounds for the frequency spec

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


%calculate the conf limits throughout the frequency spectrum
EMEPradialfreqUP=EMEPradialfreq*coeffUp;
EMEPradialfreqLOW=EMEPradialfreq*coeffLow;

%__________________________________________________________________________

%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
    %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);
h1=plot(freqs,EMEPradialfreq,'b');
hold on
h2=plot(freqs,EMEPradialfreqUP,'b--');
h2a=plot(freqs,EMEPradialfreqLOW,'b--');
legend([h1,h2],'EMEP radial','95% confidence limits','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);