clear
L=100;
xmin=-L;
xmax=L;
n=90;
m=6;
ne=2*(n-2-1)*(m-1)+2*(m-1);
a=zeros(m*(n-2)+2);
b=zeros(m*(n-2)+2);
xy=zeros(m*(n-2)+2,2);
x=zeros(1,m*(n-2)+2);
y=zeros(1,m*(n-2)+2);
x(1)=xmin;
y(1)=uu(xmin);
for k=1:(n-1-1)
    for i=1:m
        x(i+1+(k-1)*m)=xmin+k*(xmax-xmin)/(n-1);
        t=x(i+(k-1)*m+1);
        y(i+1+(k-1)*m)=uu(t)+(i-1)*(vv(t)-uu(t))/(m-1);
    end
end
x(m*(n-2)+2)=xmax;
y(m*(n-2)+2)=uu(xmax);
clear k i
ne=2*(n-2-1)*(m-1)+2*(m-1);
a(1,1)=m;
a(m*(n-2)+2,m*(n-2)+2)=m;
a(1+m,1+m)=3;
a(m*(n-3)+2,m*(n-3)+2)=3;
for i=1:n-3
    a(2+m*(i-1),2+m*(i-1))=4;
    a(2*m+1+(i-1)*m,2*m+1+(i-1)*m)=4;
end
for j=1+1+1:1+m-1
    a(j,j)=5;
    a(m*(n-2)+3-j,m*(n-2)+3-j)=5;
end
clear i j
for i=1:n-4
    for j=1:m-2
        a(m+2+j+(i-1)*m,m+2+j+(i-1)*m)=6;
    end
end
a=sparse(a);
for i=1:m
    b(1,i+1)=-1;
    b(m*(n-2)+2-i,m*(n-2)+2)=-1;
end
clear i
for i=2:(n-2)*m+2-1-m
    b(i,i+m)=-1;
end
clear i
for i=2:n-1
    for j=1:m-1
        b(j+1+(i-2)*m,j+1+1+(i-2)*m)=-1;
    end
end
clear i j
for i=2:n-2
    for j=1:m-1
        b(j+1+(i-2)*m,j+1+1+m+(i-2)*m)=-1;
    end
end
clear i j
b=sparse(b);
b=a+b'+b;
clear a
a=sparse(b);
clear b
xy=[x',y'];
gplot(a,xy);
axis equal

e1=zeros(1,3);
u1=[0,0,0];
for k=1:m-1
    e1=[1,k+1,k+1+1];
    u1=[u1;e1];
end
clear k
e2=zeros(2,3);
for k=1:(n-3)
    for i=1:(m-1)
        e2=[i+1+(k-1)*m,i+1+k*m,k*m+1+1+i;i+1+(k-1)*m,k*m+1+1+i,i+1+(k-1)*m+1];
        u1=[u1;e2];
    end
end
clear k i 
for k=1:m-1
    e1=[(n-2)*m+2-m+k-1,(n-2)*m+2,(n-2)*m+2-m+1+k-1];
    u1=[u1;e1];
end
u1(1,:)=[];
clear k e1
N=u1'
min=0.01;
j=sqrt(-1);
ur=10;
er=1;
k0=8;
z0=15;
jz=7;
alphax=1/ur*ones(1,ne);
alphay=1/ur*ones(1,ne);
beta=-(k0).^2*er*ones(1,ne);
ff=-j*k0*z0*jz*ones(1,ne);

nd=zeros(1,n);
p=ones(1,n);
nd(1)=1;
for i=2:n
   nd(i)=(i-2)*m+2;
end
clear i
for i=1:n
    p(i)=exp(j*10*x(nd(i)));
end
clear i
ns=zeros(2,n-1);
for i=1:n-1
    ns(1,i)=1+m*(i-1);
    ns(2,i)=1+m*i;
end
ns(2,n-1)=ns(2,n-2)+1;
clear i
gamma=zeros(1,n-1);
clear i
for i=1:n-1
    gamma(i)=1/ur*(j*k0+k(ns(1,i))/2)
end
clear i 
q=zeros(1,n-1);
for i=1:n-1
    
    q(i)=gamma(i)*exp(j*k0*x(ns(1,i)))+1/ur*exp(j*k0*x(ns(1,i)))*j*k0*x(ns(1,i))/sqrt(1+4*(1-(x(ns(1,i))).^2/10000)/(21.1089.^2));
end
clear i

K=zeros(m*(n-2)+2,m*(n-2)+2);
B=zeros(1,m*(n-2)+2);
for e=1:ne
    i=N(1,e);
    j=N(2,e);
    mm=N(3,e);
    be(1)=y(j)-y(mm);
    be(2)=y(mm)-y(i);
    be(3)=y(i)-y(j);
    ce(1)=x(mm)-x(j);
    ce(2)=x(i)-x(mm);
    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 mm ii jj
