Clear Filters
Clear Filters

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

3 views (last 30 days)
%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 Comments
Image Analyst
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?
Please read this link:
Priyank Pathak
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.
%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);

Sign in to comment.

Answers (0)

Categories

Find more on Get Started with MATLAB in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!