% plot 2d spec for Delilah 13 Oct 1000

togo='test163';
eval([' cd E:\data\models\InWave\readswan\delilah08\',togo])

printfig=0 %0;  %1
%
netcdf_load('ocean_delilah_stations.nc')
psdfig=figure;
for mm=1:10
  clear Sf_lo Sf_hi Fn_lo Fn_hi ff
  if (mm<10);   subplot(3,3,mm); end
  hold on
  nn=mm;  %+9;
  if (mm==10); nn=19; end
 zetas=squeeze(zeta(nn,18001:36000));
  zvar = zetas- mean(zetas);
  zvar=detrend(zvar);
  Fn_lo = fft(zvar)/length(zvar);
  Fn_hi = Fn_lo;
%
  dt = 0.2; % 5 Hz
  nt = length(zetas); % length of time series
  df = 1/(nt*dt);
  for i = 1:nt
      ff (i) = i*df;
  end
  T = 1./ff;
  nf1 = length (ff);
%
  %%%% band passing higher frequencies
  for j = 1:nf1
      if ff(j) > 0.05 % I think this is the cutoff Ap used
          Fn_lo(j) = 0;
      end
      if ff(j) < 0.005 % I think this is the cutoff Ap used
          Fn_lo(j) = 0;
      end
     if ff(j) < 0.05
          Fn_hi(j) = 0;
      end
      if ff(j) > 1   % remove right half of the mirrorred spectra
          Fn_hi(j) = 0;
      end
  end
  for j = 1:nf1
      Sf_lo(j) =2*abs(Fn_lo(j))^2/df;
%     Sf_hi(j) =2*abs(Fn_hi(j))^2/df;
  end
  %jcw do an avg
  for j = 8:nf1-7
      Sf_lo(j) =sum(Sf_lo(j-7:j+7))/15;
 %    Sf_hi(j) =sum(Sf_hi(j-7:j+7))/15;
  end
 %interpolate to coarser
  ff_new=[ff(1):0.0078:0.05];
  Sf_lo_new=interp1(ff,Sf_lo,ff_new);
  if (mm<10); plot(ff_new,Sf_lo_new); end
  sigma2(mm) = sum(Sf_lo_new(:)*0.0078);
  Hsig_mod_lo(mm) = 4*sqrt(sigma2(mm));
  set(gca,'xlim',[0 0.05],'ylim',[0 0.8])
end

cd E:\data\models\InWave\readswan\Projects\Delilah\data\gages\rawTS
%load delilah_13octhr10.mat
load delilah_13octhr16.mat
paros=0; strain=1;
clear ff
if paros
  dt=1;
  nt = length(paros_time); % length of time series
elseif strain
  dt = 1/8; % 8 Hz
  nt = length(strain_time); % length of time series
% nt = length(strain_time)/4; % length of time series
end
df = 1/(nt*dt);
for i = 1:nt
    ff (i) = i*df;
end
T = 1./ff;
nf1 = length (ff);
figure(psdfig);
for mm=1:9
  clear Sf_lo Sf_hi Fn_lo Fn_hi
  subplot(3,3,mm)
  hold on
  if strain
   if (mm==1); zetas=squeeze(site10_strain); end
   if (mm==2); zetas=squeeze(site20_strain); end
   if (mm==3); zetas=squeeze(site30_strain); end
   if (mm==4); zetas=squeeze(site40_strain); end
   if (mm==5); zetas=squeeze(site50_strain); end
   if (mm==6); zetas=squeeze(site50_strain)*nan; end
   if (mm==7); zetas=squeeze(site70_strain); end
   if (mm==8); zetas=squeeze(site80_strain); end
   if (mm==9); zetas=squeeze(site90_strain); end
  elseif paros
   if (mm==1); zetas=squeeze(site10_paros); end
   if (mm==2); zetas=squeeze(site20_paros); end
   if (mm==3); zetas=squeeze(site30_paros); end
   if (mm==4); zetas=squeeze(site40_paros); end
   if (mm==5); zetas=squeeze(site50_paros); end
   if (mm==6); zetas=squeeze(site50_paros)*0; end
   if (mm==7); zetas=squeeze(site70_paros); end
   if (mm==8); zetas=squeeze(site50_paros)*0; end
   if (mm==9); zetas=squeeze(site50_paros)*0; end
  end
  zvar=detrend(zetas);
  zvar = zvar- mean(zvar);   %jcw just add
  Fn_lo = fft(zvar)/length(zvar);
  Fn_hi=Fn_lo;

% depth atten
%  hh(mm)=4.25;
  z=[0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ]; %  use mean of z!!!
% deep=[-0.28 -1.39 -0.87 -1.32 -1.88 -2.33 -3.03 -3.68 -4.25];  %water depth relative to MSL NGVD
  deep=[-0.57 -1.70 -1.92 -1.56 -1.65 -2.35 -2.97 -3.71 -4.24];  %water depth relative to MSL, based on grd4 per
% gage=[ 0.07 -0.77 -0.66 -0.82 -0.98 -1.02 -1.05 -1.66 -1.68];   %depth of sensor below MSL- this was current meters.
  gage=[-0.65 -0.86 -0.75 -0.86 -0.98 -1.22 -1.20 -3.04 -3.06];   %depth of sensor below MSL- strain gages
% did this after test 91
 
if(1)
  for j = 1:nf1
    D=z-deep;
    k (mm,j) = qkhfs(2.*3.14/T(j),D(mm))/D(mm); 
    atten (mm,j) = cosh (k(mm,j)*D(mm)) / cosh(k(mm,j)*(-deep(mm) + gage(mm)));
    Fn_lo(j) = Fn_lo(j)*atten(mm,j);
    Fn_hi(j) = Fn_hi(j)*atten(mm,j);
  end
end

  %%%% band passing higher frequencies
  for j = 1:nf1
      if ff(j) > 0.05 % I think this is the cutoff Ap used
          Fn_lo(j) = 0;
      end
      if ff(j) < 0.005 % I think this is the cutoff Ap used
          Fn_lo(j) = 0;
      end
      if ff(j) < 0.05
          Fn_hi(j) = 0;
      end
      if ff(j) > 0.4   % remove right half of the mirrorred spectra
          Fn_hi(j) = 0;
      end
  end
  for j = 1:nf1
      Sf_lo(j) =2*abs(Fn_lo(j))^2/df;
      Sf_hi(j) =2*abs(Fn_hi(j))^2/df;
  end
  for j = 16:nf1-15
      Sf_lo(j) =sum(Sf_lo(j-15:j+15))/31;
      Sf_hi(j) =sum(Sf_hi(j-15:j+15))/31;
  end

  ff_new=[ff(1):0.0078:ff(end)];
  Sf_lo_new=interp1(ff,Sf_lo,ff_new);
  Sf_hi_new=interp1(ff,Sf_hi,ff_new);
  sigma2(mm) = sum(Sf_lo_new(:)*0.0078);
  Hsig_obs_lo(mm) = 4*sqrt(sigma2(mm));
  sigma2(mm) = sum(Sf_hi_new(:)*0.0078);
  Hsig_obs_hi(mm) = 4*sqrt(sigma2(mm));

%
% plot(ff,Sf_lo,'r','linewidth',2)
  plot(ff_new,Sf_lo_new,'r','linewidth',2)
  set(gca,'xlim',[0 0.05],'ylim',[0 0.8])
  text(0.03,0.6,['C',num2str(mm),'0'])
  xlabel('f, Hz')
  ylabel('Energy spectra, m^2/Hz')
end
subplot(3,3,2)
if (~printfig)
 title(['Blue = InWave, Black = Observed  at hr 1600 on Oct 13, 1990 ', togo])
end
legend('Model', 'Observations')
set(gcf,'pos',[1051         118         787         578]);
%print -dpng delilah_13Oct1600_fft.png
%print -depsc E:\data2\papers\inwave\figures\delilah_13Oct1600_fft.eps


%now do X-shore Hsig profile
eval(['cd E:\data\models\InWave\readswan\delilah08\',togo])   % this is the one to use normlally

sites=['CM10 CM20 CM30 CM40 CM50 CM60 CM70 CM80 CM90'];
xloc=[125.06 144.99 169.97 188.94 207.41 229.25 245.00 295.21 370.12];
yloc=[985.95 985.94 985.61 985.95 985.88 986.08 985.91 985.97 986.11];
%%yloc=[1200.0 1200.0 1200.0 1200.0 1200.0 1200.0 1200.0 1200.0 1200.0];
zloc=[-0.28 -1.39 -0.87 -1.32 -1.88 -2.33 -3.03 -3.68 -4.25];
x_rho=ncread('ocean_delilah_his_00001.nc','x_rho');
y_rho=ncread('ocean_delilah_his_00001.nc','y_rho');
for mm=1:9
  xlocr(mm)=floor(interp1(x_rho(:,1),[1:1:size(x_rho,1)],xloc(mm)));
  ylocr(mm)=floor(interp1(y_rho(1,:),[1:1:size(y_rho,2)],yloc(mm)));
end
fname='ocean_delilah_avg.nc';
ot=ncread(fname,'ocean_time')/3600/24;
%for del08, this was 4 min avg, so use 99:105
Hwave2=squeeze(ncread(fname,'Hwave2',[1 ylocr(1) 1],[Inf 1 Inf]));
Hwave2=Hwave2(:,15:30);
% hwave_high  cross shore profile
figure
hold on
scale=1; %sqrt(2); %Hmo    use 1/sqrt(2) for Hrms
plot(xloc,Hsig_obs_lo,'ro')  % this is from above computation of strain gage data
%plot([xloc 950],[Hsig_obs_hi 1.823],'go')  % this is from above computation of strain gage data
plot([xloc 950],[Hsig_obs_hi 2.1],'go')  % this is from above computation of strain gage data
plot([xloc 950],Hsig_mod_lo*scale,'k.')  % this
plot(x_rho(:,ylocr(1)),squeeze(sqrt(mean(Hwave2'))),'k','linewidth',2)
%add grey lines to show full variations
plot(x_rho(:,ylocr(1)),squeeze(sqrt((Hwave2'))),'k--','linewidth',0.5)
%title(['Hmo\_high = sqrt(<Hwave^2>),  Dots are Observations during the hour, lines are model 10 min avgs'; ...
%       '                    1990 Oct 13, 1600                                                           '])
if (~printfig)
  title(['Hmo\_high and Hmo\_lo 1990 Oct 13, 1600 ', togo])
else
%  title(['Hmo\_high and Hmo\_lo on Oct 13, 1990 hr 1600 (no Doppler) '])
  if (strcmp(togo,'test62'))
    title(['Hmo\_high and Hmo\_lo on Oct 13, 1990 hr 1600 (no Doppler) '])
  end
% if (strcmp(togo,'test67'))
  if (strcmp(togo,'test51'))
    title(['Hmo\_high and Hmo\_lo on Oct 13, 1990 hr 1600 (yes Doppler) '])
  end
end
xlabel('Cross-shore distance, m ')
ylabel('Hmo, m')
legend('Hmo\_lo obs','Hmo\_hi obs','Hmo\_lo model','Hmo\_hi model')
text(xloc(1)-20,2.25,'C10')
text(xloc(2)-20,2.4,'C20')
text(xloc(3)-20,2.25,'C30')
text(xloc(4)-20,2.4,'C40')
text(xloc(5)-20,2.25,'C50')
text(xloc(6)-20,2.4,'C60')
text(xloc(7)-20,2.25,'C70')
text(xloc(8)-20,2.4,'C80')
text(xloc(9)-20,2.25,'C90')
text(940-60,2.25,'8m array')
%plot([xloc 950],[Hsig_obs_hi 1.823],'go')  % this is from above computation of strain gage data
set(gca,'xlim', [50 950])
grid on
set(gcf,'pos',[777         733        1007         313]);

%print -dpng delilah_grid.png


%%   add skill + bias
yi=Hsig_obs_lo;
xi=Hsig_mod_lo(1:end-1);
yi=[yi(1:5) yi(7:end)];
xi=[xi(1:5) xi(7:end)];

bias = sum(yi-xi)/length(yi)           =0.0850734
skill = 1- sum((yi-xi).^2)/sum((xi).^2)   =  0.9159114




