function pile21

clear all;   
close all;

% global  y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 m11 c1 k1 P0 M0

r0=0.18;l=10; Ep=2.7e5;rm=20*r0;nth=100;v= 0.2;M0= 0.*r0;Ap=pi.*(r0.^2);
N=100;nn=-0.5;L=l;Ip=(pi./4).*r0^4;AA=600;P0=30;A=Ep.*Ip;
G01=load('axial1.txt');
G02=load('axial2.txt');
G03=load('axial3.txt');
G04=load('axial4.txt');

z=zeros(1,N);
z0=0; z11=3;z22=4; z33=8 ;z44=l; n1=N; n2=2*N; n3=3*N; n4=4*N;dz=L/(4*N);
eps=zeros(4,N);
phi=zeros(4,N);

r = logspace(log10(r0),log10(rm),N);
z1=linspace(0,z11,N);
z2=linspace(z11,z22,N);
z3=linspace(z22,z33,N);
z4=linspace(z33,z44,N);
rt=r;
%coefficients for eps 
eps1=0.000001*ones(1,N);
eps2=0.000001*ones(1,N);
eps3=0.000001*ones(1,N);
eps4=0.000001*ones(1,N);

t1=dfdx(eps1(1,:),z1,N);
t2=dfdx(eps2(1,:),z2,N);
t3=dfdx(eps3(1,:),z3,N);
t4=dfdx(eps4(1,:),z4,N);
%phi(1)=phi_r,phi(2)=phi_theta,phi(3)=diff(phi_r),phi(4)=diff(phi_theta),

y1new=ones(1,N); y2new=ones(1,N); y3new=ones(1,N); y4new=ones(1,N);....
    y5new=ones(1,N);y6new=ones(1,N);y7new=ones(1,N);y8new=ones(1,N);
y1old=0.01*ones(1,N);y2old=0.01*ones(1,N);y3old=0.01*ones(1,N);....
    y4old=0.01*ones(1,N);y5old=0.01*ones(1,N);y6old=0.01*ones(1,N);....
    y7old=0.01*ones(1,N);y8old=0.01*ones(1,N);
tol=0.001;   %tolerance
while(norm((y1old-y1new)./y1old))>tol||(norm((y2old-y2new)./y2old))>tol....
        ||(1/sqrt(N)*norm((y3old-y3new)./y3old))>tol||(1/sqrt(N)*...
        norm((y4old-y4new)./y4old))>tol....
        ||(1/sqrt(N)*norm((y5old-y5new)./y5old))>tol||(1/sqrt(N)*...
        norm((y6old-y6new)./y6old))>tol...
        ||(1/sqrt(N)*norm((y7old-y7new)./y7old))>tol||(1/sqrt(N)*...
        norm((y8old-y8new)./y8old))>tol;
y1=y1new;y2=y2new;y3=y3new;y4=y4new;y5=y5new;y6=y6new;y7=y7new;y8=y8new;
[1/sqrt(N)*norm((y1old-y1new)./y1old),1/sqrt(N)*...
    norm((y2old-y2new)./y2old),1/sqrt(N)*norm((y3old-y3new)./y3old),....
    1/sqrt(N)*norm((y4old-y4new)./y4old),1/sqrt(N)*...
    norm((y5old-y5new)./y5old),1/sqrt(N)*norm((y6old-y3new)./y6old),...
1/sqrt(N)*norm((y7old-y7new)./y7old),1/sqrt(N)*norm((y8old-y8new)./y8old)];
%--------------------------------------------------------------------------
solinit=bvpinit(logspace(log10(r0),log10(rm),N),[0.0001 1 -1 1]);
sol = bvp4c(@(r,phi)phiode(r,phi,rt,y1,y2,y3,y4,y5,y6,y7,y8)...
    ,@phibc,solinit);
phi = deval(sol,r);
%--------------------------------------------------------------------------
m11=Ep.*Ip;
A=m11;
dtheta=2*pi/nth;
theta=-dtheta/2;

