clear;
vr = [6;8;0]; vr0 = 10*[1;1;sqrt(2)]; k=3,
//vr = [6;8;0]; vr0 = 4.5*[1;1;sqrt(2)]; k=3,
r = sqrt(vr'*vr), r0 = sqrt(vr0'*vr0), rr0=sqrt((vr0+vr)'*(vr0+vr));
x2r=0.5:(2*k*r+0.5);
//x2r=0.5:(1.5*k*r+0.5);
//
j2r = sqrt(%pi/(2*k*r))*besselj(x2r, k*r);
j2r0 = sqrt(%pi/(2*k*r0))*besselj(x2r, k*r0);
//
y2r0 = sqrt(%pi/(2*k*r0))*bessely(x2r, k*r0);
//
vectx = x2r; hankel = j2r0 + %i * y2r0; bessel = j2r; 
//
// Polynomes de Legendre de degres 0:m en x.
//
m = length(x2r)-1;
cosx = (vr'*vr0)/(r*r0); 
lege(1) = 1; lege(2) = cosx;
for i=1:m-1
  lege(i+2) = ((2*i+1)*cosx*lege(i+1) - i * lege(i)) /(i+1);
end
//
G(1) = hankel(1)*bessel(1)*lege(1);
MoinsUn = 1;
for i=2:m+1
  MoinsUn = - MoinsUn;
  G(i) = G(i-1) + MoinsUn * (2*i-1) * hankel(i)*bessel(i)*lege(i); 
end
G = ( (%i * k)/(4*%pi) ) * G ;
//
format("v",16);
Gexact = exp(%i * k * rr0) / (4 * %pi * rr0) ,
Gappr = G($),
han0 = ( (%i * k)/(4*%pi) ) * (  sqrt(%pi/(2*k*rr0)) * besselj(0.5,k*rr0)+%i * sqrt(%pi/(2*k*rr0)) * bessely(0.5,k*rr0)),
ErreurRel = abs(G-Gexact); 
xbasc; plot2d(vectx', ErreurRel, logflag="nl") , xtitle('GreenApprox');
//xbasc; plot2d(vectx', log(ErreurRel), logflag="nn") , xtitle('GreenApprox');
//