% 1-20-2012

%   TransferFunctionMoore

%       ESPM 228, Advanced Topics in Biomet and Micromet
%       Dennis Baldocchi
%       University of California, Berkeley

%  This matlab code computes the spectral and integrated transfer functions
%  that apply to eddy covariance measurements

%  This program is adapted from the the Basic code of 
% ' CJ Moore. 1986. Boundary Layer Meteorology.37, 17-35
% 

%  The conversion to Matlab includes performing Matrix operation on the spectra
%  We replaced use of the Kaimal spectra with spectra models produced by
%  Detto, Baldocchi and Katul. 2010 Boundary Layer Meteorology, 136: 407-430 
%  that were derived from thousands of runs over a flat peatland pasture, for
%  U, T, H2O, CO2 and CH4 sensors and for stable, near neutral and unstable
%  thermal stratification
% '
% ' **************************************************************
% '
% '
% '
% ' This program corrects for:

%  1) signal attenuation due to a finite path of the sonics, hygrometer and CO2 irgas.
%  2) dynamic frequency response limitations sensors with finite time response, like
%  bead thermistors and the finite residence time air as is passes through in the cell
%  of the closed path irga.
%  3) the effects of sensor separation between w and x. 
%  4) Signal attenuation effects due to digital high and low pass filtering and aliasing
%  
% ' *************************************************************
% '
% '       Z, height above zero plane, s
% '       U, wind speed, m/s
% '       ZL, Monin Obukhov stability parameter
% '       D1, thermistor time constant, s
% '       D2, sonic anemometer time constant, s.
% '           (note in Moore's program this variable is
% '           distance constant (m) of a propellor
% '           anemometer).
% '       D3, time constant of low pass electronic filter, s
% '       D4, time constant of open path CO2/H2O LI 7500
% '       D5, residence time in closed path irga, s
% '       PQ, path length of q sensor, m
% '       PCO, path length of open CO2 sensor, m
% '       PCC, path length of closed CO2 sensor, m
% '       PW, path length of sonic anemo w, m
% '       PU, path length of sonic anemo u, m
% '       XT, separation distance of w and T, m
% '       XQ, separation distance of w and Q, m
% '       XU, separation distance of w and u sensors, m
% '       XCO, separation distance of w and open CO2 sensor
% '       XCC, separation distance of w and closed CO2 sensor
% '       DC, time constant of digital high pass filter, s
% '       NS, sampling frequency, HZ
% '       DCO2, diffusivity of CO2
% '       XL, length of the irga sampling tube
% '       UL, flow velocity thru tube
% '       RAD, radius of tube
% '
% '       covariance      transfer function      correction factor
% '        wt                 G1                     S1
% '        wq                 G2                     S2
% '        uw                 G3                     S3
% '        tt                 G4                     S4
% '        qq                 G5                     S5
% '        uu                 G6                     S6
% '        ww                 G7                     S7
% '        wc open            G8                     S8
% '        cc open            G9                     S9
% '        wc closed          G10                    S10
% '        cc closed          G11                    S11

% '
% '************************************************************
% '
% '       set up measurement parameters
% '

clear all;

        in.ht = 1;              % canopy height, m 
        in.zmeas = 3;           % measurement height, m
        in.NS = 20;             % sampling frequency        
        in.d = .6 * in.ht;      % zero plane displacement, ~ 60% of canopy height
        in.z = in.zmeas - in.d; %  z-d
        in.u = 8;               % wind velocity, m/s
        in.B = 1.667;           %  5/3
        in.ZL = +.3;           % Monin-Obukhov Length, z/L
        in.D1 = .01;           % thermometer time constant, s
        in.D2 = .01;           % sonic time constant: over sampling is done
        in.D3 = 1000;      % sampling time constant low pass filter
        in.D4 = 1/ in.NS;       %  time constant for CO2 and H2O, LI 7500
        in.D5 = .25;            %     1/e residence time of closed path irga
        in.PQ = .15 ;           %    'path size of water vapor sensor, m
        in.PW = .15;            %    'path of w and T, m
        in.PU = .15;            %    'path of u, m
        in.PCO = .15;           %    'path of open CO2/H2O, m
        in.PCC = .15;           %    'path of closed CO2/H2O, m
        in.XT = .01;    %   separation distance between sonic and temperature
        in.XQ = .1;     %      'IRGA is 20 cm from sonic
        in.XCO = .1;    %  separation distance between sonic and open co2
        in.XCC = .01;    % separation distance between sonic and closed co2
        in.XU = .01;     % separation between w and u
        in.DC = 400;    % time constant digital high pass filter
        
        %  setup on IRGA Sampling through a tube
        
        in.DCO2 = 14.7; % '  co2 diffusivity mm^2/s
        in.XL = 3400;   %   ' mm       tube length
        in.UL = 1130;   %   ' mm/s     flow velocity through tube 
        in.RAD = 2.25;  %  ' mm       tube radius
        
        in.bins=49;     % number of bins
        
        % compute time constant of sonic anemometer to be used in Transfer