m11=-Ep.*Ip;
C1=0;
k1=0;
z0=0;
%coefficients for eps 
for i=1:nth
    theta=theta+dtheta;
    mu1=Gs11(z0,z1,rt,theta,phi,eps1,t1,AA,nn,G01(i)) 
    size(mu1)
    lamda1=2.*mu1.*v./(1-2.*v);
   C1=trapz(r,(r.*mu1.*(phi(1,:).^2.*(cos(theta)).^2+phi(2,:).^2.*...
       (sin(theta)).^2))).*dtheta+C1;
   k1=trapz(r,((lamda1+2.*mu1).*(phi(3,:)).^2.*(cos(theta)).^2....
        +2.*lamda1.*1./r.*phi(3,:).*(phi(1,:)-phi(2,:)).*(cos(theta)).^2....
        +(2.*mu1+lamda1).*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(cos(theta)).^2....
        +mu1.*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(sin(theta)).^2....
        +2.*mu1.*1./r.*phi(4,:).*(phi(1,:)-phi(2,:)).*(sin(theta)).^2....
        +mu1.*(phi(4,:)).^2.*(sin(theta)).^2).*r).*dtheta+k1;
end

C2=0;
k2=0;
for i=1:nth
    theta=theta+dtheta;
    mu2=Gs12(z11,z2,rt,theta,phi,eps2,t2,AA,nn,G02(i));
    lamda2=2.*mu2.*v./(1-2.*v);
   C2=trapz(r,(r.*mu2.*(phi(1,:).^2.*(cos(theta)).^2+phi(2,:).^2.*...
       (sin(theta)).^2))).*dtheta+C2;
   k2=trapz(r,((lamda2+2.*mu2).*(phi(3,:)).^2.*(cos(theta)).^2....
        +2.*lamda2.*1./r.*phi(3,:).*(phi(1,:)-phi(2,:)).*(cos(theta)).^2....
        +(2.*mu2+lamda2).*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(cos(theta)).^2....
        +mu2.*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(sin(theta)).^2....
        +2.*mu2.*1./r.*phi(4,:).*(phi(1,:)-phi(2,:)).*(sin(theta)).^2....
        +mu2.*(phi(4,:)).^2.*(sin(theta)).^2).*r).*dtheta+k2;
end
C3=0;
k3=0;
for i=1:nth
    theta=theta+dtheta;
    mu3=Gs13(z22,z3,rt,theta,phi,eps3,t3,AA,nn,G03(i));
    lamda3=2.*mu3.*v./(1-2.*v);
   C3=trapz(r,(r.*mu3.*(phi(1,:).^2.*(cos(theta)).^2+phi(2,:).^2.*...
       (sin(theta)).^2))).*dtheta+C3;
   k3=trapz(r,((lamda3+2.*mu3).*(phi(3,:)).^2.*(cos(theta)).^2....
        +2.*lamda3.*1./r.*phi(3,:).*(phi(1,:)-phi(2,:)).*(cos(theta)).^2....
        +(2.*mu3+lamda3).*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(cos(theta)).^2....
        +mu3.*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(sin(theta)).^2....
        +2.*mu3.*1./r.*phi(4,:).*(phi(1,:)-phi(2,:)).*(sin(theta)).^2....
        +mu3.*(phi(4,:)).^2.*(sin(theta)).^2).*r).*dtheta+k3;
end
C4=0;
k4=0;
for i=1:nth
    theta=theta+dtheta;
    mu4=Gs14(z33,z4,rt,theta,phi,eps4,t4,AA,nn,G04(i));
    lamda4=2.*mu4.*v./(1-2.*v);
   C4=trapz(r,(r.*mu4.*(phi(1,:).^2.*(cos(theta)).^2+phi(2,:).^2.*...
       (sin(theta)).^2))).*dtheta+C4;
   k4=trapz(r,((lamda4+2.*mu4).*(phi(3,:)).^2.*(cos(theta)).^2....
        +2.*lamda4.*1./r.*phi(3,:).*(phi(1,:)-phi(2,:)).*(cos(theta)).^2....
        +(2.*mu4+lamda4).*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(cos(theta)).^2....
        +mu4.*1./r.^2.*(phi(1,:)-phi(2,:)).^2.*(sin(theta)).^2....
        +2.*mu4.*1./r.*phi(4,:).*(phi(1,:)-phi(2,:)).*(sin(theta)).^2....
        +mu4.*(phi(4,:)).^2.*(sin(theta)).^2).*r).*dtheta+k4;
end
%phi coefficients 

ms11=zeros(1,N);ms12=zeros(1,N);ms13=zeros(1,N);ns11=zeros(1,N);...
    ns12=zeros(1,N);

