arctica1963

10 Reputation

3 Badges

8 years, 245 days

MaplePrimes Activity


These are replies submitted by arctica1963

Here are the four working plots and the final one I just cannot see where it is going wrong. It is like the vertical scal has gone wacky - the min/max should be around ~ -25 to +40 or so, clearly well within the plot range.

FAA_isostatic1d.mw

Here are the four working plots and the final one I just cannot see where it is going wrong. It is like the vertical scal has gone wacky - the min/max should be around ~ -25 to +40 or so, clearly well within the plot range.

FAA_isostatic1d.mw

@arctica1963 Quick update. I have added in all the other equations and generated the individual plot functions. However, only two plots worked (uncompensated and Airy), the others (Pratt, Flexure and Subloads) give empty plots and errors. The equations are all correct so it should work fine. Attached the new version.

Also, I wanted to try and display the curves on a single plot and tried the "display" function but that did not work, eg

A:=plot(f(x), etc)

B:=plot(f(k), etc)

display({A,B})

Maybe I missed something obvious?

Hopefully it will be clear from the files what is missing or in error. Thanks for the help. I have also attached the SCILAB code to show (basically a Matlab clone really).


restart

``

Te := 10000

10000

(1)

g := 9.81

9.81

(2)

upsilon := .25

.25

(3)

E := 10^11

100000000000

(4)

``

G := 6.67*10^(-11)

0.6670000000e-10

(5)

npts := 256

256

(6)

`ρload` := 2800

2800

(7)

`ρwater` := 1030

1030

(8)

`ρmantle` := 3330

3330

(9)

`ρinfill` := `ρload`

2800

(10)

`ρcrust` := `ρload`

2800

(11)

thick := 7100

7100

(12)

meandepth := 4500

4500

(13)

DC := 11370

11370

(14)

Rg := E*Te^3/(12*(1-upsilon^2))

0.8888888889e22

(15)

 

ZL := 35000

35000

(16)

XKINT := 6.259*10^(-6)

0.6259000000e-5

(17)

Zuncomp := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end proc

(18)

 

Zairy := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-k*XKINT*thick)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-k*XKINT*thick)) end proc

(19)

``

Zpratt := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-(1/2)*k*XKINT*DC)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-(1/2)*k*XKINT*DC)) end proc

(20)

``

flexp := 1/(Rg*(k*XKINT)^4/(g*(`ρmantle`-`ρinfill`))+1)

1/(0.2623749728e-2*k^4+1)

(21)

``

Zflexure := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-flexp*exp(-k*XKINT*thick)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-flexp*exp(-k*XKINT*thick)) end proc

(22)

kscale := proc (k) options operator, arrow; k*0.1e-2/XKINT end proc;

proc (k) options operator, arrow; k*0.1e-2/XKINT end proc

(23)

r1 := `ρcrust`-`ρwater`

1770

(24)

r2 := `ρmantle`-`ρcrust`

530

(25)

r3 := `ρmantle`-`ρwater`

2300

(26)

Zsubload := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(r1+r2*exp(-k*XKINT*thick)-r3*exp(-k*XKINT*ZL)/flexp) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(r1+r2*exp(-k*XKINT*thick)-r3*exp(-k*XKINT*ZL)/flexp) end proc

(27)

``

