function pile21

clear all;   
close all;

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

% r0=0.353;
r0=0.45;
l=16;rm=20.*r0;L=l;% nn=-0.1;
% G0=3.3e4;
G01=1.67e3;
G02=1.67e3;
G03=1.67e3;
G04=1.67e3;
M0= 0.*r0;
% Ep=16590; 
% A=922;
Ep=1.18e5;
% Ep=2.2e6; 
Ip=pi.*(r0^4)/4;nth=100;
% Ep=A/Ip;
N=20;
P0=600;
% AA=1000;
 AA=30;
nn=-0.5;
% G0=10000;
% AA=100;
% nn=-0.5;

v= 0.2;tol=0.001; iteration=0;A=Ep*Ip;dz=L/(4*N);

eps=zeros(4,N);
phi=zeros(4,N);

z0=0; z11=l/4;z22=l/2; z33=3*l/4; z44=l;






n1=N, n2=2*N; n3=3*N; n4=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);
eps5=0.000001*ones(N,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);
%--------------------------------------------------------------------------
dtheta=2*pi/nth;
theta=-dtheta/2;

m11=Ep.*Ip;
dtheta=2*pi/nth;
theta=-dtheta/2;


C1=0;
k1=0;
z0=0;
for i=1:N
    theta=theta+dtheta;
    mu1=Gs11(z0,z1,rt,theta,phi,eps1,t1,AA,nn,G01);
    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;

%coefficients for eps2 
for i=1:N
    theta=theta+dtheta;
    mu2=Gs12(z11,z2,rt,theta,phi,eps2,t2,AA,nn,G02);
    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
%coefficients for eps3 
for i=1:N
    C3=0;
k3=0;
    theta=theta+dtheta;
    mu3=Gs13(z22,z3,rt,theta,phi,eps3,t3,AA,nn,G03);
    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

for i=1:N
    C4=0;
k4=0;
    theta=theta+dtheta;
    mu4=Gs14(z33,z4,rt,theta,phi,eps4,t4,AA,nn,G04);
    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;
        mu1=Gs21(rt(i),rt,theta,phi,eps1,t1,AA,nn,G01);
        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;
        mu2=Gs22(rt(i),rt,theta,phi,eps2,t2,AA,nn,G02);
        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);
        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);
        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=ms21+ms22+ms32+ms42;
ms3=ms31+ms32+ms33+ms43;
ns1=ns11+ns21+ns31+ns41;
ns2=ns21+ns22+ns32+ns42;


B=zeros(n4,n4);
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'



B=B(1:n4,1:n4);
f(1,1)=P0; f(2: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

figure;
plot(eps1(1,:),z1,eps2(1,:),z2,eps3(1,:),z3,eps4(1,:),z4);           %plotting epsi

legend('\psi(z)');
xlabel('z');
ylabel('\psi(z)');

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);
b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;
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,t2,eps2,AA,nn,G02)

b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

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;

% findG12=G02.*(epsilon_q2./epsilon_q0).^b;
% mm=find(findG12>G0);
% findG12(mm)=G0;
%     lamda2=2.*G0.*v./(1-2.*v);
% segma_r2=lamda2.*(epsilon_r2+epsilon_theta2)+2.*G0.*epsilon_r2;
% segma_theta2=lamda2.*(epsilon_r2+epsilon_theta2)+2.*G0.*epsilon_theta2;
% segma_z2=lamda2.*(epsilon_r2+epsilon_theta2);
% segma_main2=(segma_r2+segma_theta2+segma_z2)./3;
% p2=segma_main2;
% G0=pa.*D.*((p2/pa).^n).*R0.^m; 
end

function findG13=Gs13(z22,z3,rt,theta,phi,eps3,t3,AA,nn,G03)
b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

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;

% 
% 
% findG13=G0.*(epsilon_q3./epsilon_q0).^b;
% mm=find(findG13>G0);
% findG13(mm)=G0;
%     lamda3=2.*G0.*v./(1-2.*v);
% segma_r3=lamda3.*(epsilon_r3+epsilon_theta3)+2.*G0.*epsilon_r3;
% segma_theta3=lamda3.*(epsilon_r3+epsilon_theta3)+2.*G0.*epsilon_theta3;
% segma_z3=lamda3.*(epsilon_r3+epsilon_theta3);
% segma_main3=(segma_r3+segma_theta3+segma_z3)./3;
% p3=segma_main3;
% G0=pa.*D.*((p3/pa).^n).*R0.^m; 





end
function findG14=Gs14(z33,z4,rt,theta,phi,t4,eps4,AA,nn,G04)
b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

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;



