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).

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).

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

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

Download FAA_isostatic1a.mw