function fig228dd

%  generate data for fig 228 

ep=0.05; b=0;  a=1;

% set parameters
nr = 200;

alpha=1;  beta=1;
chi=0.5*sqrt(alpha^2 + beta^2 - 4*ep)/ep;

nn=40;
pi2=2*pi;
b0 = (b-a)*ifc(0,pi2,0,ep,alpha,beta)/pi2;
bess0 = besseli(0,chi);
for in=1:nn
    ap(in)=(b-a)*ifs(0,pi2,in,ep,alpha,beta)/pi;
    bp(in)=(b-a)*ifc(0,pi2,in,ep,alpha,beta)/pi;
    bess(in)=besseli(in,chi);
end;

r=linspace(-1,1,nr);
phi=pi/4;

fid = fopen('asydata.txt', 'w');

for ir=1:nr
	x=r(ir)*cos(pi/4);
	y=r(ir)*sin(pi/4);
	rp=chi*r(ir);
	s=b0*besseli(0,rp)/bess0;
	for in=1:nn
		t1=besseli(in,rp)*(ap(in)*sin(in*phi)+bp(in)*cos(in*phi))/bess(in);
                s=s+t1;
	end;
	u(ir)=a+s*exp(-0.5*(alpha*x+beta*y)/ep);

fprintf(fid, '%12.6f %12.6f\n', r(ir), u(ir));

end;

% plot(r,u,'-','Linewidth',1)

fclose(fid);



function f=fs(phi,n,ep,alpha,beta)
f = exp(0.5*(alpha*cos(phi)+beta*sin(phi))/ep)*sin(n*phi);

function f=fc(phi,n,ep,alpha,beta)
f = exp(0.5*(alpha*cos(phi)+beta*sin(phi))/ep)*cos(n*phi);

function s=ifs(a,b,n,ep,alpha,beta);
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fs(xd(1),n,ep,alpha,beta);
for j=2:2:N
	f2=fs(xd(j+1),n,ep,alpha,beta);
	s=s+4*fs(xd(j),n,ep,alpha,beta)+2*f2;
end;
s=h*(s-f2)/3;

function s=ifc(a,b,n,ep,alpha,beta);
%  simpson's method
N=max(100,100*n);
xd=linspace(a,b,N+1);
h=xd(2)-xd(1);
s=fc(xd(1),n,ep,alpha,beta);
for j=2:2:N
	f2=fc(xd(j+1),n,ep,alpha,beta);
	s=s+4*fc(xd(j),n,ep,alpha,beta)+2*f2;
end;
s=h*(s-f2)/3;