for i=1:N
    theta=-dtheta/2;
    for j=1:nth
        theta=theta+dtheta;
%         G01=load('axial1.txt');

      mu1=Gs21(rt(i),rt,theta,phi,eps1,t1,AA,nn,G01(i));
        lamda1=2.*mu1.*v./(1-2.*v);
      ms11(i)=trapz(z1,(lamda1+2*mu1).*eps1(1,:).^2.*(cos(theta)).^2.*...
            rt(i)).*dtheta+ms11(i);
      ms12(i)=trapz(z1,mu1.*eps1(1,:).^2.*(sin(theta)).^2.*rt(i)).*...
            dtheta+ms12(i);
      ms13(i)=trapz(z1,lamda1.*eps1(1,:).^2.*(cos(theta)).^2.*rt(i)).*...
            dtheta+ms13(i);
      ns11(i)=trapz(z1,mu1.*t1.^2.*(cos(theta)).^2.*rt(i)).*dtheta+ns11(i);
      ns12(i)=trapz(z1,mu1.*t1.^2.*(sin(theta)).^2.*rt(i)).*dtheta+ns12(i);
    end
end

ms21=zeros(1,N);ms22=zeros(1,N);ms23=zeros(1,N);ns21=zeros(1,N);...
    ns22=zeros(1,N);

for i=1:N
    theta=-dtheta/2;
    for j=1:nth
        theta=theta+dtheta;
%            G02=load('axial2.txt');

      mu2=Gs22(rt(i),rt,theta,phi,eps2,t2,AA,nn,G02(i));
        lamda2=2.*mu2.*v./(1-2.*v);
      ms21(i)=trapz(z2,(lamda2+2*mu2).*eps2(1,:).^2.*(cos(theta)).^2.*...
          rt(i)).*dtheta+ms21(i);
      ms22(i)=trapz(z2,mu2.*eps2(1,:).^2.*(sin(theta)).^2.*rt(i)).*...
          dtheta+ms22(i);
      ms23(i)=trapz(z2,lamda2.*eps2(1,:).^2.*(cos(theta)).^2.*rt(i)).*...
          dtheta+ms23(i);
      ns21(i)=trapz(z2,mu2.*t2.^2.*(cos(theta)).^2.*rt(i)).*dtheta+ns21(i);
      ns22(i)=trapz(z2,mu2.*t2.^2.*(sin(theta)).^2.*rt(i)).*dtheta+ns22(i);
    end
end

ms31=zeros(1,N);ms32=zeros(1,N);ms33=zeros(1,N);ns31=zeros(1,N);...
    ns32=zeros(1,N);

for i=1:N
    theta=-dtheta/2;
    for j=1:nth
        theta=theta+dtheta;
  
      mu3=Gs23(rt(i),rt,theta,phi,eps3,t3,AA,nn,G03(i));
        lamda3=2.*mu3.*v./(1-2.*v);
      ms31(i)=trapz(z3,(lamda3+2*mu3).*eps3(1,:).^2.*(cos(theta)).^2.*...
            rt(i)).*dtheta+ms31(i);
      ms32(i)=trapz(z3,mu3.*eps3(1,:).^2.*(sin(theta)).^2.*rt(i)).*...
            dtheta+ms32(i);
      ms33(i)=trapz(z3,lamda3.*eps3(1,:).^2.*(cos(theta)).^2.*rt(i)).*...
            dtheta+ms33(i);
      ns31(i)=trapz(z3,mu3.*t3.^2.*(cos(theta)).^2.*rt(i)).*dtheta+ns31(i);
      ns32(i)=trapz(z3,mu3.*t3.^2.*(sin(theta)).^2.*rt(i)).*dtheta+ns32(i);
    end
end
ms41=zeros(1,N);ms42=zeros(1,N);ms43=zeros(1,N);ns41=zeros(1,N);...
    ns42=zeros(1,N);

