function profile=vad_analysis(lidar)
%************************************************************************
% MATLAB m-file for calculating Fourier coefficients of radial Velocity
%    function that can be used to obtain the kinematic properties of the
%    local wind field.
%  VAD_analysis_on_sweeps
%				
%				Melvin Felton
%				Army Research Laboratory, Adelphi, MD
%				mfelton@arl.army.mil
%				Nov 2003
%
%	INPUT:	rws: radial wind as a na x nr array
%           azimuth: azimuth of scan na x 1 array
%           elevation: elevation of scan na x 1 array
%           range: range position 1 x nr array;
%
%	OUTPUT:	Fourier coefficients of radial velocity function, VAD for
%	        radial velocity and best fit curve, VAD parameters (velocity
%	        and azimuth of min), vertical profiles of wind direction,
%	        speed, and deformation.
%	        
%*************************************************************************
rws=lidar.RWS;
azimuth=lidar.azimuth;
elevation=lidar.elevation;
range=sqrt(lidar.xd.^2+lidar.yd.^2+lidar.zd.^2);
snr=lidar.SNR
[na,nr]=size(rws);
f1=ones(na,nr);
f2=sind(azimuth);
f3=cosd(azimuth);
f4=sind(2*azimuth);
f5=cosd(2*azimuth);
f6=sind(3*azimuth);
f7=cosd(3*azimuth);
 
disp(' ')
disp('Performing VAD analysis...')
disp(' ')
disp(' ')

%initailize arrays
b = zeros(7,nr);                              %fourier analysis
a = zeros(7,nr);                              %fourier coefficients    
A=  zeros(7,7,nr);                            %fourier analysis
s_squared= zeros(1, nr);                      %estimated variance about regression
V= zeros(na,nr);                         %radial velocity
Azimuths= V;                                  %azimuths
Elevations= V;                                %elevations
Hor_ranges= V;                                %horizontal ranges
V_hat= V;                                     %best fit curve of radial velocity

