function [ out ] = transfer_matrix_Detto( in )
%
%        constants
% 
%         omega, 2 pi n
%         '
        w = 6.283 .* in.n;  % convert frequency to radians
%         '
%         '  wavenumber, nu = n/u
%         '
        NU = in.n ./ in.u;
%         '
%         ' 2 pi n/u
%         '
        W_U = w ./ in.u;
%         '
%         ' omega squared, (2 pi n)^2
%         '
        W2 = w .* w;
        
        % ' coefficient power, 5/3
% '
    B = 1.667;
    
    
%************************************************************    
% '
% '       data acquisition sampling transfer function
% '
% '       T(n) = 0 if n > ns/2

%         else
% 
%        T(n) = 1 +(n/(ns-n)^b if n <= ns/2
% 
%       
%         ' Tranfer function for sampling and accounting for
%         ' the effects of aliasing.
%         '
%         ' if f < ns/2
%         '
%         '   Ga = 1 +(n/(n-ns))^3
%         '
        ga=zeros(1,in.bins);
        
        X = in.n ./ (in.NS - in.n);
        ga = 1 + X .* X .* X;
        
        ga(in.bins)=0;  % value for f = NS
%         '
%         '
%         Gain function for averaging..it is effectively a high pass filter
%         and removes low frequency trends

%         we are approximating this file with a digital filter filter
%         This is the filter we used in the early days of eddy covariance
%         when we did not have suitable storage space to  store raw data.
%         Instead we computed covariances on the fly by computing running
%         means with a digital filter
%
%         the running mean <x> is computed as:

%         <x>_i = alpha <x>_i-1 + (1-alpha) x_i

%         alpha = exp(-dt/Tau)...dt is the time interval for sampling, 1/NS
%         and Tau is the filter time constant
%         '
%         '   Gd(n) = 2 pi n tau/(1 + (2 pi n tau)^2/alpha)^.5
%         '
        GD = zeros(1,in.bins);
        GD(:)=1.;
%         '
%      
%         '
        X = w .* in.DC;
        ALPHA = exp(-1 / (in.NS * in.DC));
        
        X2 = X .* X; 
        GD = X ./ ((1 + X2) ./ ALPHA) .^0.5;
        
         
        
		
%         '
%         'NO DC FILTER APPLIED for electronic noise
 
%         '
%         ' Electronic filter transform, Butterworth Filter
%         '
%         '  (1 + (2 pi n tau)^2)^.5
%         '
            GV=zeros(1,in.bins);
            GV(:) = 1;

%         GV = W2 * in.VC;
%         GV = 1 + GV * GV;
%         '
%         ' no filter is applied
%         '
            
%         '
%         ' gx is the Transfer function of data acquisition system
%         '
%         '  Ga is the sampling rate and aliasing transform
%         '  Gd is the digital high pass filter transform,
%         '       as applied for mean removal.
%         '  Gv is the electronic filter transfor
%            GD is applied twice because each signal is filtered

        gx = ga .* GD .* GD .* GV .* GV;

        
%***************************************************************
%         ' first order gain for dynamic frequency
%         ' response of temperature probe with thermal inertia
%         '
%         '   G(n) = TR(n)^-.5
%         '
%         '  TR(n) = 1 + (2 pi n tau)^2
%         '
        Gtr = (1 + W2 .* in.TC .* in.TC) .^ -.5;
        
        
%         ' first order gain for dynamic frequency
%         ' response of residence time in open path irga
%         '
%         '   G(n) = TR(n)^-.5
%         '
%         '  TR(n) = 1 + (2 pi n tau)^2
%         '
        Gcr_op = (1 + W2 .* in.QC .* in.QC) .^ -.5;   % note need dot carrot
%         '
%         '
%         ' first order gain for dynamic frequency
%         ' response of residence time in closed path irga
%         '
%         '   G(n) = TR(n)^-.5
%         '
%         '  TR(n) = 1 + (2 pi n tau)^2
%         '
        Gcr = (1 + W2 .* in.QCC .* in.QCC) .^ -.5;   % note need dot carrot
%         '
%         'First order gain for dynamic frequency response of
%         'anemometer
%         '
        Gur = (1 + W2 .* in.UC .* in.UC) .^ -.5;
        
        
%********************************************************************        
%         '
%         ' Transform function for line averaging scalars
%         '
%         '  Tp(n) = 3+exp(-x) - 4(1-exp(-x)/x)
%         '
%         '  x =2 pi f, f = n p/u

       
%         '
%         ' for the u component of the sonic
%         '
        up=zeros(1,in.bins);
        up(:) = 1;
        
        X = W_U .* in.PU;  % 2 pi/U * pathlength
        
              
        Y = exp(-X);
        upp = 1 + (Y - 3 .* (1 - Y) ./ X) ./ 2;
        up = 4 .* upp ./ X;
     
        
%         '
%         ' for the w component of the sonic
%         '
        WP=zeros(1,in.bins);
        WP(:) = 1;
        
        X = W_U .* in.PW;  % 2 pi n/U * path length
        
        Y = exp(-X);
        WPP = 1 + (Y - 3 .* (1 - Y) ./ X) ./ 2;
        WP = 4 .* WPP ./ X;
        
                
      
%         '
%         'for the open path CO2 sensor
%         '
        cop=zeros(1,in.bins);
        cop(:) = 1;
        
        X = W_U .* in.PCO;
        
        
        Y = exp(-X);
        copp = 3 + Y - 4 .* (1 - Y) ./ X;
        cop = copp ./ X;
        
            
        
       
%         '
%         'for the closed path CO2 sensor
%         '
        ccp=zeros(1,in.bins);
        ccp(:) = 1;
        
        X = W_U .* in.PCC;
        
        Y = exp(-X);
        ccpp = 3 + Y - 4 .* (1 - Y) ./ X;
        ccp = ccpp ./ X;
        
        %         '
%         ' for the krypton hygrometer
%         '
        qp=zeros(1,in.bins);
        qp(:) = 1;
        X = W_U .* in.PQ;
        
             
         Y = exp(-X);
         qpp = 3 + Y - 4 .* (1 - Y) ./ X;
         qp = qpp ./ X;
         
         %         '
%         ' for the sonic anemometer temperature
%         '
        TP=zeros(1,in.bins);
        TP(:) = 1;
        X = W_U .* in.PW;
        
             
         Y = exp(-X);
         TPP = 3 + Y - 4 .* (1 - Y) ./ X;
         TP = TPP ./ X;
           
%***************************************************************       
%         '
%         ' transfer function for lateral and longitudinal
%         ' separation between a sensor and w.
%         '
%         '  for the temperature sensor
%         '
        TS=zeros(1,in.bins);
        TS(:) = 1;
        
        X = NU .* in.XT;
        
        TS = exp(-9.9 .* X .^ 1.5);
        
       
       
%         '
%         '  for the open path IRGA humidity sensor
%         '
        QS=zeros(1,in.bins);
        QS(:) = 1;
        
        X = NU .* in.XQ;
        
        QS = exp(-9.9 .* X .^ 1.5);
        
             
        
%         '
%         '  for the open path CO2 sensor
%         '
        XCOS=zeros(1,in.bins);
        XCOS(:) = 1;
        
        X = NU .* in.XCO;
        XCOS = exp(-9.9 .* X .^ 1.5);
        
%         '
%         '  for the closed path CO2 sensor
%         '
        CCS=zeros(1,in.bins);
        CCS(:) = 1;
        
        X = NU .* in.XCC;
        CCS = exp(-9.9 .* X .^ 1.5);
        
                
        US=zeros(1,in.bins);
        US(:)=1;
        
        
% '*********************************************************************
% '       Compute fluctuation dampening via sampling thru a tube
% '
        CCD = exp(-.021 .* w .* w .* in.RAD .* in.RAD .* in.XL ./ (in.DCO2 .* in.UL));
%         '
% **************************************************************
%         ' compute the gains...double check when they should be squared, 
%           multiplied or the square root applied..is confusin

%        debug...adding time response factor for the LI-7500 cr_op
%         

%       Gain_wc = Gsampling * Tsensor separation * Tpath_w * Tpath_c
%       * Tresponse_w * Tresponse_c

%       Gain_cc = Gsampling * Tpath_c * Tpath_c
%       * Tresponse_c * Tresponse_c

%       if closed path is used then T_tube is needed, too

%         out.G1 = gx .* TS .* sqrt(WP .* TP) .* Gtr;  % wt
%         out.G2 = gx .* QS .* sqrt(WP .* qp) .* Gcr_op;  % wq
%         out.G3 = gx .* US .* sqrt(WP .* up) .* Gur; % uw
%         out.G4 = gx .* Gtr .* Gtr; %tt
%         out.G5 = gx .* qp .* Gcr_op .* Gcr_op; % qq
%         out.G6 = gx .* up .* Gur .* Gur;  % uu
%         out.G7 = gx .* WP .* Gur .* Gur; %ww
%         out.G8 = gx .* XCOS .* sqrt(WP .* cop)  .* Gcr_op; % wc or wq open
%         out.G9 = gx .* cop .* Gcr_op .* Gcr_op; % cc or qq open
%         out.G10 = gx .* CCS .* sqrt(WP .* ccp .* CCD) .* Gcr ; % wc closed
%         out.G11 = gx .* ccp .* Gcr .* Gcr.* CCD;  % % cc closed
%         

%  The sensors dont have a finite time response so I dropped Gxr

        out.G1 = gx .* TS .* sqrt(WP .* TP) ;  % wt
        out.G2 = gx .* QS .* sqrt(WP .* qp);  % wq
        out.G3 = gx .* US .* sqrt(WP .* up) ; % uw
        out.G4 = gx .* WP; %tt
        out.G5 = gx .* qp ; % qq
        out.G6 = gx .* up ;  % uu
        out.G7 = gx .* WP ; %ww
        out.G8 = gx .* XCOS .* sqrt(WP .* cop)  ; % wc or wq open
        out.G9 = gx .* cop ; % cc or qq open
        out.G10 = gx .* CCS .* sqrt(WP .* ccp .* CCD)  ; % wc closed
        out.G11 = gx .* ccp .* CCD;  % % cc closed
        
        out.WP=WP;
        out.TP=TP;
        out.qp=qp;
        out.cop=cop;
        out.ccp=ccp;
        out.TS=TS;
        out.QS=QS;
        out.XCOS=XCOS;
        out.CCS=CCS;
        
        out.gx=gx;
        out.ga=ga;
        out.GD=GD;
       
        

end