for i=1:N
    theta=-dtheta/2;
    for j=1:nth
        theta=theta+dtheta;
 
      mu4=Gs24(rt(i),rt,theta,phi,eps4,t4,AA,nn,G04(i));
        lamda4=2.*mu4.*v./(1-2.*v);
      ms41(i)=trapz(z4,(lamda4+2*mu4).*eps4(1,:).^2.*(cos(theta)).^2.*...
            rt(i)).*dtheta+ms41(i);
      ms42(i)=trapz(z4,mu4.*eps4(1,:).^2.*(sin(theta)).^2.*rt(i)).*...
            dtheta+ms42(i);
      ms43(i)=trapz(z4,lamda4.*eps4(1,:).^2.*(cos(theta)).^2.*rt(i)).*...
            dtheta+ms43(i);
      ns41(i)=trapz(z4,mu4.*t4.^2.*(cos(theta)).^2.*rt(i)).*dtheta+ns41(i);
      ns42(i)=trapz(z4,mu4.*t4.^2.*(sin(theta)).^2.*rt(i)).*dtheta+ns42(i);
    end
end


ms1=zeros(1,N);ms2=zeros(1,N);ms3=zeros(1,N);ns1=zeros(1,N);ns2=zeros(1,N);
ms1=ms11+ms21+ms31+ms41;
ms2=ms12+ms22+ms32+ms42;
ms3=ms13+ms23+ms33+ms43;
ns1=ns11+ns21+ns31+ns41;
ns2=ns12+ns22+ns32+ms42;

B=zeros(n2,n2);
B(1,1)=(2*A/(dz.^4))+(2*C1/(dz.^2))+k1;  %equ1
B(1,2)=((-4*A/(dz.^4))-(2*C1/(dz.^2)));
B(1,3)=(2*A/(dz.^4));

B(2,1)=(-2*A/(dz.^4))-(C1/(dz.^2)); %equ2
B(2,2)=(5*A/(dz.^4))+(2*C1/(dz.^2))+k1;
B(2,3)=(-4*A/(dz.^4))-(C1/(dz.^2));
B(2,4)=(A/(dz.^4));

B(n1-1,n1-3)=(A/(dz.^4)); %equ5
B(n1-1,n1-2)=(-4*A/(dz.^4))-(C1/(dz.^2));
B(n1-1,n1-1)=(6*A/(dz.^4))+(2*C1/(dz.^2))+k1;
B(n1-1,n1)=((-4*A/(dz.^4))-(C1/(dz.^2)));
B(n1-1,n1+1)=((-4*A/(dz.^4))-(C1/(dz.^2)));
B(n1-1,n1+2)=(A/(dz.^4));
B(n1,n1)=1; %equ6
B(n1,n1+1)=-1;
%--------------------------------------------------------------------------
B(n1+1,n1-1)=(A/(dz.^4)); %equ7
B(n1+1,n1)=((-4*A/(dz.^4))-(C2/(dz.^2)));
B(n1+1,n1+1)=((-4*A/(dz.^4))-(C2/(dz.^2)));
B(n1+1,n1+2)=((6*A/(dz.^4))+(2*C2/(dz.^2))+k2);
B(n1+1,n1+3)=((-4*A/(dz.^4))-(C2/(dz.^2)));
B(n1+1,n1+4)=(A/(dz.^4));

B(n2-1,n2-3)=A/(dz.^4); %equ10
B(n2-1,n2-2)=(-4.*A/(dz.^4))-(C2/(dz.^2));
B(n2-1,n2-1)=(6.*A/(dz.^4))+(2.*C2./(dz.^2))+k2;
B(n2-1,n2)=(-4.*A/(dz.^4))-(C2./(dz.^2));
B(n2-1,n2+1)=(-4.*A/(dz.^4))-(C2./(dz.^2));
B(n2-1,n2+2)=A/(dz.^4);

B(n2,n2)=1; %equ11
B(n2,n2+1)=-1;
%--------------------------------------------------------------------------
B(n2+1,n2-1)=A/(dz.^4); %equ12
B(n2+1,n2)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n2+1,n2+1)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n2+1,n2+2)=(6.*A/(dz.^4))+(2.*C3./(dz.^2))+k3;
B(n2+1,n2+3)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n2+1,n2+4)=A/(dz.^4);

B(n3-1,n3-3)=A/(dz.^4);%equ15
B(n3-1,n3-2)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n3-1,n3-1)=(6.*A/(dz.^4))+(2.*C3./(dz.^2))+k3;
B(n3-1,n3)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n3-1,n3+1)=(-4.*A/(dz.^4))-(C3./(dz.^2));
B(n3-1,n3+2)=A/(dz.^4);

