%% Atm S 547 HW1

%% Plot time series of turbulence data

  rho = 1.21; % [kg/m3]
  Cp = 1004; % [J/kg-K]
  L = 2.5e6; % [J/kg]
  g = 9.8;   % [m/s2]
  U0 = 100;  % [m/s]

% Load time series and normalize time to start at 0

  load('rf18L1.txt');
  t = rf18L1(:,1);
  t = t - t(1);
  u = rf18L1(:,2); 
  v = rf18L1(:,3); 
  w = rf18L1(:,4); 
  T = rf18L1(:,5); 
  q = rf18L1(:,6);
  p = rf18L1(:,7); 

%    ... and plot velocity data

  subplot(2,2,1)
  plot(t,u,'g.',t,v,'r.',t,w,'b.');
  xlabel('Time (s)')
  ylabel('m s^{-1}')
  title('Unfiltered velocity components')
  legend('u','v','w')

%%  2. Spectral analysis of w. Outputs power spectral density Pww as function
%     of frequency f (Hz). See p. 3-10, Matlab Signal Processing Toolbox Manual

  rate = 25;                   % Data rate (Hz)
  nfft = 4096;                 % Samples per data window
  window = hanning(nfft);      % Use Hanning data taper function
  noverlap = nfft/2;           % Overlap between periodogram windows
  [Pww,f] = pwelch(w,window,noverlap,nfft,rate);
  [Puu,f] = pwelch(u,window,noverlap,nfft,rate);

  C = 0.03  % Eyeballed constant in fit Pww  = C*f^-5/3 in inertial range
  subplot(2,2,2)
  p=loglog(f,f.*Pww,f,f.*Puu,'g',f,C*f.^(-2/3),'r');
  xlabel('Frequency [s^{-1}]')
  ylabel('f*(spectral power) [m^2 s^{-2}]')
  title('Power spectra') 
  axis([1e-2,20,3e-3,3])
  hold on
  plot([.3 .3],[3e-3 3],'k--')
  text(0.5,5*C,'Inertial range','Color','Red')
  hold off
  legend(p,'w','u','Cf^{-2/3}')   % Only include lines in p in legend

%%  3. Turbulent dissipation  E(k) = eps^2/3 k^-5/3, k = 2*pi*f/U0
                           
  diss = (2*pi/U0)*(1.25*C)^1.5
  
%%  4. High pass filtering (passband f > 0.05) and std of w 

  cutoff = 0.05*2/rate;         % Cutoff frequency in Nyquist units
  [b, a] = butter(4,cutoff,'high'); 
  whi = filtfilt(b,a,w);
  wstd = std(whi)

%%  5. Autocorrelation timescale tau for w

  whi_lag_corr = xcov(whi,'coeff');
  w_lag_corr = xcov(w,'coeff');
  Tmax = 10;
  lag_range = length(whi) + (0:(Tmax*rate));
  lag = (0:(Tmax*rate))/rate;
  subplot(2,2,4);
  p=plot(lag,w_lag_corr(lag_range),'g',...
         lag,whi_lag_corr(lag_range),'b');
  title('Autocorrelation of w and high-pass w')
  xlabel('Lag (s)')
  hold on
  plot([0 10],[0 0],'k--')  % Zero lagged autocorrelation line
  plot([0 10],exp(-1)*[1 1],'k--') % Intercept for determining tau
  tau_w = interp1(lag,w_lag_corr(lag_range),exp(-1));
  tau_whi = interp1(lag,whi_lag_corr(lag_range),exp(-1));
  text(0.75,0.45,['\tau = ' num2str(round(tau_w,2)) ' s'],'Color','g')
  text(0.75,0.29,['\tau_{hi} = ' num2str(round(tau_whi,2)) ' s'],'Color','b')
  legend(p,'w','whi')

  hold off

%%  6. High pass filtering of other variables

  uhi = filtfilt(b,a,u);
  vhi = filter(b,a,v);
  Thi = filtfilt(b,a,T);
  qhi = filtfilt(b,a,q);

%  Plot high-pass velocity time series

  subplot(2,2,3)
  plot(t,uhi,'g.',t,vhi+5,'r.',t,whi+10,'b.');
  xlabel('Time (s)')
  ylabel('m s^{-1}')
  title('High-pass velocity components')
  legend('u','v+5','w+10')

  ustd = std(uhi)
  vstd = std(vhi)
  Tstd = std(Thi)
  qstd = std(qhi)
  wu = cov(whi,uhi);   
  wu = wu(1,2);            
  wv = cov(whi,vhi);
  wv = wv(1,2);            
  wT = cov(whi,Thi);  
  wT = wT(1,2);            
  wq = cov(whi,qhi);
  wq = wq(1,2);            
  wucorr = wu/(wstd*ustd)
  wvcorr = wv/(wstd*vstd)
  wTcorr = wT/(wstd*Tstd)
  wqcorr = wq/(wstd*qstd)

  uflux = rho*wu
  vflux = rho*wv
  sensflux = rho*Cp*wT % Sensible heat flux
  wq = wq/1000;        % Convert q from g/kg to kg/kg to compute flux
  latflux = rho*L*wq   % Latent heat flux
  buoyflux = g*(wT/mean(T+273) + 0.61*wq)  % Buoyancy flux
