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
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:
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 CenterFile Exchange에서 Get Started with MATLAB에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by