function OK = pea_simpson_et_al_1990_sigma_test
%pea_simpson_et_al_1990_sigma_test unit test for pea_simpson_et_al_1990_sigma
%
% For the 3 analytical test euqtions see appendix 1 of
%
% * F.C. Groenendijk, 1992. Distribution in time and space of fresh water
%   and suspended matter in the Dutch coastal zone. Directoraat-Generaal 
%   Rijkswaterstaat, Dienst Getijdewateren, 1992. <a href="http://books.google.com/books?id=SGAhGwAACAAJ">Google books</a>.
%
%See also: pea_simpson_et_al_1990_sigma, pea_simpson_et_al_1990

   drho               = 10;
   H                  = 20;
   zcen_intface       = 0:-.1:-H;
   h                  = abs(zcen_intface(1)-zcen_intface(end));
   zcen_centres       = corner2center1(zcen_intface);

   zcen_intface2      = repmat(zcen_intface,[2 1]);
   
   dsigma             = abs(diff(zcen_intface)./H);
   
   OK = [];

%% LINEAR testcase Groenendijk, 1992

   L.n                  = length(zcen_centres);
   L.density            = 0 - drho.*zcen_centres./H;
   L.density2           = repmat(L.density     ,[2 1]);
   L.densitymean        = meanweighted(L.density,2,abs(diff(zcen_intface)));
   
   L.phi  = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H, flipdim(L.density ,2));
   L.phi1 = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H, flipdim(L.density2,2));
   
   L.phi_ok  = 9.81.*drho.*H./12;
   OK(end+1) = abs(L.phi  -L.phi_ok ) < L.phi./1e3;
   OK(end+1) = all(L.phi1 ==L.phi);
   
%% TWO-LAYER testcase Groenendijk, 1992

   T.density            = zeros(size(zcen_centres));
   T.density(zcen_centres < -H/2) = drho;
   T.density2           = repmat(T.density     ,[2 1]);
   T.densitymean        = meanweighted(T.density,2,abs(diff(zcen_intface)));
   
   T.phi  = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H, flipdim(T.density ,2));
   T.phi1 = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H, flipdim(T.density2,2));
   
   T.phi_ok  = 9.81.*drho.*H./8;
   OK(end+1) = abs(T.phi -T.phi_ok ) < T.phi./1e3;
   OK(end+1) = all(T.phi1==T.phi);
   
%% SINE-LIKE testcase Groenendijk, 1992

   S.n                  = length(zcen_centres);
   S.density            = 0 + drho./2.*(1-cos(zcen_centres*pi./H));
   S.density2           = repmat(S.density     ,[2 1]);
   S.densitymean        = meanweighted(S.density,2,abs(diff(zcen_intface)));
   
   S.phi  = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H,flipdim(S.density ,2));
   S.phi1 = pea_simpson_et_al_1990_sigma(dsigma(end:-1:1), 0, -H,flipdim(S.density2,2));
   
   S.phi_ok  = 9.81.*drho.*H./pi^2;
   OK(end+1) = abs(S.phi - S.phi_ok ) < S.phi./1e3;
   OK(end+1) = all(S.phi1==S.phi);

   OK = all(OK);
   
   if ~OK
       L
       T
       S
   end
