% Plots of some RBF
close all
clear all


% s=input('give me the dimension = '); 
% n=input('give me the order = ');
e=3;
neval = 100; 
gr = linspace(-1,1,neval);
[x,y]=meshgrid(gr,gr);
xy=sqrt(x.^2+y.^2);

% Gaussian-Laguerre
gl11 = inline('(3/2-(e*r).^2).*exp(-(e*r).^2)'); 
gl21 = inline('(2-(e*r).^2).*exp(-(e*r).^2)');
gl12 = inline('(15/8-5/2*(e*r).^2+1/2*(e*r).^4).*exp(-(e*r).^2)'); 
gl22 = inline('(3-3*(e*r).^2+1/2*(e*r).^4).*exp(-(e*r).^2)');

rbf11=feval(gl11,e,gr);
rbf12=feval(gl21,e,gr);

rbf21=feval(gl12,e,xy);
rbf22=feval(gl22,e,xy);

figure(1);
plot(gr, rbf11,'r-');
grid

figure(2);
mesh(x,y,rbf21);

figure(3);
plot(gr,rbf12,'r-');
grid

figure(4);
mesh(x,y,rbf22);

close all

% Poisson

J0=besselj(0,10*xy);
J1=sqrt(2/pi)*sin(10*xy)./(10*xy+eps);
J2=besselj(1,10*xy)./(10*xy+eps);

figure(1);
mesh(x,y,J0); 

figure(2) 
mesh(x,y,J1);
figure(3) 
mesh(x,y,J2);

close all

% Matern
mB1=exp(-5*xy);
mB2=(1+5*xy).*mB1;
mB3=(3+15*xy+25*xy.^2).*mB1;

figure(1);
mesh(x,y,mB1); 

figure(2) 
mesh(x,y,mB2);
figure(3) 
mesh(x,y,mB3);

close all

% Multiquadriche inverse
IM1=(1+25*xy.^2).^(-1/2);
IM2=1./(1+25*xy.^2);

figure(1);
mesh(x,y,IM1); 

figure(2) 
mesh(x,y,IM2);

close all
% Truncated powers

TP1=max((1-xy),0).^2;
TP2=max((1-xy),0).^4;
figure(1);
mesh(x,y,TP1); 

figure(2) 
mesh(x,y,TP2);
close all

% Whittaker potentials

W02=1-xy+xy.*exp(-1./xy);

W12=1-2*xy+(1+2*xy).*exp(-1./xy);

figure(1);
mesh(x,y,W02); 

figure(2) 
mesh(x,y,W12);

