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

2 次查看（过去 30 天）
Priyank Pathak
Priyank Pathak 2021-9-16
编辑： darova 2021-9-17
%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.
%thanks in Advance
%link for data
% https://drive.google.com/drive/folders/1fmQHWyp_cvBs9ULuWUpL3uCuSAykJbZl?usp=sharing
E=xlsread('toplatlongecs','c1:c292681');
N=xlsread('geoidlatlongecs','c1:c292681');
% syms N E real
L0=2400; %load in metres
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);
  2 个评论
Image Analyst
Image Analyst 2021-9-16
  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?
Please read this link:
TUTORIAL: How to ask a question (on Answers) and get a fast answer
Priyank Pathak
Priyank Pathak 2021-9-16
编辑：darova 2021-9-17
@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.
%link for data: LINK
% not getting required result.
E=xlsread('toplatlongecs','c1:c292681');
N=xlsread('geoidlatlongecs','c1:c292681');
% syms N E real
L0=2400; %load in metres
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);

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 2-D and 3-D Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by