function prof=vvp_regression(vad)

[na nr]=size(vad.RWS);
ind=1;
oops=(na-1)*nr;
xd=reshape(vad.xd(1:na-1,:),oops,1);
yd=reshape(vad.yd(1:na-1,:),oops,1);
zd=reshape(vad.zd(1:na-1,:),oops,1);
vr=reshape(vad.RWS(1:na-1,:),oops,1);
beta=reshape(vad.Beta(1:na-1,:),oops,1);
snr=reshape(vad.SNR(1:na-1,:),oops,1);
el=reshape(vad.elevation(1:na-1,:),oops,1);
az=reshape(vad.azimuth(1:na-1,:),oops,1);
drange=double(vad.ranges(2)-vad.ranges(1));
prof.heights=[0.5:3:50];
dheight=drange;
xpos=0.;
ypos=0;
rho=1000;
nheights=length(prof.heights);
clear vad
for j=1:length(prof.heights)
    zpp=prof.heights(j);
    ind=find(sqrt((xpos-xd).^2+(ypos-yd).^2)<rho & (abs(zd-prof.heights(j))<=dheight/2) & (~isnan(vr)) );
    np=6;
    na=length(ind);
    if (na>np)
        clear Cook H resid X vel ss DFFITS
        prof.beta(j)=mean(beta(ind));
        prof.avgsnr(j)=mean(snr(ind));
        prof.ninda(j)=length(ind);
        xdt=xd(ind);
        ydt=yd(ind);
        zdt=zd(ind);
        elt=el(ind);
        azt=az(ind);
        vrt=vr(ind);
        %X=pred(xdt,ydt,zdt,xpos,ypos,zpp,elt*pi/180,azt*pi/180,np);
        X=VVP_basis_dal(xdt,ydt,zdt,xpos,ypos,zpp,elt*pi/180,azt*pi/180);
        [nn,pp]=size(X);
        if nn==np
            X=X';
        end
        vel=pinv(X'*X)*X'*vrt;
        resid=(vrt-X*vel);
        ss=sqrt(var(resid));
        
        prof.vela(j,:)=vel;
        prof.turba(j)=ss;

        H=X*pinv(X'*X)*X';
        for m=1:length(resid)
            resid(m)=resid(m)/ss/sqrt(1-H(m,m));
        end
        C=pinv(X'*X)*(resid'*resid)/(na-np);
        [nn mm]=size(H);
        for n=1:nn
            Cook(n)=(resid(n).^2/np)*(H(n,n)/(1-H(n,n)));
            sn=((nn-pp)*ss - resid(n)^2/(1-H(n,n)))/(nn-pp-1);
            DFFITS(n)=sqrt(H(n,n)/(1-H(n,n)))*(resid(n)/sn/sqrt(1-H(n,n)));

        end

        %calculate scaled Z=XS^-1;
        inda=find(Cook<4/nn & DFFITS<2*sqrt(pp/nn));
        clear Cook H resid X vel ss DFFITS
        if length(inda)>np
            clear xdtt ydtt zdtt eltt aztt vrtt
            prof.nindb(j)=length(inda);
            xdtt=xdt(inda);
            ydtt=ydt(inda);
            zdtt=zdt(inda);
            eltt=elt(inda);
            aztt=azt(inda);
            vrtt=vrt(inda);

            X=VVP_basis_dal(xdtt,ydtt,zdtt,0,0,zpp,eltt*pi/180,aztt*pi/180);
            [nn,pp]=size(X);
            if nn==np
                X=X';
            end
            vel=pinv(X'*X)*X'*vrtt;
            resid=(vrtt-X*vel);
            ss=sqrt(var(resid));

            prof.velb(j,:)=vel;
            prof.turbb(j)=ss;

            H=X*pinv(X'*X)*X';
            for m=1:length(resid)
                resid(m)=resid(m)/ss/sqrt(1-H(m,m));
            end
            C=pinv(X'*X)*(resid'*resid)/(na-np);
            [nn mm]=size(H);
            clear Cook DFFITS
            for n=1:nn
                Cook(n)=(resid(n).^2/np)*(H(n,n)/(1-H(n,n)));
                sn=((nn-pp)*ss - resid(n)^2/(1-H(n,n)))/(nn-pp-1);
                DFFITS(n)=sqrt(H(n,n)/(1-H(n,n)))*(resid(n)/sn/sqrt(1-H(n,n)));
            end
            %calculate scaled Z=XS^-1;
            indb=find(Cook<4/length(Cook) & DFFITS<2*sqrt(pp/nn));
            clear Cook H resid X vel ss DFFITS
            if length(indb)>np
                prof.nind(j)=length(indb);
                xdttt=xdtt(indb);
                ydttt=ydtt(indb);
                zdttt=zdtt(indb);
                elttt=eltt(indb);
                azttt=aztt(indb);
                vrttt=vrtt(indb);

                X=VVP_basis_dal(xdttt,ydttt,zdttt,0,0,zpp,elttt*pi/180,azttt*pi/180);
                [nn,pp]=size(X);
                if nn==np
                    X=X';
                end
                vel=pinv(X'*X)*X'*vrttt;
                resid=(vrttt-X*vel);
                ss=sqrt(var(resid));

                prof.vel(j,:)=vel;
                prof.turb(j)=ss;
                
                H=X*pinv(X'*X)*X';
                for m=1:length(resid)
                    resid(m)=resid(m)/ss/sqrt(1-H(m,m));
                end
                C=pinv(X'*X)*(resid'*resid)/(na-np);
                [nn mm]=size(H);
                for n=1:nn
                    Cook(n)=(resid(n).^2/np)*(H(n,n)/(1-H(n,n)));
                    sn=((nn-pp)*ss - resid(n)^2/(1-H(n,n)))/(nn-pp-1);
                    DFFITS(n)=sqrt(H(n,n)/(1-H(n,n)))*(resid(n)/sn/sqrt(1-H(n,n)));
                end
            else
                prof.vel(j,:)=prof.velb(j,:);
                prof.turb(j)=prof.turbb(j);
                 clear Cook H resid X vel ss DFFITS
            end
        else
            prof.vel(j,:)=prof.vela(j,:);
            prof.turb(j)=prof.turba(j);
             clear Cook H resid X vel ss DFFITS

        end
    else
        prof.vel(j,:)=ones(1,3)*NaN;
        prof.turb(j)=NaN;
        prof.beta(j)=NaN;
        prof.avgsnr(j)=NaN;
        prof.ninda(j)=NaN;
         clear Cook H resid X vel ss DFFITS

    end

end
ind=find(abs(prof.vel)>30);
prof.vel(ind)=NaN;
theta=atand(abs(prof.vel(:,1))./abs(prof.vel(:,2)));
for j=1:length(prof.heights)
if (prof.vel(j,1)>0 & prof.vel(j,2)>0)
theta(j)=theta(j);
elseif (prof.vel(j,1)>0 & prof.vel(j,2)<0)
theta(j)=180-theta(j);
elseif (prof.vel(j,1)<0 & prof.vel(j,2)<0)
theta(j)=theta(j)+180;
elseif (prof.vel(j,1)<0 & prof.vel(j,2)>0)
theta(j)=360-theta(j);
end
end

prof.theta=theta;
prof.u=prof.vel(:,1);
prof.v=prof.vel(:,2);
prof.w=prof.vel(:,3);
prof.Vh=sqrt(prof.u.^2+prof.v.^2);
     
        