end;
    
  coef=lhs\rhs;
  
  z0=coef(1);
  ap=(coef(2:(1+mu))-i*coef((2+mu):(1+2*mu)))/2;  % a+ amplitudes
  am=(coef(2:(1+mu))+i*coef((2+mu):(1+2*mu)))/2;  % a- amplitudes
  if secular(1:3)=='lin',
    dz0=coef(end);
  else
    dz0=0;
  end; 
  
  xout=xin; % Copies over NaN   
  if secular(1:3)=='lin',
    for j1=1:nsub:nobs
      j2=min(j1 + nsub - 1,nobs);
      E=[ones(j2-j1+1,1) cos((2*pi)*t(j1:j2)*fu') sin((2*pi)*t(j1:j2)*fu') t(j1:j2)*(2/dt/nobsu)];
      xout(j1:j2)=E*coef;
    end;
  else  
    for j1=1:nsub:nobs
      j2=min(j1 + nsub - 1,nobs);
      E=[ones(j2-j1+1,1) cos((2*pi)*t(j1:j2)*fu') sin((2*pi)*t(j1:j2)*fu')];
      xout(j1:j2)=E*coef;
    end;
  end;

end;

   
 
%----------------------------------------------------------------------
% Check variance explained (but do this with the original fit).

xres=xin-xout; % and the residuals!

if isreal(xin),    % Real time series
  varx=cov(xin(gd));varxp=cov(xout(gd));varxr=cov(xres(gd));
  fprintf('   percent of var residual after lsqfit/var original: %5.2f %%\n',100*(varxr/varx));  
else               % Complex time series
  varx=cov(real(xin(gd)));varxp=cov(real(xout(gd)));varxr=cov(real(xres(gd)));
  fprintf('   percent of X var residual after lsqfit/var original: %5.2f %%\n',100*(varxr/varx));

  vary=cov(imag(xin(gd)));varyp=cov(imag(xout(gd)));varyr=cov(imag(xres(gd)));
  fprintf('   percent of Y var residual after lsqfit/var original: %5.2f %%\n',100*(varyr/vary));
end;


%---------- Correct for prefiltering-----------------------------------

corrfac=interp1(corr_fs,corr_fac,fu);
% To stop things blowing up!
corrfac(corrfac>100 | corrfac <.01 | isnan(corrfac))=1;

ap=ap.*corrfac;
am=am.*conj(corrfac);

%---------------Nodal Corrections-------------------------------------- 						   
% Generate nodal corrections and calculate phase relative to Greenwich. 						   
% Note that this is a slightly weird way to do the nodal corrections,							   
% but is 'traditional'.  The "right" way would be to change the basis							   
% functions used in the least-squares fit above.									   

if ~isempty(lat) & ~isempty(stime),   % Time and latitude								   

  % Get nodal corrections at midpoint time.										   
  [v,u,f]=t_vuf(ltype,centraltime,[ju;jinf],lat);									   

  vu=(v+u)*360; % total phase correction (degrees)									   
  nodcor=['Greenwich phase computed with nodal corrections applied to amplitude \n and phase relative to center time'];    
elseif ~isempty(stime),    % Time only  										   
  % Get nodal corrections at midpoint time										   
  [v,u,f]=t_vuf(ltype,centraltime,[ju;jinf]);										   
  vu=(v+u)*360; % total phase correction (degrees)									   
  nodcor=['Greenwich phase computed, no nodal corrections'];								   
else   % No time, no latitude												   
  vu=zeros(length(ju)+length(jinf),1);  										   
  f=ones(length(ju)+length(jinf),1);											   
   nodcor=['Phases at central time'];											   
end															   
fprintf(['   ',nodcor,'\n']);												   


