function [fsupsd,fsucor,fsufreq]=... 
   fsupsdgenerate(samplefreq,N,z,Xmean,sigmasquared,errorcheck)
%
% function [fsupsd,fsucor,fsufreq]= ...
%    fsupsdgenerate(samplefreq,N,z,Xmean,sigmasquared,errorcheck)
%
%  This function generates a fft-ordered spectrum associated
%  with the FSU model.  The Flat Smooth Uniform model is used to describe
%  the spectrum of boundary layer turbulence.
%
%  inputs:
%
%  samplefreq  Equivalent sampling frequency (in Hertz) in the returned fsu.
%
%  N  length of desired FSU  (must be even!!!!!)
%
%  z  height
%
%  Xmean velocity in (units of z)/sec
%
%  sigmasquared  desired E(X(t)X(t))
%
%  errorcheck  Set to 1 to check for errors
%
%  outputs:
%
%  fsupsd  power-spectral density from FSU with deltaf=samplefreq/N
%
%  fsupsd(1)=power at 0 Hertz
%  fsupsd(2)=power at deltaf
%  fsupsd(3)=power at 2*deltaf
%  etc. with the normal conjugate-structure for the last half of
%  fsu.
%
%  fsucor fsucorrelation function with deltat=1/samplefreq
%
%  fsucor(1)=E(X(t)X(t))
%  fsucor(2)=E(X(t)X(t+deltat))
%  fsucor(3)=E(X(t)X(t+2*deltat)) 
%  etc. with the normal periodic-fft structure for the last half.
%
%  To re-order fsucor, use
%
%  [(fsucor((N/2+2):N)') fsucor(1:(N/2+1))']
%
%  This puts time=0 an location N/2, with negative times to the left
%  and positive to the right.
%
%  fsufreq  frequencies at which fsupsd are calculated.  Note that only the 
%  positive frequencies are used.
%

if errorcheck==1
   if ~(floor(N/2)*2==N);
      error('Error in fsupsdgenerate:  N should be even'); end
end % end of errorcheck


deltaf=samplefreq/N;
freq=(0:(N/2))*deltaf;
freq(1)=0.000001*deltaf;
fsufreq=freq;

fsupsd=zeros([N 1]);

%  Create positive-frequency info

fsupsd(1:(N/2+1))=20.53*z/Xmean*freq./...
   (1+475.1*(z/Xmean*freq).^(5/3) )*sigmasquared./freq/2;

%  Now negative-frequency stuff

fsupsd((N/2+2):N)=conj(flipud(fsupsd(2:(N/2))));

fsucor=zeros([N 1]);
fsucor=real(ifft(fsupsd))*samplefreq;