%calculate Fourier coefficients b and A based upon method outlined in Matejka et
%  al., 1991.  analysis excludes NaN data points. 
for i=1:nr
    m = find(~isnan(rws(:,i)));                
    D = rws(m,i);
    azim= azimuth(m,i);
    elev=elevation(m,i);
 %  V(i,m)= D(m,i);
 %  Azimuths(i,m)= azimuth(m,i);
 %  Elevations(i,m)= elevation(m,i);
 %  Hor_ranges(i,m)= cos(Elevations(i,m)*pi/180).*range(i,m);
    %=[ones(size(azim)), sind(azim), cosd(azim), sind(2*azim), cosd(2*azim), sind(3*azim), cosd(3*azim)]';
    %=T.*D;
    for j=1:length(D)

        b(1,i)= b(1,i)+ D(j);
        b(2,i)= b(2,i)+ (sin(azim(j)*pi/180) * D(j));
        b(3,i)= b(3,i)+ (cos(azim(j)*pi/180) * D(j));
        b(4,i)= b(4,i)+ (sin(2*azim(j)*pi/180) * D(j));
        b(5,i)= b(5,i)+ (cos(2*azim(j)*pi/180) * D(j));
        b(6,i)= b(6,i)+ (sin(3*azim(j)*pi/180) * D(j));
        b(7,i)= b(7,i)+ (cos(3*azim(j)*pi/180) * D(j));
            
        x= azim(j)*pi/180;
        functions= [[1] [sin(x)] [cos(x)] [sin(2*x)] [cos(2*x)] [sin(3*x)] [cos(3*x)]];
        A(:,:,i)= A(:,:,i) + (functions' * functions);
        
    end
end


%solve for Fourier coefficients a
for i=1:nr
    a(:,i) = A(:,:,i)\b(:,i);
end
    

%calculate best-fit points and calculate estimated variance about
%  the regression 
for i=1:nr
    m = find(~isnan(rws(:,i)));  
    D = rws(m,i); azim= azimuth(m,i);
    for j=1:length(D)
        x= azim(j)*pi/180;
        functions= [[1] [sin(x)] [cos(x)] [sin(2*x)] [cos(2*x)] [sin(3*x)] [cos(3*x)]];
        V_hat(i,j)= sum(a(:,i)'.*functions);
    end
    
    s_squared(i)= (1/(length(D)-7))*sum((D(m)- V_hat(i,m)').^2);
    
%calculate variance
    for j= 1:7
        var(j,i)= A(j,j,i)\s_squared(i);
    end
    
end


for i=1:num_gates
%eliminate dummy values
    for j=1:length(D_sweep)
       if V(i,j)==0
         V_hat(i,j)=0;
         Hor_ranges(i,j)=0;
         Elevations(i,j)=0;
       end
%eliminate ground clutter
       if (abs(V(i,j))< 1)
            V(i,j)=0;
            Hor_ranges(i,j)=0;
            Elevations(i,j)=0;
        end
%eliminate outliers
         if (sqrt(s_squared(i))>=1.5) & (V(i,j)~=0)
            if (abs((V(i,j))-(V_hat(i,j)))>10) | (abs((V(i,j))-(V_hat(i,j)))> 2*(sqrt(s_squared(i))))
                V(i,j)=0;
                Hor_ranges(i,j)=0;
                Elevations(i,j)=0;
            end
        end  
    end
end
    

%reiterate regression
t= find(sqrt(s_squared)>= 1.5);
for r=1:2
    for i=1:length(t)
        for j=1:length(V(t(i),:))
        
            b(1,t(i))= b(1,t(i))+ V(t(i),j);
            b(2,t(i))= b(2,t(i))+ (sin(Azimuths(t(i),j)*pi/180) * V(t(i),j));
            b(3,t(i))= b(3,t(i))+ (cos(Azimuths(t(i),j)*pi/180) * V(t(i),j));
            b(4,t(i))= b(4,t(i))+ (sin(2*Azimuths(t(i),j)*pi/180) * V(t(i),j));
            b(5,t(i))= b(5,t(i))+ (cos(2*Azimuths(t(i),j)*pi/180) * V(t(i),j));
            b(6,t(i))= b(6,t(i))+ (sin(3*Azimuths(t(i),j)*pi/180) * V(t(i),j));
            b(7,t(i))= b(7,t(i))+ (cos(3*Azimuths(t(i),j)*pi/180) * V(t(i),j));
            
            x= Azimuths(t(i),j)*pi/180;
            functions= [[1] [sin(x)] [cos(x)] [sin(2*x)] [cos(2*x)] [sin(3*x)] [cos(3*x)]];
            A(:,:,t(i))= A(:,:,t(i)) + (functions' * functions);
        end
    
        a(:,t(i)) = A(:,:,t(i))\b(:,t(i));
    
        for j=1:length(V(t(i),:))
            x= Azimuths(t(i),j)*pi/180;
            functions= [[1] [sin(x)] [cos(x)] [sin(2*x)] [cos(2*x)] [sin(3*x)] [cos(3*x)]];
            V_hat(t(i),j)= sum(a(:,t(i))'.*functions);
        end
    
        o= find(V(t(i),:)==0);
        n= find(V(t(i),:)~=0);
        s_squared(t(i))= (1/(length(V(t(i),:))-7-(length(V(t(i),:))-length(o))))*sum((V(i,n)- V_hat(i,n)).^2);
         
    
        for j=1:length(V(t(i),:))
            if (abs((V(t(i),j))-(V_hat(t(i),j)))>10) | (abs((V(t(i),j))-(V_hat(t(i),j)))> 2*(sqrt(s_squared(t(i)))))
                V(t(i),j)=0;
                Hor_ranges(t(i),j)=0;
                Elevations(t(i),j)=0;
            end
        end
                   
    end
    t= find(sqrt(s_squared)>= 1.5);
end 



%create VAD plots for radial velocity and best-fit curve
VAD_plot= input('Would you like to create a VAD plot?  Type 1 for yes or 2 for no.\n');
disp('------------------------------------------------------')
while (VAD_plot < 1) | (VAD_plot > 2)
    VAD_plot= input('Invalid entry.  Please type 1 for yes or 2 for no.\n');
    disp('------------------------------------------------------')
end

if (VAD_plot == 1)
    range_gate= input('Enter a range gate between 1 and 80.\n');
    disp('------------------------------------------------------')
    while (range_gate < 1) | (range_gate > num_gates)
        range_gate= input('Invalid entry.  Please enter a number between 1 and 80.\n');
        disp('------------------------------------------------------')
    end
    
    figure;

    o= find(V(range_gate,:)~=0);
    %subplot(1,2,1)
    plot(Azimuths(range_gate,o), V(range_gate,o), 'k.')
    xlabel('Azimuth (degrees)')
    ylabel('Radial Velocity (m/s)')
    title('VAD')
    %subplot(1,2,2)
    hold on
    plot(Azimuths(range_gate,:), V_hat(range_gate,:), 'r.')
    %xlabel('Azimuth (degrees)')
    %ylabel('Radial Velocity (m/s)')
    %title('Best fit VAD')
    
    %figure;
    %plot(Azimuths(range_gate,o), abs(V(range_gate,o)-V_hat(range_gate,o)), 'k.')
    
end




%create vertical profiles of hieght vs. wind speed, deformation, and direction.
%  the shear is calculated within a coordinate system in which the x axis is
%  drawn parallel to the wind direction.
for i=1:num_gates
    
    n= find(V(i,:)~=0);
    
[C1,I1] = min(V_hat(i,:));
[C2,I2] = max(V_hat(i,:));

wind_field(i, 1)= abs(C1)+abs(C2);
%wind_field(i, 1)= -(a(2,i)^2+a(3,i)^2)^1/2/cos(mean(Elevations(i,n)*pi/180));    %magnitude of wind speed
wind_field(i, 2)= (C2-C1)/2*cos(Data{2,1}.El(i)*pi/180);   %horizontal wind speed
wind_field(i, 3)= (C2+C1)/2*sin(Data{2,1}.El(i)*pi/180);   %vertical wind speed
wind_field(i, 4)= Azimuths(i, I1);             %azimuth of min
wind_field(i, 5)= Azimuths(i, I2);             %azimuth of max
wind_field(i, 6)= Header{2,1}.Ranges(i);       %range to target
wind_field(i, 7)= (sin(El_sweep(1)*pi/180)).*Header{2,1}.Ranges(i)/1000;   %elevation of target

if (a(2,i) < 0) 
    theta(i)= (pi/2)-atan((a(3,i)/a(2,i))*pi/180);
end

if (a(2,i) >= 0) 
    theta(i)= (3*pi/2)-atan((a(3,i)/a(2,i))*pi/180);
end

if a(4,i)<0
    gamma(i)= (pi/4)-1/2*atan(((a(5,i)*cos(2*theta(i)*pi/180)+a(4,i)*sin(2*theta(i)*pi/180))/(a(4,i)*cos(2*theta(i)*pi/180)+a(5,i)*sin(2*theta(i)*pi/180)))*pi/180);
end

if a(4,i)>=0
    gamma(i)= 3*pi/4-1/2*atan(((a(5,i)*cos(2*theta(i)*pi/180)+a(4,i)*sin(2*theta(i)*pi/180))/(a(4,i)*cos(2*theta(i)*pi/180)+a(5,i)*sin(2*theta(i)*pi/180)))*pi/180);
end

deform(i)= -2*((a(5,i)*cos(2*theta(i)*pi/180)+a(4,i)*sin(2*theta(i)*pi/180))^2+(a(4,i)*cos(2*theta(i)*pi/180)+a(5,i)*sin(2*theta(i)*pi/180))^2)^1/2/mean(Hor_ranges(i,n))*cos(mean(Elevations(i,n))*pi/180);

E= 6.371*10^6;
diverg(i)= 2*a(1,i)/mean(Hor_ranges(i,n))*cos((mean(Elevations(i,n))*pi/180)+(1/E + (-4*10^-8))*mean(Hor_ranges(i,n)));

end
height=wind_field(:,7);
windspeed=wind_field(:,1);
winddirection=wind_field(:,4);
save 'okcprofile.dat' height windspeed winddirection -ascii
%plot vertical profiles
figure; 
subplot(1,4,1)
plot(wind_field(:, 4), wind_field(:,7), 'k-')
xlabel('direction (degrees)')
ylabel('height (km)')
title('Vertical profile of wind direction')
subplot(1,4,2)
plot(wind_field(:,1), wind_field(:,7), 'k-')
%hold on
%plot(deform, wind_field(:,7), 'r--')
xlabel('speed (m/s)')
ylabel('height (km)')
title('Vertical profile of wind speed') 
subplot(1,4,3)
plot(deform, wind_field(:,7), 'r-')
xlabel('resultant deformation (1/s)')
ylabel('height (km)')
title('Vertical Profile of resultant deformation')
%subplot(1,5,4)
%plot(wind_field(:,2), wind_field(:,7), 'k-')
%xlabel('horizontal wind speed (m/s)')
%ylabel('height (km)')
%title('Vertical profile of horizontal wind speed')
subplot(1,4,4)
plot(diverg(2:num_gates), wind_field(2:num_gates,7), 'k-')
xlabel('horizontal divergence (1/s)')
ylabel('height (km)')
title('Vertical profile of horizontal divergence')
      
