function prof=vvp_regression(vad)

[na nr]=size(vad.RWS);
ind=1;
oops=na*nr;
xd=reshape(vad.xd,oops,1);
yd=reshape(vad.yd,oops,1);
zd=reshape(vad.zd,oops,1);
vr=reshape(vad.RWS,oops,1);
beta=reshape(vad.Beta,oops,1);
snr=reshape(vad.SNR,oops,1);
el=reshape(vad.elevation,oops,1);
az=reshape(vad.azimuth,oops,1);
drange=double(vad.ranges(2)-vad.ranges(1));
prof.heights=[15:1:1000];
dheight=16*drange;
nheights=length(prof.heights);
clear vad
for j=1:length(prof.heights)
    zpp=prof.heights(j);
    ind=find((abs(zd-prof.heights(j))<=dheight/2) & (~isnan(vr)) );
    np=5;
    na=length(ind);
    if (na>np+10)
        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,0,0,zpp,elt*pi/180,azt*pi/180,np);
        [nn,pp]=size(X);
        if nn==np
            X=X';
        end
        vel=pinv(X'*X)*X'*vrt;
        resid=(vrt-X*vel);
        ss=sqrt(var(resid));
        vela(j,:)=vel;
        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/length(Cook) & DFFITS<2*sqrt(pp/nn));
        clear Cook H resid X vel ss DFFITS
        if length(inda)>np+5
            prof.nindb(j)=length(inda);
            xdtt=xdt(inda);
            ydtt=ydt(inda);
            zdtt=zdt(inda);
            eltt=elt(inda);
            aztt=azt(inda);
            vrtt=vrt(inda);

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

            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
            
                velb(:)=vel;
                turbb=ss;
            %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=pred(xdttt,ydttt,zdttt,0,0,zpp,elttt*pi/180,azttt*pi/180,np);
                [nn,pp]=size(X);
                if nn==np
                    X=X';
                end
                vel=pinv(X'*X)*X'*vrttt;
                resid=(vrttt-X*vel);
                ss=sqrt(var(resid));
                velc(j,:)=vel;
                turbc(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
                velc(:)=vel;
                turbc=ss;
                %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(indc)>np
                    prof.nind(j)=length(indc);
                    xdtttt=xdttt(indc);
                    ydtttt=ydttt(indc);
                    zdtttt=zdttt(indc);
                    eltttt=elttt(indc);
                    aztttt=azttt(indc);
                    vrtttt=vrttt(indc);

                    X=pred(xdtttt,ydtttt,zdtttt,0,0,zpp,eltttt*pi/180,aztttt*pi/180,np);
                    [nn,pp]=size(X);
                    if nn==np
                        X=X';
                    end
                    vel=pinv(X'*X)*X'*vrtttt;
                    resid=(vrtttt-X*vel);
                    ss=sqrt(var(resid));
                    prof.vel(j,:)=vel;
                    prof.turb(j)=ss;
                else
                    prof.vel(j,:)=velc(j,:);
                    prof.turb(j)=turbc(j);
                end
            else
                prof.vel(j,:)=velb(j,:);
                prof.turb(j)=turbb(j);
            end
        else
            prof.vel(j,:)=prof.vela(j,:);
            prof.turb(j)=prof.turba(j);

        end
    else
        prof.vel(j,:)=ones(1,np)*NaN;
        prof.turb(j)=NaN;
        prof.beta(j)=NaN;
        prof.avgsnr(j)=NaN;

    end

end

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)+180;
elseif (prof.vel(j,1)>0 & prof.vel(j,2)<0)
theta(j)=360-theta(j);
elseif (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)=360-theta(j)-180;
end
end
prof.theta=theta;
prof.u=prof.vel(:,1);
prof.v=prof.vel(:,2);
%prof.w=prof.vel(:,6);
prof.w=ones(size(prof.vel(:,2)))*NaN;
prof.Vh=sqrt(prof.u.^2+prof.v.^2);
     
        