function [Vrm,scale,ee,UU,VV,WW]=vpair_2(isc,sca,SS1,SS2,sc1,xd1,zd1,sc2,xd2,zd2,xd,zd,az,el,Vreff,thetaprime,downwash)

if xd1>xd2
    aa=[xd1,zd1];
    xd1=xd2; zd1=zd2;
    xd2=aa(1);zd2=aa(2);
end
xdp=xd-xd1*ones(size(xd));
zdp=zd-zd1*ones(size(zd));
%Left vortex (clockwise rotation)
Xp=complex(xdp,zdp);
        rhop=abs(Xp);
        ang=angle(Xp);
        Vtheta=(1-exp(-1.26*((rhop/SS1).^2)))./(2*pi*rhop);
        up=Vtheta.*sin(ang);
        wp=-Vtheta.*cos(ang);
        
UUL=cosd(thetaprime)*up;
VVL=-sind(thetaprime)*up;
WWL=wp;

xdp=xd-xd2*ones(size(xd));
zdp=zd-zd2*ones(size(zd));

Xp=complex(xdp,zdp);
        rhop=abs(Xp);
        ang=angle(Xp);
        Vtheta=(1-exp(-1.26*((rhop/SS2).^2)))./(2*pi*rhop);
        up=-Vtheta.*sin(ang);
        wp=Vtheta.*cos(ang);
        
UUR=cosd(thetaprime)*up;
VVR=-sind(thetaprime)*up;
WWR=wp;

UU=sc1*UUL+sc2*UUR;
VV=sc1*VVL+sc2*VVR;
WW=sc1*WWL+sc2*WWR;
Vrm=sind(az).*cosd(el).*UU+cosd(az).*cosd(el).*VV+sind(el).*(WW+downwash*ones(size(WW)));
in=find(~isnan(Vreff));
if isc==0
scale=mean(abs(Vreff(in))./abs(Vrm(in)));
%scale=std(Vreff(in))/std(Vrm(in));
ee=pcorr(Vreff(in),Vrm(in));
else
    scale=sca;
    ee=sqrt(sum((Vreff(in)-Vrm(in)*scale).^2)/length(Vrm(in)));
    ee=1/ee;

end

%scale=std(Vreff(in))/std(Vrm(in));
return
end