% gdama02
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Pythion, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-02, Edit and integration with text
% gdama02_01
%
% straight line fit to 1965-2016 global temperature data
 
% Hansen, J., Mki. Sato, R. Ruedy, K. Lo, D.W. Lea, and M. Medina-Elizade,
% 2006: Global temperature change. Proc. Natl. Acad. Sci., 103, 14288-14293,
% doi:10.1073/pnas.0606291103.
 
clear all;
fprintf('gdama02_01\n');
gdama02_01
 
% load data from text file
D=load('../data/global_temp.txt');
t=D(:,1);
d=D(:,2);
N=length(d);
 
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1965, 2020, -0.5, 1.0] );
plot(t,d,'r-','LineWidth',3);
plot(t,d,'ro','LineWidth',3);
xlabel('calendar year');
ylabel('temperature anomaly, deg C');
 
% set up data kernel
M=2;
G=[ones(N,1), t];
 
% least squares solution and predicted data
mest = (G'*G)\(G'*d);
dpre = G*mest;
 
% plot straight line fit
plot(t,dpre,'b-','LineWidth',3);
fprintf('Caption: Straight line fit (blue line) to temperature anomaly data (red curve)\n');
Caption: Straight line fit (blue line) to temperature anomaly data (red curve)
 
% variance calculation
e = (d - dpre);
E = e'*e;
s2d = E/N;
C = s2d * inv(G'*G);
sm = sqrt(C(2,2));
 
% display slope
fprintf('slope: %f +/- %f deg/yr\n', mest(2), 2*sm);
slope: 0.017534 +/- 0.001677 deg/yr
% gda01_02
%
% setup for parabola
% Supports Section 1.3.2
 
% Hansen, J., Mki. Sato, R. Ruedy, K. Lo, D.W. Lea, and M. Medina-Elizade,
% 2006: Global temperature change. Proc. Natl. Acad. Sci., 103, 14288-14293,
% doi:10.1073/pnas.0606291103.
 
clear all;
 
D=load('../data/global_temp.txt');
to=D(:,1);
d=D(:,2);
N=length(d);
 
% set up data kernel
% note time shift of 1965, necessary to avoid squaring very big numbers
t = to - 1965;
M=3;
G=[ones(N,1), t, t.^2];
 
% write out some of G
fprintf('First few rows of data kernel G\n');
First few rows of data kernel G
G(1:10,:)
ans = 10×3
1 0 0 1 1 1 1 2 4 1 3 9 1 4 16 1 5 25 1 6 36 1 7 49 1 8 64 1 9 81
 
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1965, 2020, -0.5, 1.0] );
plot(to,d,'r-','LineWidth',3);
plot(to,d,'ro','LineWidth',3);
xlabel('calendar year');
ylabel('temperature anomaly, deg C');
 
% least squares solution and predicted data
mest = (G'*G)\(G'*d);
dpre = G*mest;
 
% plot parabolic fit
plot(to,dpre,'b-','LineWidth',3);
fprintf('Caption: Quadratic fit (blue line) to temperature anomaly data (red curve)\n');
Caption: Quadratic fit (blue line) to temperature anomaly data (red curve)
 
