function prof=vvp_regression_RHI(rhod,zd,rws,snr,ranges,az,el,zh)

[na nr]=size(rws);
ind=1;
oops=(na-1)*nr;
rhod=reshape(rhod(1:na-1,:),oops,1);
zd=reshape(zd(1:na-1,:),oops,1);
vr=reshape(rws(1:na-1,:),oops,1);
snr=reshape(snr(1:na-1,:),oops,1);
el=reshape(el(1:na-1,:),oops,1);
az=reshape(az(1:na-1,:),oops,1);
drange=double(ranges(2)-ranges(1));
prof.heights=zh;
dheight=5*drange;
xpos=0.;
ypos=0;
rho=10000;
nheights=length(prof.heights);
for j=1:length(prof.heights)
    zpp=prof.heights(j);
    ind=find((abs(zd-prof.heights(j))<=dheight/2) & (~isnan(vr)) );
    np=2;
    na=length(ind);
    if (na>np)
        clear Cook H resid X vel ss DFFITS
        prof.ninda(j)=length(ind);
        rhodt=rhod(ind);
        zdt=zd(ind);
        elt=el(ind);
        azt=az(ind);
        vrt=vr(ind);
        X=pred_RHI(rhodt,zdt,elt*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);
            rhodtt=rhodt(inda);
            zdtt=zdt(inda);
            eltt=elt(inda);
            aztt=azt(inda);
            vrtt=vrt(inda);

            X=pred_RHI(rhodtt,zdtt,eltt*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);
                rhodttt=rhodtt(indb);
                
                zdttt=zdtt(indb);
                elttt=eltt(indb);
                azttt=aztt(indb);
                vrttt=vrtt(indb);

                X=pred_RHI(rhodttt,zdttt,elttt*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,np)*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;
prof.rho=prof.vel(:,1);
prof.w=prof.vel(:,2);
     
        