function [F,aF] = calcRF3(k0,ep1,mu1,sp,dvec,op)
% computes [del][del]{dvec gR}

k = k0*sqrt(ep1*mu1);

R = [ (op(1)-sp(1)) (op(2)-sp(2)) (op(3)+sp(3)) ];
r = sqrt( R(1)*R(1) + R(2)*R(2) + R(3)*R(3) );

tm1 = 1/r/r;
tm2 = (i*k-1/r);
tm3 = exp(i*k*r)/r;

F = (R*dvec')*tm3/r*(tm1 + tm2*(tm2 - 1/r))*R/r + tm2*tm3/r*dvec;
aF = (R*dvec')*tm3/r*tm2;