for s=1:n-1
    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/6;
    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 s
for i=1:n
    for j=1:m*(n-2)+2
        if nd(i)==j
       K(j,j)=1;
      
        else 
        B(j)=B(j)-K(j,nd(i))*p(i);    
        K(j,nd(i))=0;
        K(nd(i),j)=0;
        end
    end
    
end
clear i j
for i=1:n
    B(nd(i))=p(i); 
end   
clear i
ans=inv(K)*B';





k0=2*pi; 
p=size(nd);
sizend=p(1,2);
dfaiin=zeros(1,2*sizend-2);
NN=2*sizend-2;
for i=1:sizend
   yx=x(nd(i));
    dfaiin(i)=j*k0*yx.*exp(j*k0*yx).*((51401035164868039/57646075230342348800-501956252569373/900719925474099200*i)*pi*exp(-14/25*i*yx*pi)-(71303374314635473/14411518807585587200+225974581611468529/57646075230342348800*i)*pi*exp(-47/50*i*yx*pi)+(243612760602038057/57646075230342348800-257323462866191897/115292150460684697600*i)*pi*exp(-29/50*i*yx*pi)-(73810570878776379/23058430092136939520+12018273693734983/11529215046068469760*i)*pi*exp(-41/50*i*yx*pi)+(143825810653222071/3689348814741910323200-29514238382438997/3689348814741910323200*i)*pi*exp(-33/50*i*yx*pi)+(-132568101814513339/28823037615171174400+92048547277357051/14411518807585587200*i)*pi*exp(-23/50*i*yx*pi)+(2923444717529347/115292150460684697600-45841459011290567/922337203685477580800*i)*pi*exp(-19/50*i*yx*pi)+(-3725285452569567/1801439850948198400+3121656575622129/288230376151711744*i)*pi*exp(-9/25*i*yx*pi)+(2749134074891505/1152921504606846976-750517040282643/720575940379279360*i)*pi*exp(-3/5*i*yx*pi)-(184534285534777799/922337203685477580800+68364810912748083/737869762948382064640*i)*pi*exp(-43/50*i*yx*pi)+(44225487806572671/14411518807585587200+89998036734315377/230584300921369395200*i)*pi*exp(-19/25*i*yx*pi)-(620822392532635/1152921504606846976+1607658398685779/11529215046068469760*i)*pi*exp(-4/5*i*yx*pi)+(16215417343644687/11529215046068469760+80326447161567309/92233720368547758080*i)*pi*exp(-9/10*i*yx*pi)-(58912950671988917/7205759403792793600+8475007924910387/1801439850948198400*i)*pi*exp(-7/50*i*yx*pi)+(-86019887384572427/28823037615171174400+1846353542106957/720575940379279360*i)*pi*exp(-13/25*i*yx*pi)-(897194693816681/1441151880758558720+1877624937212079/3602879701896396800*i)*pi*exp(-4/25*i*yx*pi)+(7003458733176275/9223372036854775808-2162520956497183/1152921504606846976*i)*pi*exp(-i*yx*pi)-(7380610983079133/2882303761517117440+583751867597161/180143985094819840*i)*pi*exp(-1/5*i*yx*pi)-(4505729151403781/900719925474099200+1880408665234849/7205759403792793600*i)*pi*exp(-2/25*i*yx*pi)+(20790613126386509/23058430092136939520-48199248319412701/737869762948382064640*i)*pi*exp(-7/10*i*yx*pi)+(5793042521749937/115292150460684697600-17605040462203393/46116860184273879040*i)*pi*exp(-17/50*i*yx*pi)+(-10004115476036157/922337203685477580800+40439415968764057/922337203685477580800*i)*pi*exp(-7/25*i*yx*pi)+(3761818683648417/900719925474099200+243003718595605463/922337203685477580800*i)*pi*exp(-37/50*i*yx*pi)+(2804299174760577/9223372036854775808+2201243511199179/720575940379279360*i)*pi*exp(-3/10*i*yx*pi)+(2623808684734677/720575940379279360+3654744229610619/450359962737049600*i)*pi*exp(-6/25*i*yx*pi)-(3098559313368035/2305843009213693952+834699494959099/1152921504606846976*i)*pi*exp(-22/25*i*yx*pi)+(-23607621562529007/28823037615171174400+140214190194774783/230584300921369395200*i)*pi*exp(-27/50*i*yx*pi)+(187494191214660291/57646075230342348800+144483985123546419/230584300921369395200*i)*pi*exp(-39/50*i*yx*pi)+(35913980564367939/57646075230342348800-11729657348275053/11529215046068469760*i)*pi*exp(-11/25*i*yx*pi)+(30115192046563793/14411518807585587200+47356976902472393/7205759403792793600*i)*pi*exp(-13/50*i*yx*pi)+(-2255330697246181/1152921504606846976+8958049265962987/4611686018427387904*i)*pi*exp(-1/2*i*yx*pi)-(2007659643590319/720575940379279360+2685993331110257/2882303761517117440*i)*pi*exp(-1/10*i*yx*pi)+(10672885937047533/7205759403792793600+15057125934423981/7205759403792793600*i)*pi*exp(-3/50*i*yx*pi)+(-58828277004357/576460752303423488+535435293724375/2305843009213693952*i)*pi*exp(-2/5*i*yx*pi)+(22180207525354269/7205759403792793600+19746965816879589/7205759403792793600*i)*pi*exp(-24/25*i*yx*pi)+(15046906762443397/461168601842738790400+162469221208110383/7378697629483820646400*i)*pi*exp(-23/25*i*yx*pi)-(18258664672733373/3602879701896396800+9725862371508561/7205759403792793600*i)*pi*exp(-3/25*i*yx*pi)+(-3333998204762847/900719925474099200+1849827284259779/1801439850948198400*i)*pi*exp(-16/25*i*yx*pi)+(7967157386289333/900719925474099200-700430319186897/1801439850948198400*i)*pi*exp(-1/25*i*yx*pi)-(307457712904524399/461168601842738790400+306964522280055087/461168601842738790400*i)*pi*exp(-49/50*i*yx*pi)+(15852342663856917/3602879701896396800+33271095034182393/7205759403792793600*i)*pi*exp(-9/50*i*yx*pi)+(19804304122715547/23058430092136939520-26551347958066219/230584300921369395200*i)*pi*exp(-17/25*i*yx*pi)+(42752165923155741/23058430092136939520-15289421497320319/23058430092136939520*i)*pi*exp(-31/50*i*yx*pi)+(683964540640053/46116860184273879040+350421986275223/144115188075855872*i)*pi*exp(-8/25*i*yx*pi)+(-16197818205569457/57646075230342348800+18587092495829649/57646075230342348800*i)*pi*exp(-12/25*i*yx*pi)-(18401637477978337/5764607523034234880+76912599906339219/14411518807585587200*i)*pi*exp(-11/50*i*yx*pi)+(127187241937371693/28823037615171174400+10031589417060729/5764607523034234880*i)*pi*exp(-21/25*i*yx*pi)+(2468484588722467/112589990684262400-945363304066419/225179981368524800*i)*pi*exp(-1/50*i*yx*pi)+(50065728563628657/14411518807585587200-5224192581618927/922337203685477580800*i)*pi*exp(-18/25*i*yx*pi)+(-21800916005368899/23058430092136939520+55531851334187823/28823037615171174400*i)*pi*exp(-21/50*i*yx*pi));