B(n3,n3)=1;%equ16
B(n3,n3+1)=-1;
%--------------------------------------------------------------------------

B(n3+1,n3-1)=A/(dz.^4);%equ17
B(n3+1,n3)=(-4.*A/(dz.^4))-(C4./(dz.^2));
B(n3+1,n3+1)=(-4.*A/(dz.^4))-(C4./(dz.^2));
B(n3+1,n3+2)=(6.*A/(dz.^4))+(2.*C4./(dz.^2))+k4;
B(n3+1,n3+3)=(-4.*A/(dz.^4))-(C4./(dz.^2));
B(n3+1,n3+4)=A/(dz.^4);

B(n4-1,n4-3)=A/(dz.^4);%equ20
B(n4-1,n4-2)= (-4.*A/(dz.^4))-(C4./(dz.^2));
B(n4-1,n4-1)=(5.*A/(dz.^4))+(2.*C4./(dz.^2))+k4;
B(n4-1,n4)= (-2.*A/(dz.^4))-(C4./(dz.^2));

B(n4,n4)=1;%equ21
%--------------------------------------------------------------------------
%equ3 and %equ4

for i=3:n1-2
    for j=i-2:i-2
    B(i,j)=(A/(dz.^4));
    end
end
 for i=3:n1-2
    for j=i-1:i-1
    B(i,j)=((-4*A/(dz.^4))-(C1/(dz.^2)));
end
end

for i=3:n1-2
    for j=i:i
     B(i,j)=((6*A/(dz.^4))+(2*C1/(dz.^2))+k1);
end
end
for i=3:n1-2
    for j=i+1:i+1
        B(i,j)=((-4*A/(dz.^4))-(C1/(dz.^2)));
    end
end

for i=3:n1-2
    for j=i+2:i+2
     B(i,j)=(A/(dz.^4));
    end
end
%--------------------------------------------------------------------------
 for i=n1+2:n2-2
    for j=i-2:i-2
    B(i,j)=A./(dz.^4);
end
end

 for i=n1+2:n2-2
    for j=i-1:i-1
    B(i,j)=-4.*A./(dz.^4)-C2./(dz.^2);
end
end

for i=n1+2:n2-2
    for j=i:i
     B(i,j)=6.*A./(dz.^4)+2.*C2./(dz.^2)+k2;
end
end
for i=n1+2:n2-2
    for j=i+1:i+1
        B(i,j)=-4.*A./(dz.^4)-C2./(dz.^2);
end
end

for i=n1+2:n2-2
    for j=i+2:i+2
     B(i,j)=(A./(dz.^4));
    end
end
%--------------------------------------------------------------------------
for i=n2+2:n3-2
    for j=i-2:i-2
    B(i,j)=A/(dz.^4);
    end
end
 for i=n2+2:n3-2
    for j=i-1:i-1
    B(i,j)=-4.*A./(dz.^4)-C3./(dz.^2);
end
end

for i=n2+2:n3-2
    for j=i:i
     B(i,j)=(6.*A./(dz.^4))+2.*(C3./(dz.^2))+k3;
end
end
for i=n2+2:n3-2
    for j=i+1:i+1
        B(i,j)=(-4.*A./(dz.^4))-(C3./(dz.^2));
end
end

for i=n2+2:n3-2
    for j=i+2:i+2
     B(i,j)=(A./(dz.^4));
    end
end
%--------------------------------------------------------------------------
for i=n3+2:n4-2
    for j=i-2:i-2
    B(i,j)=A./(dz.^4);
    end
end
 for i=n3+2:n4-2
    for j=i-1:i-1
    B(i,j)=(-4.*A./(dz.^4))-C4./(dz.^2);
end
end

for i=n3+2:n4-2
    for j=i:i
     B(i,j)=(6.*A./(dz.^4))+2.*(C4./(dz.^2))+k4;
end
end
for i=n3+2:n4-2
    for j=i+1:i+1
        B(i,j)=(-4.*A./(dz.^4))-(C4./(dz.^2));
end
end

for i=n3+2:n4-2
    for j=i+2:i+2
     B(i,j)=(A./(dz.^4));
    end
end

B=B(1:n4,1:n4);

