% plot_inwave_dir_test
%
% 1000 x 1200 with waves from left
%
% jcw 9/9/2021
%

cd E:\data\models\InWave\advect_test\output_poseidon

colormap('hsv')
qq=hsv;
qq=flipud(qq);

% so we tested all at 5 dx 
% 2 5 10 20 NoVary : dtheta
% s10 and s64 spread
% dhilbert and ihilbert 
% =20 runs

%here are the 5 dx grids for advection of wave from left
%
%decide trhat the 10dx and 5 dx are about the same, but just keep with 5 dx
%x_rho10=ncread('ocean_inwave_test_dir_10dx10dtheta.nc','x_rho');
%y_rho10=ncread('ocean_inwave_test_dir_10dx10dtheta.nc','y_rho');
%
ot=ncread('ocean_inwave_test_dir_5dx20dthetas64i.nc','ocean_time');
x_rho5=ncread('ocean_inwave_test_dir_5dx20dthetas64i.nc','x_rho');
y_rho5=ncread('ocean_inwave_test_dir_5dx20dthetas64i.nc','y_rho');

Hwave5dx20dthetas64i=squeeze(ncread('ocean_inwave_test_dir_5dx20dthetas64i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas64i=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas64i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx5dthetas64i=squeeze(ncread('ocean_inwave_test_dir_5dx5dthetas64i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx2dthetas64i=squeeze(ncread('ocean_inwave_test_dir_5dx2dthetas64i.nc','Hwave',[1 1 1600],[Inf Inf 1]));

Hwave5dx20dthetas64d=squeeze(ncread('ocean_inwave_test_dir_5dx20dthetas64d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas64d=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas64d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx5dthetas64d=squeeze(ncread('ocean_inwave_test_dir_5dx5dthetas64d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx2dthetas64d=squeeze(ncread('ocean_inwave_test_dir_5dx2dthetas64d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
%
Hwave5dx20dthetas10i=squeeze(ncread('ocean_inwave_test_dir_5dx20dthetas10i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas10i=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas10i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx5dthetas10i=squeeze(ncread('ocean_inwave_test_dir_5dx5dthetas10i.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx2dthetas10i=squeeze(ncread('ocean_inwave_test_dir_5dx2dthetas10i.nc','Hwave',[1 1 1600],[Inf Inf 1]));

Hwave5dx20dthetas10d=squeeze(ncread('ocean_inwave_test_dir_5dx20dthetas10d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas10d=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas10d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx5dthetas10d=squeeze(ncread('ocean_inwave_test_dir_5dx5dthetas10d.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx2dthetas10d=squeeze(ncread('ocean_inwave_test_dir_5dx2dthetas10d.nc','Hwave',[1 1 1600],[Inf Inf 1]));

Hwave5dx10dthetas64i_novary=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas64i_novary.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas64d_novary=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas64d_novary.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas10i_novary=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas10i_novary.nc','Hwave',[1 1 1600],[Inf Inf 1]));
Hwave5dx10dthetas10d_novary=squeeze(ncread('ocean_inwave_test_dir_5dx10dthetas10d_novary.nc','Hwave',[1 1 1600],[Inf Inf 1]));

%%%%%

% hwave for dhilbert
scale=1; %/sqrt(2);
S64=figure;
colormap(qq)
for mm=1:1 %1:2:size(ot)
tidx=1%mm;
subplot(621)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas64d(:,:,tidx))*scale);colorbar
% pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(623)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas64d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(625)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas64d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(627)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas64d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(629)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas64d_novary(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
%subplot(6,1,11)
% pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx45dthetas64d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
%  caxis([0 3])
%  set(gca,'xlim',[0 1000])
%drawnow
%pause(0.2)

%%%%%%%%%%%%%%%%%%%%%%%

subplot(622)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas64i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(624)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas64i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(626)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas64i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(628)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas64i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(6,2,10)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas64i_novary(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])

%subplot(621);   title('2 dtheta s64 d')    %s64 Spectra, d/theta = 2, direct construct
%subplot(623);   title('5 dtheta s64 d')
%subplot(625);   title('10 dtheta s64 d')
%subplot(627);   title('20 dtheta s64 d')
%subplot(629);   title('novary  s64 d')
%subplot(622);   title('2 dtheta s64 i')
%subplot(624);   title('5 dtheta s64 i')
%subplot(626);   title('10 dtheta s64 i')
%subplot(628);   title('20 dtheta s64 i')
%subplot(6,2,10);title('novary  s64 i')

subplot(621);  text(1200, 1600, 'S64 Spectra')
subplot(621);   title('d\theta = 2,  direct construct'); ylabel('Distance, m')
subplot(623);   title('d\theta = 5,  direct construct')
subplot(625);   title('d\theta = 10, direct construct')
subplot(627);   title('d\theta = 20, direct construct')
subplot(629);   title('novary      , direct construct'); xlabel('Distance, m')
subplot(622);   title('d\theta = 2,  integrated construct')
subplot(624);   title('d\theta = 5,  integrated construct')
subplot(626);   title('d\theta = 10, integrated construct')
subplot(628);   title('d\theta = 20, integrated construct')
subplot(6,2,10);title('novary      , integrated construct')


end

scale=1; %/sqrt(2);
S10=figure;
colormap(qq)
for mm=1:1 %1:2:size(ot)
tidx=1%mm;
subplot(621)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas10d(:,:,tidx))*scale);colorbar
% pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(623)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas10d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(625)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(627)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(629)
 pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10d_novary(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
%subplot(6,1,11)
% pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx45dthetas64d(:,:,tidx))*scale);colorbar
%  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10d(:,:,tidx))*scale);colorbar
%  caxis([0 3])
%  set(gca,'xlim',[0 1000])
drawnow
pause(0.2)

%%%%%%%%%%%%%%%%%%%%%%%

subplot(622)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx2dthetas10i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(624)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx5dthetas10i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(626)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(628)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx20dthetas10i(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])
subplot(6,2,10)
  pcolorjw(x_rho5,y_rho5,squeeze(Hwave5dx10dthetas10i_novary(:,:,tidx))*scale);colorbar
  caxis([0 3])
  set(gca,'xlim',[0 1000])

%subplot(621);   title('2 dtheta s10 d')
%subplot(623);   title('5 dtheta s10 d')
%subplot(625);   title('10 dtheta s10 d')
%subplot(627);   title('20 dtheta s10 d')
%subplot(629);   title('novary  s10 d')
%subplot(622);   title('2 dtheta s10 i')
%subplot(624);   title('5 dtheta s10 i')
%subplot(626);   title('10 dtheta s10 i')
%subplot(628);   title('20 dtheta s10 i')
%subplot(6,2,10);title('novary  s10 i')

subplot(621);  text(1200, 1600, 'S10 Spectra')
subplot(621);   title('d\theta = 2,  direct construct'); ylabel('Distance, m')
subplot(623);   title('d\theta = 5,  direct construct')
subplot(625);   title('d\theta = 10, direct construct')
subplot(627);   title('d\theta = 20, direct construct')
subplot(629);   title('novary      , direct construct'); xlabel('Distance, m')
subplot(622);   title('d\theta = 2,  integrated construct')
subplot(624);   title('d\theta = 5,  integrated construct')
subplot(626);   title('d\theta = 10, integrated construct')
subplot(628);   title('d\theta = 20, integrated construct')
subplot(6,2,10);title('novary      , integrated construct')



end



% down here lets comute teh longshore averaged cross-shore means and stdevs
% use stations files at 1 sec dt
% decide to use his files, save every second.
%clear
%
for mm=1:20
  mm
  if mm==1; netcdf_load('ocean_inwave_test_dir_5dx2dthetas64d.nc'); end
  if mm==2; netcdf_load('ocean_inwave_test_dir_5dx5dthetas64d.nc'); end
  if mm==3; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64d.nc'); end
  if mm==4; netcdf_load('ocean_inwave_test_dir_5dx20dthetas64d.nc'); end
  if mm==5; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64d_novary.nc'); end
%
  if mm==6; netcdf_load('ocean_inwave_test_dir_5dx2dthetas64i.nc'); end
  if mm==7; netcdf_load('ocean_inwave_test_dir_5dx5dthetas64i.nc'); end
  if mm==8; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64i.nc'); end
  if mm==9; netcdf_load('ocean_inwave_test_dir_5dx20dthetas64i.nc'); end
  if mm==10; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64i_novary.nc'); end
%
  if mm==11; netcdf_load('ocean_inwave_test_dir_5dx2dthetas10d.nc'); end
  if mm==12; netcdf_load('ocean_inwave_test_dir_5dx5dthetas10d.nc'); end
  if mm==13; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10d.nc'); end
  if mm==14; netcdf_load('ocean_inwave_test_dir_5dx20dthetas10d.nc'); end
  if mm==15; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10d_novary.nc'); end
%
  if mm==16; netcdf_load('ocean_inwave_test_dir_5dx2dthetas10i.nc'); end
  if mm==17; netcdf_load('ocean_inwave_test_dir_5dx5dthetas10i.nc'); end
  if mm==18; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10i.nc'); end
  if mm==19; netcdf_load('ocean_inwave_test_dir_5dx20dthetas10i.nc'); end
  if mm==20; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10i_novary.nc'); end
%figure
%plot(squeeze(Hwave(end-1,:)))
%for mm=1100:tidxe
%  pcolorjw(Hwave(:,:,mm)')
%  caxis([1 3])
%  colorbar
%  drawnow
%  pause(0.2)
%end
  tidxs=500; tidxe=length(ocean_time);
  clear Hwavemean2 Hwavevar2
  for i=1:size(Hwave,1)
    for j=1:size(Hwave,2)
      Hwavemean2(i,j)=mean(Hwave(i,j,tidxs:tidxe).^2);
      Hwavevar2(i,j)=std(Hwave(i,j,tidxs:tidxe).^2);
    end
  end
%figure; pcolorjw(Hwavemean2'); colorbar
%figure; pcolorjw(Hwavevar2'); colorbar
%now avg in along shore
GF=Hwavevar2./Hwavemean2;
GF=mean(GF,2);
%Hwavemean2=mean(Hwavemean2,2);
%Hwavevar2=mean(Hwavevar2,2);
%
  if ismember(mm,[1 6 11 16]); GF2=GF;; end
  if ismember(mm,[2 7 12 17]); GF5=GF;; end
  if ismember(mm,[3 8 13 18]); GF10=GF;; end
  if ismember(mm,[4 9 14 19]); GF20=GF;; end
  if ismember(mm,[5 10 15 20]); GFNV=GF;; end

%figure; colormap('jet')
  if ismember(mm, [5 15])
    if ismember(mm, [5]);figure(S64);end
    if ismember(mm, [15]);figure(S10);end
    subplot(6,2,11)
    hold on
    plot(x_rho(:,1),GF2'/(GF2(1)),'g')
    plot(x_rho(:,1),GF5'/(GF5(1)),'k')
    plot(x_rho(:,1),GF10'/(GF10(1)),'r')
    plot(x_rho(:,1),GF20'/(GF20(1)))
    plot(x_rho(:,1),GFNV'/(GFNV(1)))
    xlabel('dist (m)')
    ylabel('Relative GF')
    %title('s10 d')
%   title('s64 i')
    legend('2', '5', '10', '20','NV')
    set(gca,'ylim',[0 1])
    set(gca,'xlim',[0 1000])
  end
  if ismember(mm, [10 20])
    if ismember(mm, [10]);figure(S64);end
    if ismember(mm, [20]);figure(S10);end
    subplot(6,2,12)
    hold on
    plot(x_rho(:,1),GF2'/(GF2(1)),'g')
    plot(x_rho(:,1),GF5'/(GF5(1)),'k')
    plot(x_rho(:,1),GF10'/(GF10(1)),'r')
    plot(x_rho(:,1),GF20'/(GF20(1)))
    plot(x_rho(:,1),GFNV'/(GFNV(1)))
    xlabel('dist (m)')
    ylabel('Relative GF')
    %title('s10 d')
%   title('s64 i')
    legend('2', '5', '10', '20','NV')
    set(gca,'ylim',[0 1])
    set(gca,'xlim',[0 1000])
  end
%  set(gcf,'pos',[988   139   760   878])
%set(gcf, 'DefaultFigureRenderer', 'painters');
%set(gcf, 'DefaultFigureRenderer', 'OpenGL');
%set(gcf, 'DefaultFigureRenderer', 'zbuffer');
%  print -dpng figX_s10_groupiness.png
%  print -dpng figX_s64_groupiness.png
%  print -depsc figX_s10_groupiness.eps
%  print -depsc figX_s64_groupiness.eps
%  new jcw zbuffer way::::
%  exportgraphics(gcf,'figX_s64_groupiness.eps', 'ContentType', 'vector')
%  exportgraphics(gcf,'figX_s10_groupiness.eps', 'ContentType', 'vector')
end
%   save inwave_dir_test_10may2022.mat

%%%% here do a quick fig of the spectra 10 and 64
dirs=[180:10:360];
ndirs=length(dirs);
freqs=[0.04443
0.0542
0.06396
0.07373
0.0835
0.09326
0.10303
0.11279
0.12256
0.13232
0.14209
0.15186
0.16162
0.17139
0.18115
0.19092
0.20068
0.21045
0.22021
0.22998
0.23975
0.24951
0.25928
0.26904
0.27881
0.28857
0.29834
0.30811
0.31787];
nfreqs=length(freqs);
Tp=10;
Hs=2;
Ohm=2*pi./freqs;
[ Sf, Amp, Phase ] = jonswap_lip( Ohm, Hs, Tp);
%make dir spreading
%https://www.sciencedirect.com/topics/engineering/directional-spreading
%  for each spread s=10 for 24 deg spread sea, and s=64 for 10deg spread swell 
s=64;  % 64;   %10;
S=zeros(nfreqs,ndirs);
for mm=1:ndirs
  S(:,mm)=Sf*(cosd(abs(dirs(mm)-dirs(12))/2)^(2*s));
end
% compute Hsig
Sdir = sum(S,2)*abs(dirs(2)-dirs(1))*pi/180;
m0=sum(Sdir.*gradient(freqs));
hm0gew=4.004*sqrt(m0);
scale=Hs^2/16/m0;
S64=S.*scale;
%
s=10;  % 64;   %10;
S=zeros(nfreqs,ndirs);
for mm=1:ndirs
  S(:,mm)=Sf*(cosd(abs(dirs(mm)-dirs(12))/2)^(2*s));
end
% compute Hsig
Sdir = sum(S,2)*abs(dirs(2)-dirs(1))*pi/180;
m0=sum(Sdir.*gradient(freqs));
hm0gew=4.004*sqrt(m0);
scale=Hs^2/16/m0;
S10=S.*scale;
%
figure
colormap('jet')
subplot(121)
pcolorjw(freqs,dirs,S10'); caxis([0 4]); colorbar
ylabel('Direction, deg')
xlabel('frequency, Hz')
title('Energy spectra (m^2/deg/Hz) for spread s = 10')
subplot(122)
pcolorjw(freqs,dirs,S64'); caxis([0 8]); colorbar
xlabel('frequency, Hz')
title(['Energy spectra (m^2/deg/Hz) for spread s = 64', ...
       '                                             '])
set(gcf,'pos',[771         537        1090         318])
% print -dpng figX_s10s64_Espectra.png
%
% new way here
figure
%colormap('jet')
%cmocean('turbid')
%cmocean('balance')
%cmocean('deep')
%cmocean('thermal')
subplot(121)
polarPcolor(freqs',dirs,S10,'Nspokes',13);
caxis([0 4]);
cmocean('tempo')
%ylabel('Direction, deg')
%xlabel('frequency, Hz')
title(['Directional spectra (m^2/deg/Hz) for spread s = 10';...
       '                                                  '])
subplot(122)
polarPcolor(freqs',dirs,S64,'Nspokes',13);
caxis([0 8]);
cmocean('tempo')
%xlabel('frequency, Hz')
title(['Directional spectra (m^2/deg/Hz) for spread s = 64';...
       '                                                  '])
set(gcf,'pos',[771         537        1090         318])
set(gcf, 'DefaultFigureRenderer', 'painters');
% print -dpng fig3_s10s64_DirSpec.png
% print -depsc figX_s10s64_DirSpec.eps

%directional spectra
 

%%%%%%%%%%%
% donw here we compute the Hsig low = IG wave energy
for mm=1:20
  mm
  if mm==1; netcdf_load('ocean_inwave_test_dir_5dx2dthetas64d.nc'); end
  if mm==2; netcdf_load('ocean_inwave_test_dir_5dx5dthetas64d.nc'); end
  if mm==3; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64d.nc'); end
  if mm==4; netcdf_load('ocean_inwave_test_dir_5dx20dthetas64d.nc'); end
  if mm==5; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64d_novary.nc'); end
%
  if mm==6; netcdf_load('ocean_inwave_test_dir_5dx2dthetas64i.nc'); end
  if mm==7; netcdf_load('ocean_inwave_test_dir_5dx5dthetas64i.nc'); end
  if mm==8; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64i.nc'); end
  if mm==9; netcdf_load('ocean_inwave_test_dir_5dx20dthetas64i.nc'); end
  if mm==10; netcdf_load('ocean_inwave_test_dir_5dx10dthetas64i_novary.nc'); end
%
  if mm==11; netcdf_load('ocean_inwave_test_dir_5dx2dthetas10d.nc'); end
  if mm==12; netcdf_load('ocean_inwave_test_dir_5dx5dthetas10d.nc'); end
  if mm==13; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10d.nc'); end
  if mm==14; netcdf_load('ocean_inwave_test_dir_5dx20dthetas10d.nc'); end
  if mm==15; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10d_novary.nc'); end
%
  if mm==16; netcdf_load('ocean_inwave_test_dir_5dx2dthetas10i.nc'); end
  if mm==17; netcdf_load('ocean_inwave_test_dir_5dx5dthetas10i.nc'); end
  if mm==18; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10i.nc'); end
  if mm==19; netcdf_load('ocean_inwave_test_dir_5dx20dthetas10i.nc'); end
  if mm==20; netcdf_load('ocean_inwave_test_dir_5dx10dthetas10i_novary.nc'); end

  tidxs=500; tidxe=length(ocean_time);
%  tidxs=200; tidxe=length(ocean_time);
  clear Hsig_mod_lo
% for i=100:size(zeta,1)  %saved as workspace.mat
  for i=1:size(zeta,1)  %saved as workspace_alli.mat
   disp(['mm is ', num2str(mm), ' for i of ', num2str(i), ' to i = ', num2str(size(zeta,1))])
    for j=1:size(zeta,2)
      clear Sf_lo Fn_lo
      zetas=squeeze(zeta(i,j,tidxs:tidxe));
      zvar = zetas- mean(zetas);
      zvar=detrend(zvar);
      Fn_lo = fft(zvar)/length(zvar);
%
      if (mm==1)
        dt = 1;
        nt = length(zetas); % length of time series
        df = 1/(nt*dt);
        for ii = 1:nt
          ff (ii) = ii*df;
        end
        T = 1./ff;
        nf1 = length (ff);
      end
%
  %%%% band passing higher frequencies
      for k = 1:nf1
          if ff(k) > 0.05 % I think this is the cutoff Ap used
              Fn_lo(k) = 0;
          end
          if ff(k) < 0.005 % I think this is the cutoff Ap used
              Fn_lo(k) = 0;
          end
      end
      for k = 1:nf1
          Sf_lo(k) =2*abs(Fn_lo(k))^2/df;
      end
  %jcw do an avg
%      for k = 8:nf1-7
%          Sf_lo(k) =sum(Sf_lo(k-7:k+7))/15;
%      end
  %interpolate to coarser
%     ff_new=[ff(1):0.0078:0.05];
%     Sf_lo_new=interp1(ff,Sf_lo,ff_new);
%     sigma2 = sum(Sf_lo_new(:)*0.0078);
      sigma2 = sum(Sf_lo(:)*df);
      Hsig_mod_lo(i,j) = 4*sqrt(sigma2);
    end
  end
  eval(['IG',num2str(mm),'=mean(Hsig_mod_lo(:));'])
end
%
figure
hold on
XX=180; %60
%subplot(221)
  plot([2 5 10 20 XX],[IG1 IG2 IG3 IG4 IG5],'k-+','linewidth',2)
%  title('s64 d IG')
  xlabel('direction bin size')
  ylabel('IG Energy')
%subplot(223)
  plot([2 5 10 20 XX],[IG6 IG7 IG8 IG9 IG10],'r-+','linewidth',2)
%  title('s64 i IG')
%  xlabel('direction bin size')
%  ylabel('IG Energy')
%subplot(222)
  plot([2 5 10 20 XX],[IG11 IG12 IG13 IG14 IG15],'b-+','linewidth',2)
%  title('s10 d IG')
%  xlabel('direction bin size')
%  ylabel('IG Energy')
%subplot(224)
  plot([2 5 10 20 XX],[IG16 IG17 IG18 IG19 IG20],'g-+','linewidth',2)
%  title('s10 i IG')
%  xlabel('direction bin size')
%  ylabel('IG Energy')
legend('s 64 d','s 64 i','s 10 d','s 10 i')
set(gcf,'pos',[771         537        1090         318])
% print -dpng figX_IG_vs_bindirs.png
% print -depsc figX_IG_vs_bindirs.eps
% save load hsig_lo.mat