% '
       
        in.UC = in.D2 ;
% '
% 'square time constant of electronic RC filter
% '
        in.VC = in.D3 ;

        in.TC = in.D1 ;
        
        in.QC = in.D4 ;
        
        in.QCC = in.D5 * in.D5;
        
% Compute natural frequency classes from 10^-4 to 10^1,  to plot and assess spectra
% 
% yielding in.bin classes, gives us spectra from 0.0001 to 10 hz,

in.n(1)=0.0001;   % lowest frequency

%LF= 1.90108;  % log (NS) = log(0.0001) + NS log(a); n = 0.0001 * a^# bins

LF= 10^((log10(in.NS)-log10(0.0001))/in.bins);

% compute the spectral bins for the computations

for i=1:in.bins-1
   
    in.n(i+1)=in.n(i) * LF;
end
       
     
        
%         call SPECTRA, this subroutine computes the coefficients
%         for use in the formulae of Kaimal
        
  %      [yspec]=spectra(in);
  
  % call the spectral functions in Detto et al
  
       [yspec]=Spectra_Detto(in);  
        
%         integration of spectra with Trapezoid Rule.
%         The logarithmic spectrum is integrated between
%         the log frequency range of 10^-4 to nyquist freq (10 Hz)
%         over bins equal log(intervals).

%         do from 1 to 19: n = 0.0001 * a^(I)

% call subroutine in matrix form and then perform sums on the data in the
% loop

%[ytrans]= transfer_matrix(in, yspec);

 [ytrans]= transfer_matrix_Detto(in);




% Consider integrating the co-spectra, power spectra and transfer functions
% with Matlab functions.

% Because we have uneven steps to cover the multiple orders of magnitude
% for the spectra it is better to use the trapezoid integral even though
% Simpson's rule may be preferred.

%z=trapz(in.n, ytrans.G1);

% to reduce numerical errors I have upped the number of bins from 19 to 49

% integrate transfer functions weighted by turbulence spectra

                 
          
    % Compute the integrated correction value from the Transfer functions.
    
    % S_wq = int (T_wq(n) * S_wq(n)) dn/int(S_wq(n))dn
    
    % must apply different spectra depending whether the atmospheric
    % thermal stratification is stable, near neutral or unstable
    
    % The correction factors are compared with the measured power spectra
    % and co spectra for temperature as these have the least correction.
    
   % Stable conditions
   
    if(in.ZL > 0.05)
        indx=1;        
    end    
        
    % near neutral stratification
    
    if(in.ZL > -0.05 && in.ZL < 0.05)
        indx=2;
        
        
      

    end
    
    % unstable thermal stratification
    
    if (in.ZL < -0.05)
        indx=3;
    end
        
        
    
      % plot some power spectra and co-spectra 