f(1,1)=((2*P0/dz)-(2*M0/(dz.^2)));f(1,2)=(M0/(dz.^2)); f(3:n4)=0;
u=inv(B)*f'
eps1(1,:)=u(1:n1,1);
eps2(1,:)=u(n1+1:n2,1);
eps3(1,:)=u(n2+1:n3,1);
eps4(1,:)=u(n3+1:n4,1); 
t1=dfdx(eps1(1,:),z1,N);
t2=dfdx(eps2(1,:),z2,N);
t3=dfdx(eps3(1,:),z3,N);
t4=dfdx(eps4(1,:),z4,N);

y1new=1./ms1.*dfdx(ms1,rt,N);
y2new=1./rt.^2.*(ms1+ms2+ms3)./ms1-1./rt.*1./ms1.*dfdx(ms3,rt,N)+ns1./ms1;
y3new=(ms2+ms3)./ms1;
y4new=1./rt.^2.*(ms1+ms2+ms3)./ms1-1./rt.*1./ms1.*dfdx(ms3,rt,N);
y5new=1./ms2.*dfdx(ms2,rt,N);
y6new=1./rt.^2.*ms1./ms2+1./rt.*1./ms2.*dfdx(ms2,rt,N)+ns2./ms2;
y7new=(ms2+ms3)./ms2;
y8new=1./rt.^2.*ms1./ms2+1./rt.*1./ms2.*dfdx(ms2,rt,N);
y1old=y1;y2old=y2;y3old=y3;y4old=y4;y5old=y5;y6old=y6;y7old=y7;y8old=y8;
my=(y1old-y1new)./y1old;

end

end

%--------------------------------------------------------------------------
function  dphidr = phiode(r,phi,rt,y1,y2,y3,y4,y5,y6,y7,y8)

x1=interp1(rt,y1,r,'cubic');
x2=interp1(rt,y2,r,'cubic');
x3=interp1(rt,y3,r,'cubic');
x4=interp1(rt,y4,r,'cubic');
x5=interp1(rt,y5,r,'cubic');
x6=interp1(rt,y6,r,'cubic');
x7=interp1(rt,y7,r,'cubic');
x8=interp1(rt,y8,r,'cubic');
  dphidr = [ 
           phi(3) 
           phi(4)
           -phi(3).*x1+x2.*phi(1)+x3./r.*phi(4)-x4.*phi(2)
           -phi(4).*x5+x6.*phi(2)-x7./r.*phi(3)-x8.*phi(1)
          
            ];
end

function res = phibc(phia,phib)
  res = [ phia(1)-1 
          phib(1)
          phia(2)-1
          phib(2)
          ];
end


 %-------------------------------------------------------------------------
function finddmdr=dfdx(mm,rt,N)
        diff1=(mm(1)-mm(2))./(rt(1)-rt(2))+(mm(1)-mm(3))./(rt(1)-rt(3))-...
            (mm(2)-mm(3))./(rt(2)-rt(3));
        diffi=(mm(1:(N-2))-mm(2:(N-1)))./(rt(1:(N-2))-rt(2:(N-1)))....
            -(mm(1:(N-2))-mm(3:N))./(rt(1:(N-2))-rt(3:N))....
            +(mm(2:(N-1))-mm(3:N))./(rt(2:(N-1))-rt(3:N));
        diffn=-(mm(N-1)-mm(N-2))./(rt(N-1)-rt(N-2))+(mm(N-2)-...
            mm(N))./(rt(N-2)-rt(N))+(mm(N-1)-mm(N))./(rt(N-1)-rt(N));
finddmdr=[diff1,diffi,diffn];
end
%--------------------------------------------------------------------------

function findG11=Gs11(z0,z1,rt,theta,phi,eps1,t1,AA,nn,G01)
epsiz1=interp1(z1,eps1(1,:),z0);
depsiz1=interp1(z1,t1,z0);
epsilon_r1=cos(theta).*epsiz1.*phi(3,:);
epsilon_theta1=1./rt.*cos(theta).*epsiz1.*(phi(1,:)+phi(2,:));
epsilon_rtheta1=-epsiz1./rt./2.*sin(theta).*(phi(1,:)+phi(2,:)-rt.*...
    phi(4,:));
epsilon_thetaz1=1./2.*sin(theta).*phi(2,:).*depsiz1;
epsilon_zr1=1./2.*cos(theta).*phi(1,:).*depsiz1;
epsilon_q1=2/3.*sqrt(1/2.*((epsilon_r1-epsilon_theta1).^2+...
    epsilon_theta1.^2+epsilon_r1.^2)....
    +3.*(epsilon_rtheta1.^2+epsilon_thetaz1.^2+epsilon_zr1.^2));
