function [fit,cov_fit,data_hat,chi_2_misfit]=coord_model_hadley(time, ...
          data,C,rate_flag,yr_flag,syr_flag,offsets,epoch_tau,TAUs,VBREAKs)

%    [fit,cov_fit,data_hat,chi_2_misfit]=coord_model(time, ...
%                                        data,C,rate_flag,yr_flag,syr_flag,offsets,epoch_tau,TAUs,VBREAKs)
%
%    time        vector of x axis values (time for instance)
%    data        vector of y axis values (measured data)
%    C           full covariance matrix of measured data
%    rate_flag   [0 1] estimate a rate
%    yr_flag     [0 1] estimate a yearly signal (phase wrt Jan.1)
%                (adds 2 parameters to the model)
%    syr_flag    [0 1] estimate a semi-annual signal (phase wrt Jan.1)
%                (adds 2 parameters to the model)
%    offsets     (nx1) vector of offset times.
%    epoch_tau   (dx1) vector of epochs for exponential decays
%    TAUs        (dx1) vector of exponential decay time constants
%    VBREAKs     (vx1) vector of velocity break times.  A new rate 
%                will be estimated starting at each break.
%
%    fit         vector of [intercept; slope(s); (sine); (cosine); (offset(s)); (decay(s))]
%    cov_fit     covariance matrix of fitted parameters
%    data_hat    predicted data
%    chi_2_misfit chi square test statistic computed in the rotated space

%    Hadley Johnson, hjohnson@ucsd.edu, 23-Sep-1998 version.
%    modified by Rosanne Nikolaidis, Nov 2000 to include exponential decay and velocity breaks

%following line can be used to minimize memory usage
% global C

time=time(:);
%time1=time(1);
%time_yr = time - mod(time1,1);
format long
data=data(:);
n=length(time);
C=chol(C);
C=C';
New_d=C\data;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Build the A matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A=[ones(n,1)];

if rate_flag==1,
  A=[A time];
  nvb = length(VBREAKs);
  if nvb>0
     mask=time<VBREAKs(1);
     A(~mask,2)=VBREAKs(1)*ones(sum(~mask),1);
     for i=1:nvb-1,
%       A=[A time>VBREAKs(i)];
       mask1=time<VBREAKs(i);
       mask2=time>VBREAKs(i+1);
       A=[A [zeros(sum(mask1),1);time(~mask1&(~mask2))-VBREAKs(i); ...
          ones(sum(mask2),1)*VBREAKs(i+1)-VBREAKs(i)]];
     end
     lastvb = VBREAKs(nvb);
%     A=[A time>lastvb];
     mask=time>lastvb;
     A=[A mask.*(time-lastvb)];
  end
end

if yr_flag==1,
%  A=[A sin(time_yr*2*pi) cos(time_yr*2*pi)];
   A=[A sin(time*2*pi) cos(time*2*pi)];
end

if syr_flag==1,
%  A=[A sin(time_yr*4*pi) cos(time_yr*4*pi)];
  A=[A sin(time*4*pi) cos(time*4*pi)];
end

for i=1:length(offsets),
  A=[A time>=offsets(i)];
end

for i=1:length(TAUs),
  if TAUs(i)~=0,
     A=[A (time>=epoch_tau(i)).*exp(-(time-epoch_tau(i))./TAUs(i))];
  end
end

p=length(A(1,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

New_A=C\A;
[C,R]=qr(New_A);
R1=R(1:p,1:p);
t1=[inv(R1) zeros(p,n-p)]*C';
cov_fit=t1*t1';
fit=t1*New_d;
data_hat=A*fit;

New_d_hat=New_A*fit;
chi_2_misfit=sum((New_d-New_d_hat).^2); %computed in the rotated space where
                                        %noise is N(0,1).