% variance calculation
e = (d - dpre);
E = e'*e;
s2d = E/N;
C = s2d * inv(G'*G);
sm2 = diag(C);
sm = sqrt(sm2);
 
% display slope
fprintf('coefficients of equation d = m1 + m2*(t-1965) + m3*(t-1965)^2');
coefficients of equation d = m1 + m2*(t-1965) + m3*(t-1965)^2
fprintf('m1: %f +/- %f (95%%)\n', mest(1), 2*sm(1) );
m1: -0.084579 +/- 0.072183 (95%)
fprintf('m2: %f +/- %f (95%%)\n', mest(2), 2*sm(2) );
m2: 0.014783 +/- 0.006545 (95%)
fprintf('m3: %f +/- %f (95%%)\n', mest(3), 2*sm(3) );
m3: 0.000054 +/- 0.000124 (95%)
% gdama01_03
%
% G for various cases
 
clear all;
fprintf('gdama01_03\n');
gdama01_03
 
 
clear all;
 
z = [1, 2, 4, 8, 9, 12, 15, 20]';
N=length(z);
 
% straight line case
M=2;
G=zeros(N,M);
G=[ones(N,1), z];
fprintf('Kernel for linear case\n');
Kernel for linear case
G
G = 8×2
1 1 1 2 1 4 1 8 1 9 1 12 1 15 1 20
 
% quadratic case
M=2;
G=zeros(N,M);
G=[ones(N,1), z, z.^2];
fprintf('Kernel for quadratic case\n');
Kernel for quadratic case
G
G = 8×3
1 1 1 1 2 4 1 4 16 1 8 64 1 9 81 1 12 144 1 15 225 1 20 400
 
% acoustic tomography case
N=8;
M=16;
G=zeros(N,M);
h=1;
for i = (1:4)
for j = (1:4)
% measurements over rows
k = (i-1)*4 + j;
G(i,k)=h;
% measurements over columns
k = (j-1)*4 + i;
G(i+4,k)=h;
end
end
fprintf('Kernel for tomography case\n');
Kernel for tomography case
G
G = 8×16
1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1
 
% CAT scan case
% let the image be square, from (0 to 1) in both x and y
% let it be divided into M=KxL pixels which constitute model parameters
% the pixel array is unwrapped using k=i+(j-1)*K
% let N rays be collected with endpoits at (x1,y1) and (x2,y2).
K=10;
L=10;
M=K*L;
N=3;
Dx=1/K; % x dimension of pixel
Dy=1/L; % y dimension of pizel
% endpoints of rays
X1 = [ [0,0]', [0,1]', [0,0.5]' ];
X2 = [ [1,1]', [1,0]', [1, 0.5]' ];
Ds = 0.001;
I = zeros(K,L); % image of rays, for plotting purposes
G = zeros( N, M );
for n=(1:N)
% step along ray
x1=X1(:,n); x2=X2(:,n); % end points of ray
R = sqrt((x2-x1)'*(x2-x1)); % length of ray
t = (x2-x1)/R; % tangent to ray
Nr = floor( R/Ds ); % number of ray segments
x=x1; % starting point along ray
for ir = (1:Nr)
x = x1 + (ir-1)*Ds*t; % point along ray
% indices of pixel containing the point
px = floor( x(1) / Dx );
if( px<1 ) % check that is withing bounds
px=1;
elseif( px>K )
px=K;
end
py = floor( x(2) / Dy );
if( py<1 ) % check that is withing bounds
py=1;
elseif( py>L )
py=L;
end
% index in unwrapped vector of model parameters
k=px+(py-1)*K;
I(px,py) = I(px,py)+Ds; % make inage of rays
G(n,k)=G(n,k)+Ds; % increment data kernel
end
end
 
gda_draw(I,'caption Rays');
fprintf('Caption: Ray diagram\n');
Caption: Ray diagram
fprintf('\n');
fprintf('Transposed data kernel, G-transpose\n');
Transposed data kernel, G-transpose
full(G)'
ans = 100×3
0.2830 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1410 0 0 0.1410 0 0 0 0
% gdama01_04
%
% model Mars rover Mossbauer spectra using Lorentzian curves
 
clear all;
fprintf('gdama01_04\n');
gdama01_04
 
% load data
D=load('../data/mars_soil.txt');
v=D(:,1);
d=D(:,2);
d=d/max(d); % normalize
N=length(d);
 
% delete negative velocities
i=find(v>=0,1);
v=v(i:N);
d=d(i:N);
N=length(v);
 
% number of peaks determined by clicking on graph
disp('click on the bottom each peak');
click on the bottom each peak
disp(' click to the left of zero when done');
click to the left of zero when done
disp(' ');
 
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, 12, min(d), max(d)] );
plot(v,d,'r-','LineWidth',2);
plot(v,d,'ro','LineWidth',3);
xlabel('velocity, mm/s');
ylabel('counts');
title('click bottom of each peak, then left of axis');
 
% lorentzian curve of peak amplitude a, center velocity v0 and width c
% f(v) = a c^2 / ( (v-v0)^2 + c^2) )
% df/da = c^2 / ( (v-v0)^2 + c^2) )
% df/dv0 = 2 a c^2 (v-v0) / ( (v-v0)^2 + c^2) )^2
% df/dc = 2 a c / ( (v-v0)^2 + c^2) ) - a c^3 / ( (v-v0)^2 + c^2) )^2
% 3 model parameters per lorentzian, (a, v0, c)
 
