%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WAVES 1 - LAB 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PURPOSE : compute the frequency surface elevation spectrum from pressure time series.
% HISTORY : - origination in matlab
% - 01.2019 comments and rewriting, C. Peureux, cpeureux@ifremer.fr
% REQUIREMENTS : deep_waimea.mat

clear all ;
close all ;
clc ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1. Preliminary calculations and loading data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Parameters
doplot=0; % 0 : plots not displayed ; 1 : plots displayed

%%% Physical constants
g=9.81; % gravity acceleration in m/s^2
rho=1026.; % water density in kg / m^3

%%% Reads data
load deep_waimea.mat;
stname='WAIMEA'; year='2012';source='University of Hawaii';lato=21.671; lono=-158.117;
offbot=0.1; % height of pressure sensor over bottom [m]
zone=0;
fs=1;
nfft=40*5; % nfft specifies the FFT length that csd uses.
df=fs/nfft;                  % Frequential resolution 
numoverlap = nfft/2 ;        % numoverlap is the number of samples by which the sections overlap.
Nf=nfft/2+1;
Nfo=Nf-1;  % up to Nyquist.
freq=linspace(df,df*(Nf-1),Nf-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2. Setting NetCDF file attributes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% where the frequency spectrum will be recorded for later plotting and comparison

%%% Prepares netCDF file for recording
fileoutNC = [stname '_' year '_freq.nc']; % name of netCDF file
fana=netcdf.create(fileoutNC,'CLOBBER'); % opening the corresponding netCDF file handle using Matlab netCDF package
g1=netcdf.getConstant('GLOBAL'); % define shortcut for attribute GLOBAL for variables
NC_UNLIMITED=netcdf.getConstant('NC_UNLIMITED'); % defines shortcut for allowing unlimited variable size
today=datestr(now,'YYYY-mm-DD');
t1=datestr(time(1),'YYYY-mm-DDTHH:MM:SSZ');
t2=datestr(time(end),'YYYY-mm-DDTHH:MM:SSZ');
nt=length(time)*4;
nf=Nfo;

netcdf.putAtt(fana,g1,'id',stname);
netcdf.putAtt(fana,g1,'naming_authority','LOPS');
netcdf.putAtt(fana,g1,'title','Ocean wave spectra, bottom pressure and surface elevation');
netcdf.putAtt(fana,g1,'summary','');
netcdf.putAtt(fana,g1,'cdm_feature_type','station');
netcdf.putAtt(fana,g1,'keywords','');
netcdf.putAtt(fana,g1,'keywords_vocabulary','');
netcdf.putAtt(fana,g1,'standard_name_vocabulary','CF-1.6');
netcdf.putAtt(fana,g1,'scientific_project','IOWAGA');
netcdf.putAtt(fana,g1,'acknowledgement',[ 'Please acknowledge data source ' source ]);
netcdf.putAtt(fana,g1,'license','none');
netcdf.putAtt(fana,g1,'format_version','1.0');
netcdf.putAtt(fana,g1,'history','1.0: Processing from matlab code');
netcdf.putAtt(fana,g1,'publisher_name','LOPS');
netcdf.putAtt(fana,g1,'publisher_url','http://www.umr-lops.fr');
netcdf.putAtt(fana,g1,'publisher_email','fpaf@ifremer.fr');
netcdf.putAtt(fana,g1,'creator_name','Fabrice Ardhuin');
netcdf.putAtt(fana,g1,'creator_url','http://www.umr-lops.fr');
netcdf.putAtt(fana,g1,'creator_email','fabrice.ardhuin@ifremer.fr');
netcdf.putAtt(fana,g1,'processing_software','LOPS wave toolbox');
netcdf.putAtt(fana,g1,'processing_level','0');
netcdf.putAtt(fana,g1,'references','');
netcdf.putAtt(fana,g1,'geospatial_lat_min',lato);
netcdf.putAtt(fana,g1,'geospatial_lat_max',lato);
netcdf.putAtt(fana,g1,'geospatial_lat_units','degrees');
netcdf.putAtt(fana,g1,'geospatial_lon_min',lono);
netcdf.putAtt(fana,g1,'geospatial_lon_max',lono);
netcdf.putAtt(fana,g1,'geospatial_lon_units','degrees');
netcdf.putAtt(fana,g1,'geospatial_vertical_min','');
netcdf.putAtt(fana,g1,'geospatial_vertical_max','');
netcdf.putAtt(fana,g1,'geospatial_vertical_units','meters above mean sea level');
netcdf.putAtt(fana,g1,'geospatial_vertical_positive','up');
netcdf.putAtt(fana,g1,'time_coverage_start',t1);
netcdf.putAtt(fana,g1,'time_coverage_stop',t2);
netcdf.putAtt(fana,g1,'time_coverage_resolution','3 hourly');

%%% Define dimensions of variables
dimida3 = netcdf.defDim(fana,'frequency',nf);
dimida2 = netcdf.defDim(fana,'string16',16);
dimida1 = netcdf.defDim(fana,'station',1);
dimida0 = netcdf.defDim(fana,'time',NC_UNLIMITED);

%%% Define variables and their attributes
vara0=netcdf.defVar(fana,'time','double',[dimida0]);
netcdf.putAtt(fana,vara0,'long_name','julian day (UT)') ;
netcdf.putAtt(fana,vara0,'standard_name','time') ;
netcdf.putAtt(fana,vara0,'units','days since 1990-01-01T00:00:00Z') ;
netcdf.putAtt(fana,vara0,'conventions','Relative julian days with decimal part (as parts of the day)') ;
netcdf.putAtt(fana,vara0,'axis','T') ;

vara1=netcdf.defVar(fana,'station','int',[dimida1]);
netcdf.putAtt(fana,vara1,'axis','X') ;

vara2=netcdf.defVar(fana,'string16','int',[dimida2]);
netcdf.putAtt(fana,vara2,'axis','W') ;

vara3=netcdf.defVar(fana,'station_name','char',[dimida2 dimida1]);
netcdf.putAtt(fana,vara3,'long_name','station_name') ;
netcdf.putAtt(fana,vara3,'contents','XW') ;
netcdf.putAtt(fana,vara3,'associates','string16 station') ;

vara4=netcdf.defVar(fana,'longitude','float',[dimida1 dimida0]);
netcdf.putAtt(fana,vara4,'long_name','longitude') ;
netcdf.putAtt(fana,vara4,'units','degree_east') ;
netcdf.putAtt(fana,vara4,'valid_min',-180.) ;
netcdf.putAtt(fana,vara4,'valid_max',180.) ;
netcdf.putAtt(fana,vara4,'_FillValue',single(9.96921E+36));
netcdf.putAtt(fana,vara4,'contents','XT') ;
netcdf.putAtt(fana,vara4,'associates','station time') ;

vara5=netcdf.defVar(fana,'latitude','float',[dimida1 dimida0]);
netcdf.putAtt(fana,vara5,'long_name','latitude') ;
netcdf.putAtt(fana,vara5,'units','degree_north') ;
netcdf.putAtt(fana,vara5,'valid_min',-90.) ;
netcdf.putAtt(fana,vara5,'valid_max',90.) ;
netcdf.putAtt(fana,vara5,'_FillValue',single(9.96921E+36));
netcdf.putAtt(fana,vara5,'contents','XT') ;
netcdf.putAtt(fana,vara5,'associates','time station') ;

vara6=netcdf.defVar(fana,'frequency','double',[dimida3]);
netcdf.putAtt(fana,vara6,'long_name','frequency of center band');
netcdf.putAtt(fana,vara6,'standard_name','wave_frequency');
netcdf.putAtt(fana,vara6,'globwave_name','frequency');
netcdf.putAtt(fana,vara6,'units','s-1');
netcdf.putAtt(fana,vara6,'valid_min',0.');
netcdf.putAtt(fana,vara6,'valid_max',10.');
netcdf.putAtt(fana,vara6,'axis','Y');

vara7=netcdf.defVar(fana,'frequency1','double',[dimida3]);
netcdf.putAtt(fana,vara7,'long_name','frequency of lower band');
netcdf.putAtt(fana,vara7,'standard_name','frequency_of_lower_band');
netcdf.putAtt(fana,vara7,'globwave_name','frequency_of_lower_band');
netcdf.putAtt(fana,vara7,'units','s-1');
netcdf.putAtt(fana,vara7,'valid_min',0.');
netcdf.putAtt(fana,vara7,'valid_max',10.');
netcdf.putAtt(fana,vara7,'axis','Y');

vara8=netcdf.defVar(fana,'frequency2','double',[dimida3]);
netcdf.putAtt(fana,vara8,'long_name','frequency of upper band');
netcdf.putAtt(fana,vara8,'standard_name','frequency_of_upper_band');
netcdf.putAtt(fana,vara8,'globwave_name','frequency_of_upper_band');
netcdf.putAtt(fana,vara8,'units','s-1');
netcdf.putAtt(fana,vara8,'valid_min',single(0.));
netcdf.putAtt(fana,vara8,'valid_max',10.);
netcdf.putAtt(fana,vara8,'axis','Y');

%vara9=netcdf.defVar(fana,'ef','float',[dimida3 dimida1 dimida0]);
vara9=netcdf.defVar(fana,'ef','float',[dimida3 dimida1 dimida0]);

netcdf.putAtt(fana,vara9,'long_name','wave_elevation_spectrum');
netcdf.putAtt(fana,vara9,'standard_name','wave_elevation_spectrum');
netcdf.putAtt(fana,vara9,'globwave_name','sea_surface_variance_spectral_density');
netcdf.putAtt(fana,vara9,'units','m2 s');
netcdf.putAtt(fana,vara9,'valid_min',single(0.));
netcdf.putAtt(fana,vara9,'valid_max',single(1.e+20));
netcdf.putAtt(fana,vara9,'scale_factor',single(1.));
netcdf.putAtt(fana,vara9,'add_offset',single(0.));
netcdf.putAtt(fana,vara9,'_FillValue',single(9.96921E+36));
netcdf.putAtt(fana,vara9,'content','YXT');
netcdf.putAtt(fana,vara9,'associates','frequency station time') ;

vara10=netcdf.defVar(fana,'epf','float',[dimida3 dimida1 dimida0]);

netcdf.putAtt(fana,vara10,'long_name','seafloor_pressure_spectrum');
netcdf.putAtt(fana,vara10,'standard_name','sea_floor_wave_pressure_spectrum');
netcdf.putAtt(fana,vara10,'globwave_name','sea_floor_wave_pressure_spectrum');
netcdf.putAtt(fana,vara10,'units','m2 s');
netcdf.putAtt(fana,vara10,'valid_min',single(0.));
netcdf.putAtt(fana,vara10,'valid_max',single(1.e+20));
netcdf.putAtt(fana,vara10,'scale_factor',single(1.));
netcdf.putAtt(fana,vara10,'add_offset',single(0.));
netcdf.putAtt(fana,vara10,'_FillValue',single(9.96921E+36));
netcdf.putAtt(fana,vara10,'content','YXT');
netcdf.putAtt(fana,vara10,'associates','frequency station time') ;

vara11=netcdf.defVar(fana,'depth','float',[dimida1 dimida0]);
netcdf.putAtt(fana,vara11,'long_name','latitude') ;
netcdf.putAtt(fana,vara11,'units','m') ;
netcdf.putAtt(fana,vara11,'valid_min',0.) ;
netcdf.putAtt(fana,vara11,'valid_max',12000.) ;
netcdf.putAtt(fana,vara11,'_FillValue',single(9.96921E+36));
netcdf.putAtt(fana,vara11,'contents','XT') ;
netcdf.putAtt(fana,vara11,'associates','time station') ;


netcdf.endDef(fana); 

freq=linspace(df,df*(Nf-1),Nf-1); 
freq1=freq(1:Nfo)-df*0.5;
freq2=freq(1:Nfo)+df*0.5;

netcdf.putVar(fana,vara1,1); 
netcdf.putVar(fana,vara2,linspace(1,16,16)); 
stationname='                ';
ns=length(stname);
stationname(1:ns)=stname;
netcdf.putVar(fana,vara3,stationname); 
netcdf.putVar(fana,vara6,freq); 
netcdf.putVar(fana,vara7,freq1); 
netcdf.putVar(fana,vara8,freq2); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 3. Compute surface elevation spectrum and record it in netCDF file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nt0=3600*3; % seconds in 3 hours
nrec=209*4;
irec=0;
%%% Loop over times
while (irec < nrec) % loop over time
    jrec=floor(irec/4);
    krec=(irec-jrec*4);
    imax=min([43180,nt0*(krec+1)]);
    nt=imax-krec*nt0;
    elevation = hh(krec*nt0+1:imax,jrec+1)*0.01; %% hh : variable in deep_waimea.mat matrix ; corresponds to pressure divided by rho*g [m]. The 0.01 comes from the conversion from cm to m.
    avg=mean(elevation);

    now=time(jrec+1)+krec/8; %
    dates2=datevec(now);
    msg=sprintf('Reading data for time: %d %d %d %d %d \n', ... 
    dates2(1),dates2(2),dates2(3),dates2(4),dates2(5));
    fprintf('%s \n',msg)

    See=zeros(Nf,1);          % Power spectrum 
    Ef=zeros(Nf-1,1);        % Power spectrum density

    hanning=transpose(0.5 * (1-cos(2*pi*linspace(0,nfft-1,nfft)/(nfft-1))));
    NS1=floor(nt/nfft);
    NS=NS1*2-1;
    elevmat=zeros(nfft,NS);
    elevmat(1:nfft,1:NS1)=reshape(elevation(1:nfft*NS1),[nfft,NS1]);
    elevmat(1:nfft,NS1+1:NS)=reshape(elevation(1+nfft/2:nfft/2+nfft*(NS1-1)),[nfft,NS1-1]);
    elevmat2=detrend(elevmat,'linear');             % We remove the linear trend from each column
   
    if doplot==1
        %PLOT OF DETRENDED TIME SERIES
        figure(3)
        t=linspace(0,length(elevation)-1,length(elevation))/fs;
        plot(t,elevation,'k-')
        xlabel('time(s)');
        ylabel('pressure(m)');
        grid
        pause
        figure(4)
        plot(t(1:NS1*nfft),reshape(elevmat2(:,1:NS1),NS1*nfft,1),'k-')
        legend('p-pmean');
        xlabel('time(s)');
        grid
        pause
    end
   
    Hs1=4.*sqrt(mean(var(elevmat2)));
    Eh=ones(1,NS);
    H=hanning*Eh;
   
   NSbad=0;
   nstep=4;  % checks for bad samples
   while (nstep > 0) 
      nstep=nstep-1; 
      var2=var(elevmat2);
      meanvar=mean(var2);
      index=find(var2-3*meanvar > 0);
      taille=size(index);
      NSbad=NSbad+taille(2);
      if (NSbad == 0) 
          nstep=0;
      end
      elevmat2(:,index)=0;
   end
   elevmat2=H.*detrend(elevmat2,'constant');             % remove mean and apply window
   Hs2=4.*sqrt(mean(var(elevmat2)));
   wc2=1/mean(hanning.^2);                              % window correction factor
   fac=NS/(NS-NSbad);
   
   pspec=fft(elevmat2,nfft,1)/nfft;
  
   Sppmat=(abs(pspec).^2)*wc2/df;
   
    
    
   freq1=linspace(0,df*(nfft-1),nfft);
   %plot(freq1,transpose(Seemat(1,:)))
   dep=mean(elevmat,1)+(offbot);
  
   d=offbot;
   i=0;
   
   % converts two sided spectral density to one sided spectral density
   Sppmat=Sppmat.*2; 
   Seemat=Sppmat;
   % the depth correction is only applied when the correction coefficient
   % is less than a threshold, and only if the signal is strong enough
   indcut=1;
   while i < NS
      i=i+1;
      wn=dispNewtonTH(freq1(2:nfft),dep(i));
      cosur=transpose((cosh(wn*dep(i))./cosh(wn.*d)).^2);

      indc=find(cosur < 1000);
      indc=indc(end);
      taille=size(indc);
      if taille(1) > 0
      Seemat(2:indc,i)=Sppmat(2:indc,i).*cosur(1:indc-1);
      end
   end
   fc=freq1(indc+1);
   
   
   
   tratio=2.*pi.*freq1(2:nfft)./tanh(wn.*mean(dep));
   depth=mean(dep);
   
 


   See=mean(Seemat,2)*fac;
   Spp=mean(Sppmat,2)*fac;
 
   Ef=transpose(See(2:Nf));
   Ep=transpose(Spp(2:Nf));
   
   if doplot==1
   %MPLOT OF FREQUENCY SPECTRA
   figure(1)
   hp=loglog(freq(1:(Nf-1)),Sppmat(2:Nf,:));
      axis([0,0.5,0.0000001,100.]);
        set(hp,'LineWidth',2);
   xlabel('f (Hz)');
   ylabel('Bottom pressure spectrum (m^2/Hz)');
   grid;
   pause

   %PLOT OF MEAN FREQUENCY SPECTRUM
   figure(2)
   hp=loglog(freq(1:(Nf-1)),Ef(1:(Nf-1)), ...
          freq(1:(Nf-1)),0.001*freq(1:(Nf-1)).^-4,'r-');
      axis([0,0.5,0.0000001,100.]);
        set(hp,'LineWidth',2);
        legend('Elevation spectrum','f^{-4}');
   xlabel('f (Hz)');
   ylabel('surface elevation spectrum (m^2/Hz)');
   grid;
   pause
   end
      

Hs=4.*sqrt(sum(Ef(:))*df);
      Hsband=4.*sqrt(sum(Ef(6:(Nf-1)/2))*df);
      HsIG=4.*sqrt(sum(Ef(1:5))*df);
      [Emax,ifp]=max(Ef(1:(Nf-1)/2));


      if (ifp ~= 1)  
         fp=sum(freq(ifp-1:ifp+1).*Ef(ifp-1:ifp+1)) / sum(Ef(ifp-1:ifp+1));
         %sth1p=(2*(1 - ((sum(Ef(ifp-1:ifp+1).^2.*(a1(ifp-1:ifp+1).^2 + b1(ifp-1:ifp+1).^2))).^0.5/sum(Ef(ifp-1:ifp+1))))).^0.5*180/pi;
      else    
         fp=freq(ifp);
      end
      f01=(sum(Ef(1:Nfo).*(freq(1:Nfo).^1))/sum(Ef(1:Nfo)));
      f0m2=(sum(Ef(1:Nfo).*(freq(1:Nfo).^-2))/sum(Ef(1:Nfo)))^-0.5;
      f02=(sum(Ef(1:Nfo).*(freq(1:Nfo).^2))/sum(Ef(1:Nfo)))^0.5;;
%       msg=sprintf('irec= %d , depth, Hs :%6.2f %6.2f %6.2f  %7.2f %7.2f, errors : %d %d \n', ... 
%       irec,depth,Hs,1/fp,NSbad,fac);
%       fprintf('%s \n',msg)


      netcdf.putVar(fana,vara0,irec,1,now-datenum(1990,1,1)); 
      netcdf.putVar(fana,vara9,[0 0 irec ],[Nfo 1 1 ],Ef(1:Nfo));
      netcdf.putVar(fana,vara10,[ 0 0 irec],[Nfo 1  1 ],Ep(1:Nfo));
      netcdf.putVar(fana,vara4,[0 irec],[1 1],lono); 
      netcdf.putVar(fana,vara5,[0 irec],[1 1],lato); 
      netcdf.putVar(fana,vara11,[0 irec],[1 1],depth); 
   irec=irec+1;
   end
   
netcdf.close(fana);