figure(1)
          clf;
          SS=reshape(yspec.CoSp(1,2,:),1,in.bins);
          loglog(in.n, in.n .* SS,'+-')
          %semilogx(in.n,ytrans.WT,'+-')
          
          hold all
          
          SS=reshape(yspec.CoSp(2,2,:),1,in.bins);
          loglog(in.n, in.n .* SS,'+-')
          
          hold all
          
          SS=reshape(yspec.CoSp(3,2,:),1,in.bins);
          loglog(in.n, in.n .* SS,'+-')
          title('co-spectrum <wT>')
          legend ('stable', 'near neutral', 'unstable')
          
          
 
 % plot some transfer functions
 
          figure(2)
          clf;
          semilogx(in.n,ytrans.G1,'+-')
          
          hold all
          
          semilogx(in.n,ytrans.G2,'+-')
          
          hold all
          
          semilogx(in.n,ytrans.G8,'+-')
          title('transfer function')
          legend ('<wT>', '<wq>', '<wC>')
          
          
          % transfer function for path averaging
          figure (3)
          clf;
          semilogx(in.n,ytrans.WP,'+-')
          hold all
          semilogx(in.n,ytrans.TP,'+-')
          hold all
          semilogx(in.n,ytrans.qp,'+-')
          hold all
          semilogx(in.n,ytrans.cop,'+-')
          hold all
          semilogx(in.n,ytrans.ccp,'+-')
          title('Transfer function for path averaging')
          legend('w','T','Q','CO2 open', 'CO2 closed')
          
          % transfer function for sensor separation
          figure (4)
          clf;
                    semilogx(in.n,ytrans.TS,'+-')
          hold all
          semilogx(in.n,ytrans.QS,'+-')
          hold all
          semilogx(in.n,ytrans.XCOS,'+-')
          hold all
          semilogx(in.n,ytrans.CCS,'+-')
          title('Transfer function for path averaging')
          legend('T','Q','CO2 open', 'CO2 closed')
          
          
     % transfer function for data acquisition
          figure (5)
          clf;
          semilogx(in.n,ytrans.gx,'+-')
          hold all
          semilogx(in.n,ytrans.ga,'+-')
          hold all
          semilogx(in.n,ytrans.GD,'+-')
          
          title('Transfer function for data acquisition')
          legend('total','aliasing','high pass filter')
          
          
    
S1 = trapz(in.n,ytrans.G1 .* reshape(yspec.CoSp(indx,2,:),1,in.bins)); %wt
S2 = trapz(in.n, ytrans.G2 .* reshape(yspec.CoSp(indx,2,:),1,in.bins)); %wq
S3 = trapz(in.n, ytrans.G3 .* reshape(yspec.CoSp(indx,1,:),1,in.bins)); %wu
S4 = trapz(in.n, ytrans.G4 .* reshape(yspec.Sxx(indx,2,:),1,in.bins)); %tt
S5 = trapz(in.n, ytrans.G5 .* reshape(yspec.Sxx(indx,3,:),1,in.bins)); %qq
S6 = trapz(in.n, ytrans.G6 .* reshape(yspec.Sxx(indx,1,:),1,in.bins)); %uu
%S7 = trapz(in.n, ytrans.G7 .* reshape(yspec.WW);
S8 = trapz(in.n, ytrans.G8 .* reshape(yspec.CoSp(indx,2,:),1,in.bins)); %wc
S9 = trapz(in.n, ytrans.G9 .* reshape(yspec.Sxx(indx,2,:),1,in.bins));  %cc
S10 =trapz(in.n, ytrans.G10 .* reshape(yspec.CoSp(indx,2,:),1,in.bins)); %wc
S11 =trapz(in.n, ytrans.G11 .* reshape(yspec.Sxx(indx,2,:),1,in.bins)); %cc

DS1 = trapz(in.n, reshape(yspec.CoSp(indx,2,:),1,in.bins)); %wt
DS2 = trapz(in.n, reshape(yspec.CoSp(indx,2,:),1,in.bins));  %wq
DS3 = trapz(in.n, reshape(yspec.CoSp(indx,1,:),1,in.bins));  %wu
DS4 = trapz(in.n, reshape(yspec.Sxx(indx,2,:),1,in.bins));  %tt
DS5 = trapz(in.n, reshape(yspec.Sxx(indx,2,:),1,in.bins));  %qq
DS6 = trapz(in.n, reshape(yspec.Sxx(indx,1,:),1,in.bins));  %uu
%DS7 = trapz(in.n, reshape(yspec.WW);
DS8 = trapz(in.n, reshape(yspec.CoSp(indx,2,:),1,in.bins));  %wc
DS9 = trapz(in.n, reshape(yspec.Sxx(indx,2,:),1,in.bins));  %cc
DS10 =trapz(in.n, reshape(yspec.CoSp(indx,2,:),1,in.bins));  %wc
DS11 =trapz(in.n, reshape(yspec.Sxx(indx,2,:),1,in.bins));   %wc

 
%%%%%%%% Divided transfer function, weighted by spectra, by integral of the
%%%%%%%% spectra
        
        S_wt = S1 / DS1;
        S_wq = S2 / DS2;
        S_wu = S3 / DS3;
        S_tt = S4 / DS4;
        S_qq = S5 / DS5;
        S_uu = S6 / DS6;
        %S_ww = S7 / DS7;
        S_wc_open = S8 / DS8;
        S_cc_open = S9 / DS9;
        S_wc_closed = S10 / DS10;
        S_cc_closed = S11 / DS11;

        
    
    
    
        

    
    
        

        
        
             
        
          