function[err]= spectralekman(N)

%*****************************************************************************80
%
%% SPECTRALEKMAN ?
%
%  Discussion:
%
%    This function is discussed in chapter 8 of the reference.
%
%  Modified:
%
%    01 December 2009
%
%  Author:
%
%    Benoit Cushman-Roisin, Jean-Marie Beckers.
%
%  Reference:
%
%    Benoit Cushman-Roisin, Jean-Marie Beckers,
%    Introduction to Geophysical Fluid Dynamics,
%    Academic Press, 2010.
%
%  Parameters:
%
  nu=0.01;
  h=100;
  f=0.0001;
  tauy=0;
  taux=10^(-6)*10*10*0;
  pressy=9.81*1/1000000;
  pressx=0;
  KM=400;
  dek=sqrt(2*nu/f);

% RESOLUTION and PLOT

  for j=1:N
    s(j)=2*sqrt(2*h)/(pi*(2*j-1));
    phih(j)=sqrt(2/h)*(-1)^(j-1);
    eigen(j)=(2*j-1)^2*pi*pi*nu/(4*h*h);
    FX(j)=-pressx*s(j)+taux*phih(j);
    FY(j)=-pressy*s(j)+tauy*phih(j);
    a(j)=(eigen(j)*FX(j)+f*FY(j))/(eigen(j)^2+f^2);
    b(j)=(eigen(j)*FY(j)-f*FX(j))/(eigen(j)^2+f^2);
  end

  for k=1:KM
    U(k)=0;
    V(k)=0;
  end

  for k=1:KM
    z(k)=(k-1)*h/(KM-1);
    U(k)=0;
    V(k)=0;
    uex(k)=0;
    vex(k)=0;
    for j=N:-1:1
      U(k)=U(k)+a(j)*sqrt(2/h)*sin( (2*j-1)*pi/2*z(k)/h);
      V(k)=V(k)+b(j)*sqrt(2/h)*sin( (2*j-1)*pi/2*z(k)/h);
    end
  end

% plot(U,z,V,z)
% pause
% plot(U,V)
% axis equal
% pause
% exaxt?

  for k=1:KM
    [uex(k),vex(k)]=ekmanh(h/dek,z(k)/dek);
    %uex(k)=uex(k)*taux/nu*dek;
    %vex(k)=vex(k)*taux/nu*dek;
    uex(k)=-(1+uex(k))*pressy/f;
    vex(k)=-vex(k)*pressy/f;
  end

  plot(uex,z,vex,z,U,z,'-.',V,z,'-.')
  xlabel('velocity');
  ylabel('z');

  pause

  err=0;
  for k=1:KM
    err=err+(U(k)-uex(k))^2+(V(k)-vex(k))^2;
  end
  err=sqrt(err/KM);
  for kk=1:N;
    aa(kk)=log(a(kk)*a(kk));
    xx(kk)=log(kk);
  end
  bar(xx,aa+40);
  xlabel('log(j)');
  ylabel('log(a_j^2)');

  pause

  return
end