%
% fits the following disturbance model to the given
% time series and plots the results
%
% y = P1 + P2*( sin(2*pi*t/24 + P4*pi) +
%                P3*sin(2*pi*t/12 + (P4+P5)*pi) )
%
%    [ p, sos ] = fitdist( y, delt, pflag, [fixp, po] )
%
% where   p     = parameters fitted
%         sos   = final sum-of-squares
%         y     = time series
%         delt  = sampling time interval
%         pflag = produces plots if nonzero
%         fixp  = row vector of five flags indicating,
%                 if positive nonzero, that the corresponding 
%                 parameter should be considered fixed
%                 at the value specified in "po"
%         po    = row vector of five parameters used
%                 as fixed parameters (see "fixp") or
%                 as initial values
%
% Bob Newell, April 1996
%
%----------------------------------------------------------
function [ p, sos ] = fitdist( y, delt, pflag, fixp, po )
%
global t
global yact
global r
global pfix
global pp
%
n = length( y );
t = delt * ( 0:n-1 )'; 
yact = y(:);
%
pfix = zeros(1,5);
if nargin < 4,
  pp = zeros(1,5);
  pest = [ 1 1 1 1 1 ];
else,
  pp = po;
  pest = [];
  for i = 1:5,
    if fixp(i) > 0,
	  pfix(i) = 1;
    else,
	  pest = [ pest po(i) ];
	end;
  end;
end;
%
if length(pest) > 0,
  p = fmins( 'fitdistf', pest, [0 .1 .1] )
else,
  p = po;
end;
%
sos = fitdistf( p );
sost = ( yact - mean(yact) )' * ( yact - mean(yact) );
disp( [ 'Residuals: mean = ' num2str( mean( r ) ) ...
             ' +- ' num2str( std( r ) ) ] );
disp( [ 'Sum of squares = ' num2str( sos ) ] );
disp( [ 'Percent explained = ' ...
           num2str( 100 * ( sost - sos ) / sost ) ] );
%
if pflag,
  figure;
  plot( t, yact-r, '-', t, yact, '*' );
  xlabel( 'Time' );
  ylabel( 'Data' );
  figure;
  plot( t, r, '-' );
  xlabel( 'Time' );
  ylabel( 'Residuals' );
  ax = axis;
  line( [ ax(1) ax(2) ], [ 0 0 ] );
end;
%
return;
%----------------------------------------------------------