% findG14=G0.*(epsilon_q4./epsilon_q0).^b;
% mm=find(findG14>G0);
% findG14(mm)=G0;
%     lamda4=2.*G0.*v./(1-2.*v);
% segma_r4=lamda4.*(epsilon_r4+epsilon_theta4)+2.*G0.*epsilon_r4;
% segma_theta4=lamda4.*(epsilon_r4+epsilon_theta4)+2.*G0.*epsilon_theta4;
% segma_z4=lamda4.*(epsilon_r4+epsilon_theta4);
% segma_main4=(segma_r4+segma_theta4+segma_z4)./3;
% p4=segma_main4;
% G0=pa.*D.*((p4/pa).^n).*R0.^m; 



end
%------------------------------------------------------------------------------------------------------
function findG21=Gs21(rp,rt,theta,phi,t1,eps1,AA,nn,G01)

b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,phi(2,:),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;



% 
% findG21=G0.*(epsilon_q1./epsilon_q0).^b;
% mm=find(findG21>G0);
% findG21(mm)=G0;
% 
%     lamda1=2.*G0.*v./(1-2.*v);
% segma_r1=lamda1.*(epsilon_r1+epsilon_theta1)+2.*G0.*epsilon_r1;
% segma_theta1=lamda1.*(epsilon_r1+epsilon_theta1)+2.*G0.*epsilon_theta1;
% segma_z1=lamda1.*(epsilon_r1+epsilon_theta1);
% segma_main1=(segma_r1+segma_theta1+segma_z1)./3;
% p1=segma_main1;
% G0=pa.*D.*((p1/pa).^n).*R0.^m;







end
function findG22=Gs22(rp,rt,theta,phi,t2,eps2,AA,nn,G02)

b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,phi(2,:),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=-t2./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;






% findG22=G0.*(epsilon_q2./epsilon_q0).^b;
% mm=find(findG22>G0);
% findG22(mm)=G0;
% 
%     lamda2=2.*G0.*v./(1-2.*v);
% segma_r2=lamda2.*(epsilon_r2+epsilon_theta2)+2.*G0.*epsilon_r2;
% segma_theta2=lamda2.*(epsilon_r2+epsilon_theta2)+2.*G0.*epsilon_theta2;
% segma_z2=lamda2.*(epsilon_r2+epsilon_theta2);
% segma_main2=(segma_r2+segma_theta2+segma_z2)./3;
% p2=segma_main2;
% G0=pa.*D.*((p2/pa).^n).*R0.^m;
% 






end
function findG23=Gs23(rp,rt,theta,phi,t3,eps3,AA,nn,G03)


b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;




phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,phi(2,:),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;


% 
% findG23=G0.*(epsilon_q3./epsilon_q0).^b;
% mm=find(findG23>G0);
% findG23(mm)=G0;
% 
%     lamda3=2.*G0.*v./(1-2.*v);
% segma_r3=lamda3.*(epsilon_r3+epsilon_theta3)+2.*G0.*epsilon_r3;
% segma_theta3=lamda3.*(epsilon_r3+epsilon_theta3)+2.*G0.*epsilon_theta3;
% segma_z3=lamda3.*(epsilon_r3+epsilon_theta3);
% segma_main3=(segma_r3+segma_theta3+segma_z3)./3;
% p3=segma_main3;
% G0=pa.*D.*((p3/pa).^n).*R0.^m;









end
function findG24=Gs24(rp,rt,theta,phi,t4,eps4,AA,nn,G04)

b=-0.5;
epsilon_q0=10e-5; %maximum deviatoric strain of linear elastic at epsilon_q1=10e-5
pa=1;  %refrence pressure which is 1kpa 
n=0.8; m=0.2; R0=1; D=1000;v=0.2;

phi_r=interp1(rt,phi(1,:),rp);
phi_theta=interp1(rt,phi(2,:),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;



% 
% findG24=G0.*(epsilon_q4./epsilon_q0).^b;
% mm=find(findG24>G0);
% findG24(mm)=G0;
% 
%     lamda4=2.*G0.*v./(1-2.*v);
% segma_r4=lamda4.*(epsilon_r4+epsilon_theta4)+2.*G0.*epsilon_r4;
% segma_theta4=lamda4.*(epsilon_r4+epsilon_theta4)+2.*G0.*epsilon_theta4;
% segma_z4=lamda4.*(epsilon_r4+epsilon_theta4);
% segma_main4=(segma_r4+segma_theta4+segma_z4)./3;
% p4=segma_main4;
% G0=pa.*D.*((p4/pa).^n).*R0.^m;
% 


end