% Demo: Tidal Harmonic Analysis

% 1. load data;
 load p20005f126.txt;
 y=p20005f126;
 % t in hour unit, w in 1/hr
 
 % 2. remove low frequency sea level variation,
 
  t=y(100:600,3)*24;
 ellow=filterhan(y(:,4),51);
 p1=y(100:600,4)-ellow(100:600); 
 figure(1)
 subplot(311)
 plot(t,y(100:600,4)); legend('field data')
 subplot(312)
 plot(t,p1); legend ('only tides') 
 subplot(313)
 plot(t,ellow(100:600)); legend ('low frequency sea level variation')
 
 % Harmonic components convert Deg to Rad
 % Angular frequency: w=(360)/T (in de)gree
 w=[28.984104;28.439729;30.00000;15.041068;13.943035]*(pi/180);
 %  M2,        N2       , S2,     K1        S1
 
 %-- find harmonic component
 
 [am phase a0]=harmon(p1,t, w);
 
 % -- find predicted tide using Harmonic analysis
  
 tpred=y(100:700,3)*24;
 len_predict=length(tpred);
 p1h=zeros(len_predict,1);

 n=length(w);
  for i=1:len_predict;
    el0=0;
    for j=1:n;
       el0=am(j)*cos(w(j)*tpred(i)-phase(j))+el0;
    end
    p1h(i)=a0+el0;    
 end
 
 % --- plot the results
 figure (2)
 
 subplot(311)
 plot(t,p1,'.b',tpred,p1h,'r')
 axis([100 700 -1 1])
 legend('observations','predictions')
 
 subplot(312)
 plot(p1,p1h(1:length(t)),'.')
  legend('correlation between observations and predictions')
  % statistics
  p1h2=p1h(1:length(t));
  r=corrcoef(p1,p1h2)
  ne=length(p1);
  'correlation coef ='
 r=corrcoef(p1,p1h2)
 ne=length(p1);
 'root mean square error='
 rmse= sqrt(sum(p1-p1h2).^2/ne)
