下面小程序计算岩石圈内的温度异常（相对于古老的稳定岩石圈，如craton而言，程序说明参见Chappell and Kusznir, 2008）：

clear
Myr=1e6*365.25*24*3600;% Unit s
a=125e3;%depth of the lithospheric plate thickness (plate model), unit m
z=0:1e3:a;% depth, unit m
gama=1.5;    % streching factor
T1=1060;    % Tm-T0 (mantle temperature - surface temperature), unit: K
tao=62.8*Myr;    % lithosphere cooling thermal decay constant
t=200*Myr;    % age of a profile considered


tmp2=0;
for n=1:30
   tmp1_1=(-1)^(n+1)/n;
   tmp1_2=gama/n/pi*sin(n*pi/gama);
   tmp1_3=exp(-n^2*t/tao)*sin(n*pi*z/a);
   tmp2=tmp1_1*tmp1_2*tmp1_3+tmp2;
end
T=T1*[2/pi*tmp2];% temperature anomaly in relative to very old lithosphere
plot(T,z/1e3)
set(gca,'ydir','reverse')
xlabel('Temperature anomaly (K)')
ylabel('Depth (km)')

结果见下图：

图中，蓝色，黑色，红色，绿色线分别代表该采样剖面的地壳年龄为0, 10, 50 和200个Ma。