findG11=AA.*epsilon_q1.^nn;
mm=find(findG11>G01);
findG11(mm)=G01;
end

function findG12=Gs12(z11,z2,rt,theta,phi,eps2,t2,AA,nn,G02)
epsiz2=interp1(z2,eps2(1,:),z11);
depsiz2=interp1(z2,t2,z11);
epsilon_r2=cos(theta).*epsiz2.*phi(3,:);
epsilon_theta2=1./rt.*cos(theta).*epsiz2.*(phi(1,:)+phi(2,:));
epsilon_rtheta2=-epsiz2./rt./2.*sin(theta).*(phi(1,:)+phi(2,:)-...
    rt.*phi(4,:));
epsilon_thetaz2=1./2.*sin(theta).*phi(2,:).*depsiz2;
epsilon_zr2=1./2.*cos(theta).*phi(1,:).*depsiz2;
epsilon_q2=2/3.*sqrt(1/2.*((epsilon_r2-epsilon_theta2).^2+...
    epsilon_theta2.^2+epsilon_r2.^2)....
    +3.*(epsilon_rtheta2.^2+epsilon_thetaz2.^2+epsilon_zr2.^2));
findG12=AA.*epsilon_q2.^nn;
mm=find(findG12>G02);
findG12(mm)=G02;
end

function findG13=Gs13(z22,z3,rt,theta,phi,eps3,t3,AA,nn,G03)
epsiz3=interp1(z3,eps3(1,:),z22);
depsiz3=interp1(z3,t3,z22);
epsilon_r3=cos(theta).*epsiz3.*phi(3,:);
epsilon_theta3=1./rt.*cos(theta).*epsiz3.*(phi(1,:)+phi(2,:));
epsilon_rtheta3=-epsiz3./rt./2.*sin(theta).*(phi(1,:)+phi(2,:)-...
    rt.*phi(4,:));
epsilon_thetaz3=1./2.*sin(theta).*phi(2,:).*depsiz3;
epsilon_zr3=1./2.*cos(theta).*phi(1,:).*depsiz3;
epsilon_q3=2/3.*sqrt(1/2.*((epsilon_r3-epsilon_theta3).^2+...
    epsilon_theta3.^2+epsilon_r3.^2)....
    +3.*(epsilon_rtheta3.^2+epsilon_thetaz3.^2+epsilon_zr3.^2));
findG13=AA.*epsilon_q3.^nn;
mm=find(findG13>G03);
findG13(mm)=G03;
end
% 
function findG14=Gs14(z33,z4,rt,theta,phi,eps4,t4,AA,nn,G04)
epsiz4=interp1(z4,eps4(1,:),z33);
depsiz4=interp1(z4,t4,z33);
epsilon_r4=cos(theta).*epsiz4.*phi(3,:);
epsilon_theta4=1./rt.*cos(theta).*epsiz4.*(phi(1,:)+phi(2,:));
epsilon_rtheta4=-epsiz4./rt./2.*sin(theta).*(phi(1,:)+phi(2,:)-rt.*...
    phi(4,:));
epsilon_thetaz4=1./2.*sin(theta).*phi(2,:).*depsiz4;
epsilon_zr4=1./2.*cos(theta).*phi(1,:).*depsiz4;
epsilon_q4=2/3.*sqrt(1/2.*((epsilon_r4-epsilon_theta4).^2+...
    epsilon_theta4.^2+epsilon_r4.^2)....
    +3.*(epsilon_rtheta4.^2+epsilon_thetaz4.^2+epsilon_zr4.^2));
findG14=AA.*epsilon_q4.^nn;
mm=find(findG14>G04);
findG14(mm)=G04;
end
% %------------------------------------------------------------------------

function findG21=Gs21(rp,rt,theta,phi,t1,eps1,AA,nn,G01)
phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,t1,rp);
dphi_r=interp1(rt,phi(3,:),rp);
dphi_theta=interp1(rt,phi(4,:),rp);
epsilon_r1=cos(theta).*eps1(1,:).*dphi_r;
epsilon_theta1=1./rt.*cos(theta).*eps1(1,:).*(phi_r+phi_theta);
epsilon_rtheta1=-eps1(1,:)./rt./2.*sin(theta).*(phi_r+phi_theta-rt.*...
    dphi_theta);
