Function X,Y =Nacircleg(Bs,S)

function [x,y]=nacircleg(bs,s)
nbs=12;
h=3;
hh=0.5;
if nargin==0,x=nbs;return;end
d=[-h/2 h/2-h/2 h/2+hh h/2-h/2h/2-h/2-hh -h/2 -h/2 h/2h/2
h/2h/2+hh h/2h/2-h/2 -h/2-hh -h/2 -h/2h/2h/2-h/2 -h/2
111111 11 0000
000000 00 1111];
if nargin==1,x=d(:,bs);return;end
x=zeros(size(s));y=zeros(size(s));
[m,n]=size(bs);
if m==1&n==1,bs=bs*ones(size(s));end
ii=find(bs==1);
x(ii)=s(ii);
y(ii)=-h/2-hh;
ii=find(bs==2);
x(ii)=s(ii);
t=s(ii);
y(ii)=-sqrt(hh.^2-(t-h/2).^2)-h/2;
clear t
ii=find(bs==3);
x(ii)=h/2+hh;
y(ii)=s(ii);
ii=find(bs==4);
x(ii)=s(ii);
t=s(ii);
y(ii)=sqrt(hh.^2-(t-h/2).^2)+h/2;
clear t
ii=find(bs==5);
x(ii)=s(ii);
y(ii)=h/2+hh;
ii=find(bs==6);
x(ii)=s(ii);
t=s(ii);
y(ii)=sqrt(hh.^2-(t+h/2).^2)+h/2;
clear t
ii=find(bs==7);
x(ii)=-h/2-hh;
y(ii)=s(ii);
ii=find(bs==8);
x(ii)=s(ii);
t=s(ii);
y(ii)=-sqrt(hh.^2-(t+h/2).^2)-h/2;
clear t
ii=find(bs==9);
x(ii)=s(ii);
y(ii)=-h/2;
ii=find(bs==10);
x(ii)=h/2;
y(ii)=s(ii);
ii=find(bs==11);
x(ii)=s(ii);
y(ii)=h/2;
ii=find(bs==12);
x(ii)=-h/2;
y(ii)=s(ii);
axis equal;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[p,e,t]=initmesh('nacircleg','hmax',inf);
subplot(2,2,1), pdemesh(p,e,t)
[p,e,t]=refinemesh('nacircleg',p,e,t);
subplot(2,2,2), pdemesh(p,e,t)
[p,e,t]=refinemesh('nacircleg',p,e,t);
subplot(2,2,3), pdemesh(p,e,t)
[p,e,t]=refinemesh('nacircleg',p,e,t);
subplot(2,2,4), pdemesh(p,e,t)
axis equal
p=p';
min=0.01;
j=sqrt(-1);
x=p(:,1);
y=p(:,2);
enumber=1024
clear p
N=size(x);
ur=10;
er=1;
k0=8;
z0=15;
jz=7;
alphax=1/ur*ones(1,enumber);
alphay=1/ur*ones(1,enumber);
beta=-(k0).^2*er*ones(1,enumber);
ff=-j*k0*z0*jz*ones(1,enumber);
n=zeros(3,1024);
n(1,:)=t(1,:);
n(2,:)=t(2,:);
n(3,:)=t(3,:);
nodes=size(x);
for i=1:nodes
if (y(i)==-1.5&(x(i)<=1.5&x(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
nd=nonzeros(pp);
for i=1:size(nd)-1
for k=1:size(nd)-1
if x(nd(k))>x(nd(k+1))
tt=nd(k);
nd(k)=nd(k+1);
nd(k+1)=tt;
end
end
end
clear i k tt pp
for i=1:nodes
if (x(i)==1.5&(y(i)<=1.5&y(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
ndd=nonzeros(pp);
for i=1:size(ndd)-1
for k=1:size(ndd)-1
if y(ndd(k))>y(ndd(k+1))
tt=ndd(k);
ndd(k)=ndd(k+1);
ndd(k+1)=tt;
end
end
end
nd=[nd;ndd];
clear i k tt pp ndd
for i=1:nodes
if (y(i)==1.5&(x(i)<=1.5-3/8&x(i)>=-1.5))
pp(i)=i;
end
end
clear i
ndd=nonzeros(pp);
for i=1:size(ndd)-1
for k=1:size(ndd)-1
if x(ndd(k))<x(ndd(k+1))
tt=ndd(k);
ndd(k)=ndd(k+1);
ndd(k+1)=tt;
end
end
end
nd=[nd;ndd];
clear i k tt pp ndd
for i=1:nodes
if (x(i)==-1.5&(y(i)<=1.5-3/8&y(i)>=-1.5))
pp(i)=i;
end
end
clear i
ndd=nonzeros(pp);
for i=1:size(ndd)-1
for k=1:size(ndd)-1
if y(ndd(k))<y(ndd(k+1))
tt=ndd(k);
ndd(k)=ndd(k+1);
ndd(k+1)=tt;
end
end
end
nd=[nd;ndd];
clear i k tt pp ndd
for i=1:nodes
if (y(i)==-2&(x(i)<=1.5-3/8&x(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
nl=nonzeros(pp);
for i=1:size(nl)-1
for k=1:size(nl)-1
if x(nl(k))>x(nl(k+1))
tt=nl(k);
nl(k)=nl(k+1);
nl(k+1)=tt;
end
end
end
clear i k tt pp
for i=1:nodes
if abs((x(i)-1.5).^2+(y(i)+1.5).^2-0.25)<min
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if x(nll(k))>x(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if (x(i)==2&(y(i)<=1.5-3/8&y(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if y(nll(k))>y(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if abs((x(i)-1.5).^2+(y(i)-1.5).^2-0.25)<min
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if y(nll(k))>y(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if (y(i)==2&(x(i)<=1.5-3/8&x(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if x(nll(k))<x(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if abs((x(i)+1.5).^2+(y(i)-1.5).^2-0.25)<min
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if x(nll(k))<x(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if (x(i)==-2&(y(i)<=1.5-3/8&y(i)>=3/8-1.5))
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if y(nll(k))<y(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll
for i=1:nodes
if abs((x(i)+1.5).^2+(y(i)+1.5).^2-0.25)<min
pp(i)=i;
end
end
clear i
nll=nonzeros(pp);
for i=1:size(nll)-1
for k=1:size(nll)-1
if y(nll(k))<y(nll(k+1))
tt=nll(k);
nll(k)=nll(k+1);
nll(k+1)=tt;
end
end
end
nl=[nl;nll];
clear i k tt pp nll min
ns=zeros(2,size(nl));
for i=1:size(nl)-1
ns(2,i)=nl(i+1,1);
ns(1,i)=nl(i,1);
end
ns(1,size(nl))=nl(size(nl));
ns(2,size(nl))=nl(1);
ns(2,1)=nl(2);
clear i
pe=zeros(1,size(nd));
for i=1:size(nd)
pe(i)=exp(j*10*x(nd(i)));
end
clear i
gamma=1/ur*[j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,...
j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,...
j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,...
j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1,j*k0+1];
q=zeros(1,size(ns));
for i=1:7
q(i)=gamma(i)*exp(j*k0*x(nl(i)));
end
for i=8:24
q(i)=gamma(i)*exp(j*k0*x(nl(i)))+2*(x(i)-1.5)*j*k0*exp(j*k0*x(nl(i)));
end
for i=25:31
q(i)=gamma(i)*exp(j*k0*x(nl(i)))+j*k0*exp(j*k0*x(nl(i)));
end
for i=32:48
q(i)=gamma(i)*exp(j*k0*x(nl(i)))+2*(x(i)-1.5)*j*k0*exp(j*k0*x(nl(i)));
end
for i=48+1:48+7
q(i)=gamma(i)*exp(j*k0*x(nl(i)));
end
for i=48+8:48+24
q(i)=gamma(i)*exp(j*k0*x(nl(i)))+2*(x(i)+1.5)*j*k0*exp(j*k0*x(nl(i)));
end
for i=48+25:48+31
q(i)=gamma(i)*exp(j*k0*x(nl(i)))-j*k0*exp(j*k0*x(nl(i)));
end
for i=48+32:48+48
q(i)=gamma(i)*exp(j*k0*x(nl(i)))+2*(x(i)+1.5)*j*k0*exp(j*k0*x(nl(i)));
end
K=zeros(N,N);
B=zeros(1,N);
for e=1:enumber
i=n(1,e);
j=n(2,e);
m=n(3,e);
be(1)=y(j)-y(m);
be(2)=y(m)-y(i);
be(3)=y(i)-y(j);
ce(1)=x(m)-x(j);
ce(2)=x(i)-x(m);
ce(3)=x(j)-x(i);
deltae=0.5*(be(1)*ce(2)-be(2)*ce(1));
for ii=1:3
for jj=1:3
if(ii==jj)
del_ij=1;
else
del_ij=0;
end
ke(ii,jj)=(alphax(e)*be(ii)*be(jj)+...
alphay(e)*ce(ii)*ce(jj))/(4*deltae)+...
beta(e)*(1+del_ij)*deltae/12;
K(n(ii,e),n(jj,e))=K(n(ii,e),n(jj,e))+ke(ii,jj);
end
bbe(ii)=ff(e)*deltae/3;
B(n(ii,e))=B(n(ii,e))+bbe(ii);
end
end
clear e i j m ii jj
for s=1:size(nl)
i=ns(1,s);
j=ns(2,s);
ls=sqrt((x(i)-x(j)).^2+(y(i)-y(j)).^2);
ks(1,1)=gamma(s)*ls/3;
ks(1,2)=gamma(s)*ls/3;
ks(2,1)=ks(1,2);
ks(2,2)=ks(1,1);
for ii=1:2
for jj=1:2
K(ns(ii,s),ns(jj,s))=K(ns(ii,s),ns(jj,s))+ks(ii,jj);
end
bbs(ii)=q(s)*ls/2;
B(ns(ii,s))=B(ns(ii,s))+bbs(ii);
end
end
clear i j ii jj
for i=1:size(nd)
for j=1:size(x)
if nd(i)==j
K(j,j)=1;
else
B(j)=B(j)-K(j,nd(i))*pe(i);
K(j,nd(i))=0;
K(nd(i),j)=0;
end
end
end
clear i j
for i=1:size(nd)
B(nd(i))=pe(i);
end
clear i
inv(K)
B'
f=inv(K)*B'