clear
// Lithospheric Thermal Gravity Anomaly Correction

// Equation from Chappell and Kuznir (2008)

//Define some initial parameters
n = 0:1:20 // series in the expansion
B = 1.5 //beta factor
Tm = 1300*274.15 //base-lithosphere temperature [C]
tau = 62.8 //lithosphere cooling thermal decay constant [Ma]
a = 125000 //equilibrium lithospheric thickness [m]
alpha = 3.28e-5 //thermal expansion coefficient K^-1
// Tm needs to be converted to Kelvin - multiply by 274.15
rho = 3300
t = 0:5:150 //lithosphere age [Ma]
G = 6.67e-11

// LTGA with increasing age of oceanic lithosphere

c=8*G*alpha*rho*a*Tm/%pi
// d=2*n+1
d=(2*n)+1
// g=c*1/d^2*(B/(d*%pi))*sin(d*%pi/B)*exp(-(d^2)*t[1,1]/tau)

//      First series in expansion
g=c*(B/(%pi))*sin(%pi/B)*%e^(-t/tau)*1000

//g=c*sum(e,*)

// g=%pi*G*alpha*rho*a*Tm*(1-1/B)*exp(-t/tau)*1000

//fb = B/(n*%pi)*sin(n*%pi/B)
//        Tz=2*Tm/%pi*sum()


//  Plotting
clf()
xlabel("Age of oceanic lithosphere (Ma)");
ylabel("Lithospheric Thermal Gravity Anomaly (mGal)"); 
set(gca(),"auto_clear","off")
set(gca(),"grid",[1,1])
a=gca();
a.x_location="top";
//a.axes_reverse=["off","on","off"];
plot(t,g)