epsilon_thetaz1=1./2.*sin(theta).*phi_theta.*t1;
epsilon_zr1=1./2.*cos(theta).*phi_r.*t1;
epsilon_q1=2/3.*sqrt(1/2.*((epsilon_r1-epsilon_theta1).^2+...
    epsilon_theta1.^2+epsilon_r1.^2)....
    +3.*(epsilon_rtheta1.^2+epsilon_thetaz1.^2+epsilon_zr1.^2));
findG21=AA.*epsilon_q1.^nn;
mm=find(findG21>G01);
findG21(mm)=G01;
end


function findG22=Gs22(rp,rt,theta,phi,t2,eps2,AA,nn,G02)
phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,t2,rp);
dphi_r=interp1(rt,phi(3,:),rp);
dphi_theta=interp1(rt,phi(4,:),rp);
epsilon_r2=cos(theta).*eps2(1,:).*dphi_r;
epsilon_theta2=1./rt.*cos(theta).*eps2(1,:).*(phi_r+phi_theta);
epsilon_rtheta2=-eps2(1,:)./rt./2.*sin(theta).*(phi_r+phi_theta-rt.*...
    dphi_theta);
epsilon_thetaz2=1./2.*sin(theta).*phi_theta.*t2;
epsilon_zr2=1./2.*cos(theta).*phi_r.*t2;
epsilon_q2=2/3.*sqrt(1/2.*((epsilon_r2-epsilon_theta2).^2+...
    epsilon_theta2.^2+epsilon_r2.^2)....
    +3.*(epsilon_rtheta2.^2+epsilon_thetaz2.^2+epsilon_zr2.^2));
findG22=AA.*epsilon_q2.^nn;
mm=find(findG22>G02);
findG22(mm)=G02;
end


function findG23=Gs23(rp,rt,theta,phi,t3,eps3,AA,nn,G03)
phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,t3,rp);
dphi_r=interp1(rt,phi(3,:),rp);
dphi_theta=interp1(rt,phi(4,:),rp);
epsilon_r3=cos(theta).*eps3(1,:).*dphi_r;
epsilon_theta3=1./rt.*cos(theta).*eps3(1,:).*(phi_r+phi_theta);
epsilon_rtheta3=-eps3(1,:)./rt./2.*sin(theta).*(phi_r+phi_theta-rt.*...
    dphi_theta);
epsilon_thetaz3=1./2.*sin(theta).*phi_theta.*t3;
epsilon_zr3=1./2.*cos(theta).*phi_r.*t3;
epsilon_q3=2/3.*sqrt(1/2.*((epsilon_r3-epsilon_theta3).^2+...
    epsilon_theta3.^2+epsilon_r3.^2)....
    +3.*(epsilon_rtheta3.^2+epsilon_thetaz3.^2+epsilon_zr3.^2));
findG23=AA.*epsilon_q3.^nn;
mm=find(findG23>G03);
findG23(mm)=G03;
end

function findG24=Gs24(rp,rt,theta,phi,t4,eps4,AA,nn,G04)
phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,t4,rp);
dphi_r=interp1(rt,phi(3,:),rp);
dphi_theta=interp1(rt,phi(4,:),rp);
epsilon_r4=cos(theta).*eps4(1,:).*dphi_r;
epsilon_theta4=1./rt.*cos(theta).*eps4(1,:).*(phi_r+phi_theta);
epsilon_rtheta4=-eps4(1,:)./rt./2.*sin(theta).*(phi_r+phi_theta-rt.*...
    dphi_theta);
epsilon_thetaz4=1./2.*sin(theta).*phi_theta.*t4;
epsilon_zr4=1./2.*cos(theta).*phi_r.*t4;
epsilon_q4=2/3.*sqrt(1/2.*((epsilon_r4-epsilon_theta4).^2+...
    epsilon_theta4.^2+epsilon_r4.^2)....
    +3.*(epsilon_rtheta4.^2+epsilon_thetaz4.^2+epsilon_zr4.^2));
findG24=AA.*epsilon_q4.^nn;
mm=find(findG24>G04);
findG24(mm)=G04;
end