# my code is not running , i have attached my code.

3 views (last 30 days)
Priyank Pathak on 16 Sep 2021
Edited: darova on 17 Sep 2021
%In image, there are two equations. I put equation (3) {Zc} in equation (5) and get quadratic equation of ZL . I developed a code for solution of quadratic equation , in which i used 'if ' condition {for eq. 3 ,where E>0 rho _w =0) .Please look at it and correct it. Result is not getting ,what it is expected,may be something is wrong.
% syms N E real
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3300; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8; % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere and crust
rmc=(r_mc.^2);
r_cm=r_c-r_m;
rma=(r_ma.^2);
ur_ac=1/(r_ac);
ur_mc=1/(r_mc);
rcw=(r_cw.^2);
ZL= 120000; %in metres
ZC= 27000; % in metres
Zmax2=Zmax^2;
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
C2=3.14*6.67*10^(-11);
N0=-(C2*(E2*r_w+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;
%r_cw=r_c;
%zC=(r_a*L0+E*(r_cw+zL*r_ma)*ur_mc;
%A=(r_c*rma/rmc)+r_m-(r_m*rma/rmc)-r_a;
for i= 1:292681;
if E(i)>0
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_c+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*r_c*r_c/rmc));
D2(i)=(2*r_a*r_c*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;
N23(i)=N22(i)';
c11(i)=D5(i)-N23(i);
xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
%xL1=real(xL11);
xL1(i)=xL11(i)/1000;
xL33(i)=xL1(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
xL2(i)=xL22(i)/1000;
xL44(i)=xL2(i);
else
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_cw+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*rcw/rmc));
D2(i)=(2*r_a*r_cw*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;
N23(i)=N22(i)';
c11(i)=D5(i)-N23(i);
xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
%xL1=real(xL11);
xL1(i)=xL11(i)/1000;
xL33(i)=xL1(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
xL2(i)=xL22(i)/1000;
xL44(i)=xL2(i);
end
end;
% xL11=(-(B)+((B.^2)-4*A*c11).^(1/2))/(A2);
% %xL1=real(xL11);
% xL1=xL11;
%
% xL33=xL1/1000;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xL22=(-(B)-((B.^2)-4*A*c11).^(1/2))/(A2);
% % xL2=real(xL22);
% xL2=xL22/1000;
%
% xL44=xL2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zL2=reshape(xL33,541,541);
zL=zL2';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL334=xL33';
zC2=(r_a*L0+E*r_c+xL334*r_ma)/(r_mc);
zC3=reshape(zC2,541,541);
zC=zC3';
subplot(2,1,1)
im= imagesc([117 135], [-38 -20] ,zL);
colorbar
hold on
colorbar;
%brighten('h',);
%C=contour(A,'LineColor','black');
x=117:0.033333333:135;
y=-38:0.033333333:-20;
minx = round(min(zL(:)),1);
maxx = round(max(zL(:)),1);
levels = minx:5:maxx;
levels2=round(levels);
[X,Y] = meshgrid(x,y);
%[C,h]=contour(X,Y,csia21,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
[C,h]=contour(X,Y,zL,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);
subplot(2,1,2)
%jet map
im= imagesc([117 135], [-38 -20] ,zC);
colorbar
hold on
colorbar;
%brighten('h',);
%C=contour(A,'LineColor','black');
x=117:0.033333333:135;
y=-38:0.033333333:-20;
minx = round(min(zC(:)),1);
maxx = round(max(zC(:)),1);
levels = minx:5:maxx;
levels2=round(levels);
[X,Y] = meshgrid(x,y);
%[C,h]=contour(X,Y,csia21,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
[C,h]=contour(X,Y,zC,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);
Image Analyst on 16 Sep 2021
1. Can you attach your data here with the paperclip icon? Or is it too big (more than 5 MB)?
2. Can you format your code as code?
3. Can you describe what is wrong. I imagine that I will run it and it will run and produce some numbers but that you don't like them for some unspecified reason.
4. What is the expected answer?
Priyank Pathak on 16 Sep 2021
Edited: darova on 17 Sep 2021
@Image Analyst,data is more than 5MB. That's why i attached a link.
% i attached an image of required result. left one for zC and right one for zL.
% not getting required result.
% syms N E real
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3300; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8; % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere and crust
rmc=(r_mc.^2);
r_cm=r_c-r_m;
rma=(r_ma.^2);
ur_ac=1/(r_ac);
ur_mc=1/(r_mc);
rcw=(r_cw.^2);
ZL= 120000; %in metres
ZC= 27000; % in metres
Zmax2=Zmax^2;
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
C2=3.14*6.67*10^(-11);
N0=-(C2*(E2*r_w+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;
for i= 1:292681;
if E(i)>0
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_c+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*r_c*r_c/rmc));
D2(i)=(2*r_a*r_c*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;
N23(i)=N22(i)';
c11(i)=D5(i)-N23(i);
xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
%xL1=real(xL11);
xL1(i)=xL11(i)/1000;
xL33(i)=xL1(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
xL2(i)=xL22(i)/1000;
xL44(i)=xL2(i);
else
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_cw+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*rcw/rmc));
D2(i)=(2*r_a*r_cw*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;
N23(i)=N22(i)';
c11(i)=D5(i)-N23(i);
xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
%xL1=real(xL11);
xL1(i)=xL11(i)/1000;
xL33(i)=xL1(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
xL2(i)=xL22(i)/1000;
xL44(i)=xL2(i);
end
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zL2=reshape(xL33,541,541);
zL=zL2';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xL334=xL33';
zC2=(r_a*L0+E*r_c+xL334*r_ma)/(r_mc);
zC3=reshape(zC2,541,541);
zC=zC3';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1)
im= imagesc([117 135], [-38 -20] ,zL);
colorbar
hold on
colorbar;
x=117:0.033333333:135;
y=-38:0.033333333:-20;
minx = round(min(zL(:)),1);
maxx = round(max(zL(:)),1);
levels = minx:5:maxx;
levels2=round(levels);
[X,Y] = meshgrid(x,y);
[C,h]=contour(X,Y,zL,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2)
im= imagesc([117 135], [-38 -20] ,zC);
colorbar
hold on
colorbar;
x=117:0.033333333:135;
y=-38:0.033333333:-20;
minx = round(min(zC(:)),1);
maxx = round(max(zC(:)),1);
levels = minx:5:maxx;
levels2=round(levels);
[X,Y] = meshgrid(x,y);
[C,h]=contour(X,Y,zC,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
clabel(C,h,'FontSize',6);