%---------------Inference Corrections----------------------------------
% Once again, the "right" way to do this would be to change the basis
% functions.
ii=find(finite(jref));
if ii,
  fprintf('   Do inference corrections\n');
  snarg=nobsu*pi*(fi(ii)   -fu(jref(ii)) )*dt;
  scarg=sin(snarg)./snarg;
 
  if size(inf.amprat,2)==1,    % For real time series
    pearg=     2*pi*(vu(mu+ii)-vu(jref(ii))+inf.ph(ii))/360;
    pcfac=inf.amprat(ii).*f(mu+ii)./f(jref(ii)).*exp(i*pearg);
    pcorr=1+pcfac.*scarg;
    mcfac=conj(pcfac);
    mcorr=conj(pcorr);
  else                          % For complex time series
    pearg=     2*pi*(vu(mu+ii)-vu(jref(ii))+inf.ph(ii,1))/360;
    pcfac=inf.amprat(ii,1).*f(mu+ii)./f(jref(ii)).*exp(i*pearg);
    pcorr=1+pcfac.*scarg;
    mearg=    -2*pi*(vu(mu+ii)-vu(jref(ii))+inf.ph(ii,2))/360;
    mcfac=inf.amprat(ii,2).*f(mu+ii)./f(jref(ii)).*exp(i*mearg);
    mcorr=1+mcfac.*scarg;
  end;
    
  ap(jref(ii))=ap(jref(ii))./pcorr;   % Changes to existing constituents
  ap=[ap;ap(jref(ii)).*pcfac];        % Inferred constituents

  am(jref(ii))=am(jref(ii))./mcorr;
  am=[am;am(jref(ii)).*mcfac];

  fu=[fu;fi(ii)];
  nameu=[nameu;namei(ii,:)];
end;

% --------------Error Bar Calculations---------------------------------
%
% Error bar calcs involve two steps:
%      1) Estimate the uncertainties in the analyzed amplitude
%         for both + and - frequencies (i.e., in 'ap' and 'am').
%         A simple way of doing this is to take the variance of the
%         original time series and divide it into the amount appearing
%         in the bandwidth of the analysis (approximately 1/length).
%         A more sophisticated way is to assume "locally white"
%         noise in the vicinity of, e.g., the diurnal consistuents.
%         This takes into account slopes in the continuum spectrum.
%
%      2) Transform those uncertainties into ones suitable for ellipse
%         parameters (axis lengths, angles). This can be done 
%         analytically for large signal-to-noise ratios. However, the 
%         transformation is non-linear at lows SNR, say, less than 10
%         or so.
%

xr=fixgaps(xres); % Fill in "internal" NaNs with linearly interpolated
                  % values so we can fft things.
nreal=1;

if strmatch(errcalc(2:end),'boot'),
  fprintf('   Using nonlinear bootstrapped error estimates\n');
  
  % "noise" matrices are created with the right covariance structure
  % to add to the analyzed components to create 'nreal' REPLICATES. 
  % 

  nreal=300;             % Create noise matrices 
  [NP,NM]=noise_realizations(xr(finite(xr)),fu,dt,nreal,errcalc);
      
  % All replicates are then transformed (nonlinearly) into ellipse 
  % parameters.  The computed error bars are then based on the std
  % dev of the replicates.

  AP=ap(:,ones(1,nreal))+NP;        % Add to analysis (first column
  AM=am(:,ones(1,nreal))+NM;        % of NM,NP=0 so first column of
                                    % AP/M holds ap/m).
  epsp=angle(AP)*180/pi;            % Angle/magnitude form:
  epsm=angle(AM)*180/pi;
  ap=abs(AP);
  am=abs(AM);
elseif strmatch(errcalc,'linear'),
  fprintf('   Using linearized error estimates\n');
  %
  % Uncertainties in analyzed amplitudes are computed in different
  % spectral bands. Real and imaginary parts of the residual time series
  % are treated separately (no cross-covariance is assumed).
  %
  % Noise estimates are then determined from a linear analysis of errors,
  % assuming that everything is uncorrelated. This is OK for scalar time
  % series but can fail for vector time series if the noise is not 
  % isotropic.
  
  [ercx,eicx]=noise_stats(xr(finite(xr)),fu,dt);
  % Note - here we assume that the error in the cos and sin terms is 
  % equal, and equal to total power in the encompassing frequency bin. 
  % It seems like there should be a factor of 2 here somewhere but it 
  % only works this way! <shrug>
  [emaj,emin,einc,epha]=errell(ap+am,i*(ap-am),ercx,ercx,eicx,eicx);

  epsp=angle(ap)*180/pi;
  epsm=angle(am)*180/pi;
  ap=abs(ap);
  am=abs(am);
