всё верно. Это я упустил указать (подробнее):
r1=sqrt((xi-x[S])^2+y[S]^2);
r2=sqrt((xi-x[R])^2+y[R]^2);
вот так соответственно.
Необходимо чтобы эта вся штука приняла вид
int( exp(2*i*k*z) * f(z) ) , z=A..B );
f дальше я бы с ним разобрался