% estimate of backgroind level
A = max(d);
 
% input peaks
MAXPEAKS=100;
a = zeros(MAXPEAKS,1);
v0 = zeros(MAXPEAKS,1);
c = zeros(MAXPEAKS,1);
K=0;
for k = (1:20)
p = ginput(1);
if( p(1) < 0 )
break;
end
K=K+1;
a(K) = p(2)-A;
v0(K)=p(1);
c(K)=0.1;
end
a = a(1:K);
v0 = v0(1:K);
c = c(1:K);
 
% model parameters
M=K*3;
m = [a', v0', c']';
 
for iter=(1:10)
 
dpre = A*ones(N,1);
for i = (1:K)
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
dpre = dpre + m(i)*(m(2*K+i)^2)./temp;
end
 
% data kernel
G=zeros(N,M);
for i = (1:K)
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
G(:,i) = (m(2*K+i)^2)./temp; % d/da
G(:,K+i) = 2*m(i)*(m(2*K+i)^2)*(v-m(K+i))./(temp.^2); % d/dv0
G(:,2*K+i) = 2*m(i)*m(2*K+i)./temp - 2*m(i)*m(2*K+i)^3./(temp.^2); % d/dc
end
 
% deviations in data
dd = d - dpre;
E = dd'*dd;
 
% updated model
dm = (G'*G+0.001*eye(M))\(G'*dd);
mold=m;
m = m+dm;
 
fprintf('%f\n', E);
 
end
0.338978 0.106137 0.041103 0.038910 0.038596 0.038525 0.038508 0.038504 0.038503 0.038502
 
% plot
dpre2 = A*ones(N,1);
for i = (1:K)
temp = ((v-m(K+i)).^2+m(2*K+i)^2);
dpre2 = dpre2 + m(i)*(m(2*K+i)^2)./temp;
end
plot(v,dpre2,'b-','LineWidth',3);
clear all;
fprintf('Caption: Mossbauer spectral data (red dots) and best fit (blue curve)\n');
Caption: Mossbauer spectral data (red dots) and best fit (blue curve)
% gdama01_05
 
% plot for mixing example figure
 
clear all;
fprintf('gdama01_05\n');
gdama01_05
 
M=5; % number of elements
s1 = [1, 0, 2, 0.5, 3]'; % endmember (factor) 1
s2 = [3, 0.7, 1, 2.5, 2]'; % endmember (factor) 2
 
% sample 1
x=0.1;
sx1 = x*s2 + (1-x)*s1;
 
% sample 2
x=0.3;
sx2 = x*s2 + (1-x)*s1;
 
% sample 3
x=0.5;
sx3 = x*s2 + (1-x)*s1;
 
% sample 4
x=0.7;
sx4 = x*s2 + (1-x)*s1;
 
% sample 5
x=0.9;
sx5 = x*s2 + (1-x)*s1;
 
% make some simple vector plots
gda_draw(' ', s1, ' ', sx1, ' ', sx2, ' ', sx3, ' ', sx4, ' ', sx5, ' ', s2);
fprintf('Caption: Portrayal of linear mixing between two end members (left and right),\n');
Caption: Portrayal of linear mixing between two end members (left and right),
fprintf('each with five elements\n');
each with five elements
 
% gdama01_06
%
% filter example, using seismometer response as an exemplary filter
%
% Note: this routine uses fairly advanced techniques to
% construct a seismometer response filter that are not
% going to be explained in sufficient detail for a reader
% to follow. Sorry about that ... Note that the response
% (which is just a time series) is stored in a the vector
% uoutf1() and in the file ../data/seismometer_response.txt.
 
clear all;
fprintf('gdama01_06\n');
gdama01_06
 
% response of the Broadband East component of station J59A of
% of the US Transportable array is in this file (in inscruitable
% IRIS RESP format).
respfile1='../data/TA_J59A_BHE.txt';
 
% read the IRIS RESP file for this seismometer
[ Nzeroes, zero, Npoles, pole, F, status ] = gda_read_resp( respfile1 );
 
% set up an input ground motion with a single small spike
N=100;
Dt = 0.025;
t = Dt*[0:N-1]';
uin = zeros(N,1);
uin(1) = (1e-6);
 
% apply station 1 response, filter result, and plot
[ uout1, status ] = gda_apply_response( uin, Dt, respfile1 );
flow = 0.005;
fhigh = 1.0;
uoutf1 = gda_chebyshevfilt(uout1, Dt, flow, fhigh);
 
% save the response filter in a text file
dlmwrite('../data/seismometer_response.txt', [t, uoutf1], 'delimiter', '\t');
 
% plot a spike in ground motion
figure();
clf;
subplot(2,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, Dt*(N-1), -abs(max(uin)), abs(max(uin)) ] );
plot( Dt*[0:N-1]', uin, 'k-', 'LineWidth', 3 );
title('model parameters m');
xlabel('time (s)');
ylabel('velocity, microns');
 
% plot the response of the seismometer
subplot(2,1,2);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
axis( [0, Dt*(N-1), -max(abs(uoutf1)), max(abs(uoutf1)) ] );
plot( Dt*[0:N-1]', uoutf1, 'k-', 'LineWidth', 3 );
title('response 1');
xlabel('time (s)');
ylabel('digital counts');
fprintf('(caption: (top) Ground velocity is an impuslive time series (black curve).\n');
(caption: (top) Ground velocity is an impuslive time series (black curve).
fprintf('(right) Seismometer response is a broad pulse\n');
(right) Seismometer response is a broad pulse
 
% set up an input seismogram with several small spikes
% and a smoothly varying function
N=1024;
Dt = 0.025;
t = Dt*[0:N-1]';
uin = zeros(N,1);
uin(floor(N/4)) = (1e-6);
uin(floor(8+N/2)) = 0.5*(1e-6);
uin(floor(17+N/2)) = 0.5*(1e-6);
uin(750:850) = 0.5*(1e-6)*sin( pi*([750:850]'-750)/100 );
 
% apply the response
[ uout1, status ] = gda_apply_response( uin, Dt, respfile1 );
flow = 0.005;
fhigh = 1.0;
uoutf1 = gda_chebyshevfilt(uout1, Dt, flow, fhigh);
 
% plot the complicated ground motion
figure();
clf;
subplot(2,1,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, Dt*(N-1), -max(abs(uin)), max(abs(uin)) ] );
plot( Dt*[0:N-1]', uin, 'r-', 'LineWidth', 3 );
title('model parameters m');
xlabel('time (s)');
ylabel('velocity, microns');
 
% plot the corresponding seismometer response
subplot(2,1,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [0, Dt*(N-1), -max(abs(uoutf1)), max(abs(uoutf1)) ] );
plot( Dt*[0:N-1]', uoutf1, 'b-', 'LineWidth', 3 );
title('Seismometer recording');
xlabel('time (s)');
ylabel('digital counts');
fprintf('Caption: (top) Ground velocity. (Bottom) Seismometer recording\n');
Caption: (top) Ground velocity. (Bottom) Seismometer recording
fprintf('of ground velocity.\n');
of ground velocity.
% gdama01_07
% plot example of filter, both as wiggles and a matrix
 
clear all;
fprintf('gdama01_07\n');
gdama01_07
 
% load exemplary filter (for a seismometer)
D = load('../Data/seismometer_response.txt');
t = D(:,1);
g = D(:,2);
N = length(t);
Dt=t(2)-t(1);
 
figure();
clf;
 
% build data kernel corresponding to convolution by this filter
G = toeplitz( g, [g(1), zeros(1,N-1)] );
 
% plot filter, as wiggles
clist = [1:17:N]';
Nc = length(clist);
for i=(1:Nc)
subplot(1,Nc+2,i);
hold on;
set(gca,'LineWidth',3);
set(gca,'XAxisLocation','top');
axis ij;
axis( [-max(abs(g)), max(abs(g)), t(1), t(end)] );
plot(G(:,clist(i)),t,'k-','LineWidth',3);
xlabel(sprintf('g(t-%.1f)',(clist(i)-1)*Dt));
if( i==1 )
ylabel('time (s)');
end
end
 
% build model parameters
m = zeros(N,1);
m(1:floor(N/5)) = 0.5*1e-6*sin( pi*([1:floor(N/5)]'-1)/(N/5-1));
 
% plot input (model parameters, ground motion)
subplot(1,Nc+2,Nc+1);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
set(gca,'XAxisLocation','top');
axis ij;
axis( [-max(abs(m)), max(abs(m)), t(1), t(end)] );
plot(m,t,'b-','LineWidth',3);
xlabel('m');
 
d = G*m; % perform convolution
 
% plot output (data, seismometer response)
subplot(1,Nc+2,Nc+2);
hold on;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
set(gca,'XAxisLocation','top');
axis ij;
axis( [-max(abs(d)), max(abs(d)), t(1), t(end)] );
plot(d,t,'r-','LineWidth',3);
xlabel('d');
 
fprintf('Caption: Graphical depiction of convolution process. (left 6 plots)\n');
Caption: Graphical depiction of convolution process. (left 6 plots)
fprintf('(7th plot) Ground veocity. (8th plot) Seismometer recording.\n');
(7th plot) Ground veocity. (8th plot) Seismometer recording.
 
% plot as matrix
gda_draw(G, 'caption G', m, 'caption m', '=', d, 'caption d');
fprintf('Caption: Graphical depiction on the equation G m = d, in the\n');
Caption: Graphical depiction on the equation G m = d, in the
fprintf('case where G m represents a convolution of m with a seismometer\n');
case where G m represents a convolution of m with a seismometer
fprintf('response g\n');
response g
% gdama01_08
 
% illustrations of different types of solutions
 
clear all;
fprintf('gdama01_08\n');
gdama01_08
 
 
M = 2;
Dm = 1;
Nm = 101;
m = Dm*[0:Nm-1]';
mmin = m(1);
mmax = m(end);
mbar = [40, 60]';
sm = [5, 10];
sm2 = sm.^2;
 
figure();
clf;
 
% solution is estimate and error bars
subplot(1,3,1);
hold on;
set(gca,'LineWidth',3);
set(gca,'Fontsize',14);
axis( [mmin, mmax, mmin, mmax] );
plot( [mbar(1), mbar(1)], [mbar(2)-2*sm(2), mbar(2)+2*sm(2)], 'r-', 'LineWidth', 3 );
plot( [mbar(1)-2*sm(1), mbar(1)+2*sm(1)], [mbar(2), mbar(2)], 'r-', 'LineWidth', 3 );
plot( mbar(1), mbar(2), 'ko', 'LineWidth', 3 );
xlabel('m_1');
ylabel('m_2');
 
% solution is pdf
% use tensor product, m1 increases along rows, m2 along columns
p = exp(-0.5*(m-mbar(1)).^2/sm2(1)) * exp(-0.5*(m-mbar(2)).^2/sm2(2))';
pnorm = (Dm^2)*sum(p(:));
p = p/pnorm;
subplot(1,3,2);
hold on;
set(gca,'LineWidth',3);
set(gca,'Fontsize',14);
axis( [mmin, mmax, mmin, mmax] );
colormap('jet');
% trick for checking that orientation is right
% put a defect at m1=20, m2=50 and see if it plots correctly
% p(20,40) = max(p(:));
imagesc([0,100],[0,100],p');
xlabel('m_1');
ylabel('m_2');
colorbar;
 
N = 100;
m1real = (sm(1)*randn(N,1)+mbar(1));
m2real = (sm(2)*randn(N,1)+mbar(2));
subplot(1,3,3);
hold on;
set(gca,'LineWidth',3);
set(gca,'Fontsize',14);
axis( [mmin, mmax, mmin, mmax] );
plot( m1real, m2real, 'k.', 'LineWidth', 5 );
xlabel('m_1');
ylabel('m_2');
fprintf('Caption: Illustration of different kinds of solutions to an inverse\n');
Caption: Illustration of different kinds of solutions to an inverse
fprintf('problem. (left) Estimate (circle) and 95 pecent confidence intervals (red)\n');
problem. (left) Estimate (circle) and 95 pecent confidence intervals (red)
fprintf('(middle) Probability density function (colors). (right) Ensemble of probable\n');
(middle) Probability density function (colors). (right) Ensemble of probable
fprintf('solutions (dots).');
solutions (dots).
 
% gdama01_09
%
% range of pdfs, from simple to complicated
 
clear all;
fprintf('gdama01_09\n');
gdama01_09
 
% p axes
Dm = 0.1;
N = 101;
m1 = Dm*[0:N-1]';
mmin=0;
mmax=10;
 
% unimodal distribution
sm = 1.0;
m1bar = 5.0;
p1a = exp(-0.5*(m1-m1bar).^2/(sm^2));
norm = Dm*sum(p1a);
p1a = p1a/norm;
 
% bimodal distribution
sma = 0.7;
m1bara = 3.0;
smb = 1.4;
m1barb = 8.0;
p1b = exp(-0.5*(m1-m1bara).^2/(sma^2))+0.4*exp(-0.5*(m1-m1barb).^2/(sma^2));
norm = Dm*sum(p1b);
p1b = p1b/norm;
 
% absurdly complicated distribution
p1c = zeros(N,1);
for i=(1:10)
Ac=random('Uniform',0,1);
m1barc=random('Uniform',0,10);
smc=random('Uniform',0.05,1);
tmp=Ac*exp(-0.5*(m1-m1barc).^2/(smc^2));
p1c = p1c + tmp;
end
norm = Dm*sum(p1c);
p1c = p1c/norm;
 
figure();
set(gcf,'pos',[10, 10, 800, 250]);
clf;
 
% plot unimodal distribution
subplot(1,3,1);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, 0.5 ] );
plot(m1,p1a,'b-','LineWidth',3);
xlabel('m');
ylabel('p(m)');
 
% plot bimodal distribution
subplot(1,3,2);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, 0.5 ] );
plot(m1,p1b,'b-','LineWidth',3);
xlabel('m');
ylabel('p(m)');
 
% plot complicated distribution
subplot(1,3,3);
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [mmin, mmax, 0, 0.5 ] );
plot(m1,p1c,'b-','LineWidth',3);
xlabel('m');
ylabel('p(m)');
 
fprintf('Caption: (Left) Unimodal pdf. (Middle) Bimodal pdf. (right) Complicated pdf\n.');
Caption: (Left) Unimodal pdf. (Middle) Bimodal pdf. (right) Complicated pdf .