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'); //