end 
for i=sizend+1:2*sizend-2 
    dfaiin(i)=j*k0*x(ns(1,2*sizend-i)).*exp(j*k0*x(ns(1,2*sizend-i))).*(-522588387703071/83886080/(8709806461717850000-870980646171785*x(ns(1,2*sizend-i)).^2)^(1/2)*x(ns(1,2*sizend-i)));
end 

clear i
for i=1:NN-1
    dfaiin(i)=(dfaiin(i)+ dfaiin(i+1))/2;
end
 dfaiin(NN)=(dfaiin(NN)+ dfaiin(1))/2;
 clear i
dfai=zeros(1,NN);
gammamiddle=zeros(1,NN);
for i=1:NN-1
    gammamiddle(i)= (gamma(i)+ gamma(i+1))/2;
end
 gamma(NN)= (gamma(NN)+ gamma(1))/2;
clear i

for i=1:NN-1
    dfai(i)=dfaiin(i)-gamma(i).*(f(nl(i))+f(nl(i+1)))/2
end
 dfai(NN)=dfaiin(NN)-gamma(NN).*(f(nl(NN))+f(nl(1)))/2
clear i
dl=zeros(1,NN);
for i=1:NN-1
    dl(i)=sqrt((x(nl(i))-x(nl(i+1))).^2+(y(nl(i))-y(nl(i+1))).^2);
end
dl(NN)=sqrt((x(nl(1))-x(nl(NN))).^2+(y(nl(1))-y(nl(NN))).^2);
clear i

yinta=120*pi;
syms th
rcs=0; 

for i=1:NN-1
   ppp=dfai(i).*dl(i).*exp(j*k0*(((x(nl(i))+x(nl(i+1)))/2*cos(th)+(y(nl(i))+y(nl(i+1)))/2*sin(th))));
   rcs=rcs+ppp;
end
 ppp=dfai(NN).*dl(NN).*exp(j*k0*(((x(nl(NN))+x(nl(1)))/2*cos(th)+(y(nl(NN))+y(nl(1)))/2*sin(th))));
 rcs=rcs+ppp;
 rcs=k0*yinta*yinta/4*(abs(rcs)).^2