else
  error(['Unrecognized type of error analysis: ''' errcalc ''' specified!']);
end;

%-----Convert complex amplitudes to standard ellipse parameters--------

aap=ap./f(:,ones(1,nreal));	% Apply nodal corrections and
aam=am./f(:,ones(1,nreal));	% compute ellipse parameters.

fmaj=aap+aam;                   % major axis
fmin=aap-aam;                   % minor axis

gp=mod( vu(:,ones(1,nreal))-epsp ,360); % pos. Greenwich phase in deg.
gm=mod( vu(:,ones(1,nreal))+epsm ,360); % neg. Greenwich phase in deg.

finc= (epsp+epsm)/2;
finc(:,1)=mod( finc(:,1),180 ); % Ellipse inclination in degrees
				% (mod 180 to prevent ambiguity, i.e., 
				% we always ref. against northern 
				% semi-major axis.
	
finc=cluster(finc,180); 	% Cluster angles around the 'true' 
                                % angle to avoid 360 degree wraps.

pha=mod( gp+finc ,360); 	% Greenwich phase in degrees.

pha=cluster(pha,360);		% Cluster angles around the 'true' angle
				% to avoid 360 degree wraps.

%----------------Generate 95% CI---------------------------------------
%% For bootstrapped errors, we now compute limits of the distribution.
if strmatch(errcalc(2:end),'boot'),
     %% std dev-based estimates.
     % The 95% CI are computed from the sigmas
     % by a 1.96 fudge factor (infinite degrees of freedom).
     % emaj=1.96*std(fmaj,0,2);
     % emin=1.96*std(fmin,0,2);
     % einc=1.96*std(finc,0,2);
     % epha=1.96*std(pha ,0,2);
     %% Median-absolute-deviation (MAD) based estimates.
     % (possibly more stable?)
      emaj=median(abs(fmaj-median(fmaj,2)*ones(1,nreal)),2)/.6375*1.96;
      emin=median(abs(fmin-median(fmin,2)*ones(1,nreal)),2)/.6375*1.96;
      einc=median(abs(finc-median(finc,2)*ones(1,nreal)),2)/.6375*1.96;
      epha=median(abs( pha-median( pha,2)*ones(1,nreal)),2)/.6375*1.96;
else
   % In the linear analysis, the 95% CI are computed from the sigmas
   % by this fudge factor (infinite degrees of freedom).
   emaj=1.96*emaj;
   emin=1.96*emin;
   einc=1.96*einc;
   epha=1.96*epha;
end;
				  
if isreal(xin),
    tidecon=[fmaj(:,1),emaj,pha(:,1),epha];
else
    tidecon=[fmaj(:,1),emaj,fmin(:,1),emin, finc(:,1),einc,pha(:,1),epha];
end;

% Sort results by frequency (needed if anything has been inferred since 
% these are stuck at the end of the list by code above).
if any(finite(jref)),
 [fu,I]=sort(fu);
 nameu=nameu(I,:);
 tidecon=tidecon(I,:);
end;

snr=(tidecon(:,1)./tidecon(:,2)).^2;  % signal to noise ratio

%--------Generate a 'prediction' using significant constituents----------
xoutOLD=xout;
if synth>=0,
 if ~isempty(lat) & ~isempty(stime),
   fprintf('   Generating prediction with nodal corrections, SNR is %f\n',synth);
   xout=t_predic(stime+[0:nobs-1]*dt/24.0,nameu,fu,tidecon,'lat',lat,'synth',synth,'anal',ltype);
 elseif ~isempty(stime), 
   fprintf('   Generating prediction without nodal corrections, SNR is %f\n',synth);
   xout=t_predic(stime+[0:nobs-1]*dt/24.0,nameu,fu,tidecon,'synth',synth,'anal',ltype);
 else
   fprintf('   Generating prediction without nodal corrections, SNR is %f\n',synth);
   xout=t_predic(t/24.0,nameu,fu,tidecon,'synth',synth,'anal',ltype);
 end;
else
 fprintf('   Returning fitted prediction\n');
end;

%----------------------------------------------------------------------
% Check variance explained (but now do this with the synthesized fit).
xres=xin(:)-xout(:); % and the residuals!

%error;

if isreal(xin),    % Real time series
  varx=cov(xin(gd));varxp=cov(xout(gd));varxr=cov(xres(gd));
  fprintf('   percent of var residual after synthesis/var original: %5.2f %%\n',100*(varxr/varx));  
else               % Complex time series
  varx=cov(real(xin(gd)));varxp=cov(real(xout(gd)));varxr=cov(real(xres(gd)));
  fprintf('   percent of X var residual after synthesis/var original: %5.2f %%\n',100*(varxr/varx));

  vary=cov(imag(xin(gd)));varyp=cov(imag(xout(gd)));varyr=cov(imag(xres(gd)));
  fprintf('   percent of Y var residual after synthesis/var original: %5.2f %%\n',100*(varyr/vary));
end;


%-----------------Output results---------------------------------------

if fid>1,
 fprintf(fid,'\n%s\n',['file name: ',filen]);
elseif fid==1,
 fprintf(fid,'-----------------------------------\n');
end

if fid>0,
  fprintf(fid,'date: %s\n',date);
  fprintf(fid,'nobs = %d,  ngood = %d,  record length (days) = %.2f\n',nobs,ngood,length(xin)*dt/24);
  if ~isempty(stime); fprintf(fid,'%s\n',['start time: ',datestr(stime)]); end
  fprintf(fid,'rayleigh criterion = %.1f\n',ray);
  fprintf(fid,'%s\n',nodcor);
%  fprintf(fid,'\n     coefficients from least squares fit of x\n');
%  fprintf(fid,'\n tide    freq        |a+|       err_a+      |a-|       err_a-\n');
%  for k=1:length(fu);
%    if ap(k)>eap(k) | am(k)>eam(k), fprintf('*'); else fprintf(' '); end;
%    fprintf(fid,'%s  %8.5f  %9.4f  %9.4f  %9.4f  %9.4f\n',nameu(k,:),fu(k),ap(k),eap(k),am(k),eam(k));
%  end
  fprintf(fid,'\nx0= %.3g, x trend= %.3g\n',real(z0),real(dz0));
  fprintf(fid,['\nvar(x)= ',num2str(varx),'   var(xp)= ',num2str(varxp),'   var(xres)= ',num2str(varxr) '\n']);
  fprintf(fid,'percent var predicted/var original= %.1f %%\n',100*varxp/varx);

  if isreal(xin)
    fprintf(fid,'\n     tidal amplitude and phase with 95%% CI estimates\n');
    fprintf(fid,'\ntide   freq       amp     amp_err    pha    pha_err     snr\n');
    for k=1:length(fu);
      if snr(k)>synth, fprintf(fid,'*'); else fprintf(fid,' '); end;
      fprintf(fid,'%s %9.7f %9.4f %8.3f %8.2f %8.2f %8.2g\n',nameu(k,:),fu(k),tidecon(k,:),snr(k));
    end
  else
    fprintf(fid,'\ny0= %.3g, x trend= %.3g\n',imag(z0),imag(dz0));
    fprintf(fid,['\nvar(y)= ',num2str(vary),'    var(yp)= ',num2str(varyp),'  var(yres)= ',num2str(varyr) '\n']);
    fprintf(fid,'percent var predicted/var original= %.1f %%\n',100*varyp/vary);
    fprintf(fid,'\n%s\n',['ellipse parameters with 95%% CI estimates']);
    fprintf(fid,'\n%s\n',['tide   freq      major  emaj    minor   emin     inc    einc     pha    epha      snr']);
    for k=1:length(fu);
      if snr(k)>synth, fprintf(fid,'*'); else fprintf(fid,' '); end;
      fprintf(fid,'%s %9.7f %6.3f %7.3f %7.3f %6.2f %8.2f %6.2f %8.2f %6.2f %6.2g\n',...
	nameu(k,:),fu(k),tidecon(k,:),snr(k));
    end
    fprintf(fid,['\ntotal var= ',num2str(varx+vary),'   pred var= ',num2str(varxp+varyp) '\n']);
    fprintf(fid,'percent total var predicted/var original= %.1f %%\n\n',100*(varxp+varyp)/(varx+vary));
  end

  if fid~=1, st=fclose(fid); end
end;

xout=reshape(xout,inn,inm);
switch nargout,
  case {0,3,4}
  case {1}
   nameu = struct('name',nameu,'freq',fu,'tidecon',tidecon,'type',ltype);
  case {2}   
   nameu = struct('name',nameu,'freq',fu,'tidecon',tidecon,'type',ltype);
   fu=xout;
end;
   
%----------------------------------------------------------------------
function [nameu,fu,ju,namei,fi,jinf,jref]=constituents(minres,constit,...
                                     shallow,infname,infref,centraltime);
% [name,freq,kmpr]=constituents(minres,infname) loads tidal constituent
% table (containing 146 constituents), then picks out only the '
% resolvable' frequencies (i.e. those that are MINRES apart), base on 
% the comparisons in the third column of constituents.dat. Only 
% frequencies in the 'standard' set of 69 frequencies are actually used.
% Also return the indices of constituents to be inferred.

% If we have the mat-file, read it in, otherwise create it and read
% it in!

% R Pawlowicz 9/1/01 
% Version 1.0
%
%    19/1/02 - typo fixed (thanks to  Zhigang Xu)

% Compute frequencies from astronomical considerations.