plot(`@`(Zuncomp, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = green, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Uncompensated"], legendstyle = [location = top])

 

NULL

plot(`@`(Zairy, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = red, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Airy"], legendstyle = [location = top])

 

plot(`@`(Z*pratt, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = blue, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Pratt"], legendstyle = [location = top])

 

 

plot(`@`(Zflexure, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = red, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Flexure"], legendstyle = [location = top])

 

 

plot(`@`(Zsubload, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = green, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Buried loads"], legendstyle = [location = top])

 

 

``

``

``

``

 

 

 

 

 

NULL


Download FAA_isostatic1d.mw

Scilab code:

//Example 5.1 in Watts (2001) Isostasy and Flexure of the Lithosphere

//Define input parameters
clear
Te = 10000;          //Elastic thickness
g = 9.81;            //Gravitational acceleration
v = 0.25;            //Poisson's ratio
E = 1e11;            //Elastic modulus
G = 6.67e-11;        //Gravitational constant
npts=256 // must be a power of 2

//Define densities of the load, water, mantle, infill and crust
rho_load = 2800;
rho_w = 1030;
rho_m = 3330;
rho_infill = rho_load;
rho_c = rho_load;

//Define the mean thickness of oceanic crust
Toc = 7100
//Define the mean water depth
Zw = 4500
//Define the depth of compensation
Zcomp = 113700
D = E*Te^3/(12*(1-v^2))
//Define the depth of buried loading below mean seafloor
Zl = 35000
//Choose a wavenumber interval in 1/meters for the admittance calculations
XKint = 6.259e-6

//Calculate the gravitation admittance for different isostatic models

k=linspace(0,npts/2,(npts/2)+1)
part1=2*%pi*G*1e5* (rho_load-rho_w)*%e^(-k*XKint*Zw)*1000

//Uncompensated
Zuncomp = part1

//Airy
Zairy = part1 .*(1-%e^(-k*XKint*Toc))

//Pratt
Zpratt = part1 .*(1-%e^(-k*XKint*Zcomp/2))

//Flexure
wavek=XKint*k*1000
phik = (D*(k*XKint)^4/(g*(rho_m-rho_infill))+1)^-1
Zflex = part1 .*(1-(phik .*%e^(-k*XKint*Toc)))

//Buried
r1 = (rho_c-rho_w);
r2 = (rho_m-rho_c);
r3 = (rho_m-rho_w);
phikw = (D*((k*XKint)^4)/(g*(rho_m-rho_w))+1)^-1
Zburied = (2*%pi*G*1e5*%e^(-k*XKint*Zw)).*1000 .*(r1+(r2.*%e^(-k*XKint*Toc))-((r3.*%e^(-k*XKint*Zl))./phikw))

Wavelength=2*%pi/wavek

//Plot
scf(0)
clf(0)
xlabel("Wavenumber k (1/km)");
ylabel("Free-air admittance (mGal/km)"); 
set(gca(),"auto_clear","off")
set(gca(),"grid",[1,1])
a=gca();
a.x_location="bottom";
a.axes_reverse=["off","off","off"];
//plot(k,exp(-k*XKint*Zw))
plot(wavek,Zuncomp,'g--')
plot(wavek,Zairy,'b:')
plot(wavek,Zpratt,'r')
plot(wavek,Zflex,'c')
plot(wavek,Zburied,'m--') 

//plot(k,phik,'m--')
hl=legend(['Uncompensated';'Airy';'Pratt';'Flexure';'Buried']);

@arctica1963 Quick update. I have added in all the other equations and generated the individual plot functions. However, only two plots worked (uncompensated and Airy), the others (Pratt, Flexure and Subloads) give empty plots and errors. The equations are all correct so it should work fine. Attached the new version.

Also, I wanted to try and display the curves on a single plot and tried the "display" function but that did not work, eg

A:=plot(f(x), etc)

B:=plot(f(k), etc)

display({A,B})

Maybe I missed something obvious?

Hopefully it will be clear from the files what is missing or in error. Thanks for the help. I have also attached the SCILAB code to show (basically a Matlab clone really).


restart

``

Te := 10000

10000

(1)

g := 9.81

9.81

(2)

upsilon := .25

.25

(3)

E := 10^11

100000000000

(4)

``

G := 6.67*10^(-11)

0.6670000000e-10

(5)

npts := 256

256

(6)

`ρload` := 2800

2800

(7)

`ρwater` := 1030

1030

(8)

`ρmantle` := 3330

3330

(9)

`ρinfill` := `ρload`

2800

(10)

`ρcrust` := `ρload`

2800

(11)

thick := 7100

7100

(12)

meandepth := 4500

4500

(13)

DC := 11370

11370

(14)

Rg := E*Te^3/(12*(1-upsilon^2))

0.8888888889e22

(15)

 

ZL := 35000

35000

(16)

XKINT := 6.259*10^(-6)

0.6259000000e-5

(17)

Zuncomp := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end proc

(18)

 

Zairy := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-k*XKINT*thick)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-k*XKINT*thick)) end proc

(19)

``

Zpratt := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-(1/2)*k*XKINT*DC)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-exp(-(1/2)*k*XKINT*DC)) end proc

(20)

``

flexp := 1/(Rg*(k*XKINT)^4/(g*(`ρmantle`-`ρinfill`))+1)

1/(0.2623749728e-2*k^4+1)

(21)

``

Zflexure := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-flexp*exp(-k*XKINT*thick)) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(1-flexp*exp(-k*XKINT*thick)) end proc

(22)

kscale := proc (k) options operator, arrow; k*0.1e-2/XKINT end proc;

proc (k) options operator, arrow; k*0.1e-2/XKINT end proc

(23)

r1 := `ρcrust`-`ρwater`

1770

(24)

r2 := `ρmantle`-`ρcrust`

530

(25)

r3 := `ρmantle`-`ρwater`

2300

(26)

Zsubload := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(r1+r2*exp(-k*XKINT*thick)-r3*exp(-k*XKINT*ZL)/flexp) end proc

proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth)*(r1+r2*exp(-k*XKINT*thick)-r3*exp(-k*XKINT*ZL)/flexp) end proc

(27)

``

plot(`@`(Zuncomp, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = green, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Uncompensated"], legendstyle = [location = top])

 

NULL

plot(`@`(Zairy, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = red, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Airy"], legendstyle = [location = top])

 

plot(`@`(Z*pratt, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = blue, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Pratt"], legendstyle = [location = top])

 

 

plot(`@`(Zflexure, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = red, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Flexure"], legendstyle = [location = top])

 

 

plot(`@`(Zsubload, kscale), 0 .. .8, labels = ["Wavenumber k (1/km)", "Free-air admittance - mGal/km"], labeldirections = [horizontal, vertical], color = green, linestyle = dash, view = [0 .. .8, -30 .. 80], axes = boxed, gridlines, legend = ["Buried loads"], legendstyle = [location = top])

 

 

``

``

``

``

 

 

 

 

 

NULL


Download FAA_isostatic1d.mw

Scilab code:

//Example 5.1 in Watts (2001) Isostasy and Flexure of the Lithosphere

//Define input parameters
clear
Te = 10000;          //Elastic thickness
g = 9.81;            //Gravitational acceleration
v = 0.25;            //Poisson's ratio
E = 1e11;            //Elastic modulus
G = 6.67e-11;        //Gravitational constant
npts=256 // must be a power of 2

//Define densities of the load, water, mantle, infill and crust
rho_load = 2800;
rho_w = 1030;
rho_m = 3330;
rho_infill = rho_load;
rho_c = rho_load;

//Define the mean thickness of oceanic crust
Toc = 7100
//Define the mean water depth
Zw = 4500
//Define the depth of compensation
Zcomp = 113700
D = E*Te^3/(12*(1-v^2))
//Define the depth of buried loading below mean seafloor
Zl = 35000
//Choose a wavenumber interval in 1/meters for the admittance calculations
XKint = 6.259e-6

//Calculate the gravitation admittance for different isostatic models

k=linspace(0,npts/2,(npts/2)+1)
part1=2*%pi*G*1e5* (rho_load-rho_w)*%e^(-k*XKint*Zw)*1000

//Uncompensated
Zuncomp = part1

//Airy
Zairy = part1 .*(1-%e^(-k*XKint*Toc))

//Pratt
Zpratt = part1 .*(1-%e^(-k*XKint*Zcomp/2))

//Flexure
wavek=XKint*k*1000
phik = (D*(k*XKint)^4/(g*(rho_m-rho_infill))+1)^-1
Zflex = part1 .*(1-(phik .*%e^(-k*XKint*Toc)))

//Buried
r1 = (rho_c-rho_w);
r2 = (rho_m-rho_c);
r3 = (rho_m-rho_w);
phikw = (D*((k*XKint)^4)/(g*(rho_m-rho_w))+1)^-1
Zburied = (2*%pi*G*1e5*%e^(-k*XKint*Zw)).*1000 .*(r1+(r2.*%e^(-k*XKint*Toc))-((r3.*%e^(-k*XKint*Zl))./phikw))

Wavelength=2*%pi/wavek

//Plot
scf(0)
clf(0)
xlabel("Wavenumber k (1/km)");
ylabel("Free-air admittance (mGal/km)"); 
set(gca(),"auto_clear","off")
set(gca(),"grid",[1,1])
a=gca();
a.x_location="bottom";
a.axes_reverse=["off","off","off"];
//plot(k,exp(-k*XKint*Zw))
plot(wavek,Zuncomp,'g--')
plot(wavek,Zairy,'b:')
plot(wavek,Zpratt,'r')
plot(wavek,Zflex,'c')
plot(wavek,Zburied,'m--') 

//plot(k,phik,'m--')
hl=legend(['Uncompensated';'Airy';'Pratt';'Flexure';'Buried']);

@Carl Love The plot worked great! Thanks for the guidance on this. What I think I will do now is build the remaining functions to plot and see about re-creating the SCILAB plot.

Actually I had to recode a MathCad file to SCILBAB so a bit of a convoluted workflow!

I will get back to you about how I get on.

@Carl Love The plot worked great! Thanks for the guidance on this. What I think I will do now is build the remaining functions to plot and see about re-creating the SCILAB plot.

Actually I had to recode a MathCad file to SCILBAB so a bit of a convoluted workflow!

I will get back to you about how I get on.

@Carl Love Thanks for the pointer on the restart option; now I have the curve.

Next issue is the scaling of the x-axis. As mentioned before this is the wavenumber and the scale factor here is XKINT*1000

Should a new variable be defined like wavek:=k*XKINT*1000

Thanks for the guidance on this - still learning here

Cheers

@Carl Love Thanks for the pointer on the restart option; now I have the curve.

Next issue is the scaling of the x-axis. As mentioned before this is the wavenumber and the scale factor here is XKINT*1000

Should a new variable be defined like wavek:=k*XKINT*1000

Thanks for the guidance on this - still learning here

Cheers

Hello,

I tried the command you showed but it did not work (after a restart to clear the variables).

One thing I did spot was that the x-axis should be k*XKINT which represents the wavenumber; Zuncomp is the free-air admittance gravity anomaly.

I verified that the results of the function were correct using Zuncom(0) etc so that is working fine.

The equation is a conversion from MathCad, already have a version runninbg fine in Scilab, so all the logic is fine. The only difference is that k is incremented by 1 to npts/2, and I cannot see how that works in Maple.

Atached a revised Maple file to illustrate:

Cheers

Lester


Te := 10000

10000

(1)

g := 9.81

9.81

(2)

upsilon := .25

.25

(3)

E := 10^11

100000000000

(4)

G := 6.67*10^(-11)

0.6670000000e-10

(5)

npts := 1024

1024

(6)

`ρload` := 2800

2800

(7)

`ρwater` := 1030

1030

(8)

`ρmantle` := 3330

3330

(9)

`ρinfill` := `ρload`

2800

(10)

`ρcrust` := `ρload`

2800

(11)

thick := 7100

7100

(12)

meandepth := 4500

4500

(13)

DC := 11370

11370

(14)

Rg := E*Te^3/(12*(1-upsilon^2))

0.8888888889e22

(15)

 

ZL := 35000

35000

(16)

XKINT := 6.259*10^(-6)

0.6259000000e-5

(17)

``

Zuncomp := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end procNULL``

``

NULL

NULL

Zuncomp(10)(->)55.971

Zuncomp(0)

23.61180000*Pi

(18)

(->)

74.179

(19)

restart

plot(Zuncomp(k), k = 0 .. (1/2)*npts, labels = ["k", "Zuncomp"])

 

The Zuncomp range should be of the order 0 to 80 on the graph, with the x-axis as a function of k*XKINT*1000 (=wavenumber k) The plot below from Scilab is what I am aiming at, with the green curve the function being tested in Maple:

 

NULL

``


Download FAA_isostatic1a.mw

Hello,

I tried the command you showed but it did not work (after a restart to clear the variables).

One thing I did spot was that the x-axis should be k*XKINT which represents the wavenumber; Zuncomp is the free-air admittance gravity anomaly.

I verified that the results of the function were correct using Zuncom(0) etc so that is working fine.

The equation is a conversion from MathCad, already have a version runninbg fine in Scilab, so all the logic is fine. The only difference is that k is incremented by 1 to npts/2, and I cannot see how that works in Maple.

Atached a revised Maple file to illustrate:

Cheers

Lester


Te := 10000

10000

(1)

g := 9.81

9.81

(2)

upsilon := .25

.25

(3)

E := 10^11

100000000000

(4)

G := 6.67*10^(-11)

0.6670000000e-10

(5)

npts := 1024

1024

(6)

`ρload` := 2800

2800

(7)

`ρwater` := 1030

1030

(8)

`ρmantle` := 3330

3330

(9)

`ρinfill` := `ρload`

2800

(10)

`ρcrust` := `ρload`

2800

(11)

thick := 7100

7100

(12)

meandepth := 4500

4500

(13)

DC := 11370

11370

(14)

Rg := E*Te^3/(12*(1-upsilon^2))

0.8888888889e22

(15)

 

ZL := 35000

35000

(16)

XKINT := 6.259*10^(-6)

0.6259000000e-5

(17)

``

Zuncomp := proc (k) options operator, arrow; 200000000*Pi*G*(`ρload`-`ρwater`)*exp(-k*XKINT*meandepth) end procNULL``

``

NULL

NULL

Zuncomp(10)(->)55.971

Zuncomp(0)

23.61180000*Pi

(18)

(->)

74.179

(19)

restart

plot(Zuncomp(k), k = 0 .. (1/2)*npts, labels = ["k", "Zuncomp"])

 

The Zuncomp range should be of the order 0 to 80 on the graph, with the x-axis as a function of k*XKINT*1000 (=wavenumber k) The plot below from Scilab is what I am aiming at, with the green curve the function being tested in Maple:

 

NULL

``


Download FAA_isostatic1a.mw

Page 1 of 1