% gdama01
% W. Menke (2024), Geophysical Data Analysis and Inverse Theory
% with MATLAB(R) or Python, 5th Edition
% Code History:
% 2022-10-01, Imported from 4th Edition by W. Menke
% 2023-04-02, Edit and integration with text
% gdama01_01
% displays the current date
date
ans = '02-Apr-2023'
% gdama01_02
% examples of directory manipulation, support Section I.3
% print working directory
pwd
ans = 'D:\LiveScripts'
 
% change directory to one up
cd ..
pwd
ans = 'D:\'
 
% change directory back to LiveScripts
cd LiveScripts
pwd
ans = 'D:\LiveScripts'
 
% list files and sub-folders in the current directory
dir
. gda_read_resp.m .. gda_timedelay.m OLD gdama_01.mlx gda_Esurface.m gdama_02.mlx gda_FTFmul.m gdama_03.mlx gda_FTFrhs.m gdama_04.mlx gda_PD.m gdama_05.mlx gda_SD.m gdama_06.mlx gda_apply_response.m gdama_07.mlx gda_arrow3.m gdama_08.mlx gda_bin2int.m gdama_09.mlx gda_chebyshevfilt.m gdama_10.mlx gda_delay.m gdama_11.mlx gda_dlsmul.m gdama_12.mlx gda_dlsrhs.m gdama_13.mlx gda_draw.m gdama_14.mlx gda_filtermul.m gdama_15.mlx gda_filterrhs.m gdama_16.mlx gda_gray1table.mat gdama_17.mlx gda_gray2int.m gdastereo.m gda_int2bin.m gray1table.mat gda_int2gray.m
% gdama01_03
% example of simple algebra, c=a+b with a=4.5 and b=5.1
 
a=4.5;
b=5.1;
c=a+b;
disp('c');
c
disp(c);
9.6000
% gdama01_04
% example of powers & square root
% c=sqrt(a^2+b^2); with a=3 and b=4
 
a=6;
b=8;
c=sqrt(a^2+b^2);
disp('c');
c
disp(c);
10
% gdama01_05
% example of trig function
% c=sin(n*pi*(x-x0)/L) with n=3, x=4, x0=1; L=6
% supports Section I.3
 
n=3; x=4; x0=1; L=6;
c=sin(n*pi*(x-x0)/L);
disp('c');
c
disp(c);
-1
% gdama01_06
% example of defining vectors and matrices
 
% define row-vector, column-vector and a matrix
r = [2, 4, 6];
c = [1, 3, 5]';
disp('r');
r
disp(r);
2 4 6
disp(' ');
disp('c');
c
disp(c);
1 3 5
disp(' ');
 
% initialze vectors to zeroes
N = 3;
r = zeros(1,N);
c = zeros(N,1);
disp('r');
r
disp(r);
0 0 0
disp(' ');
disp('c');
c
disp(c);
0 0 0
disp(' ');
 
% initialze vectors to ones
N = 3;
r = ones(1,N);
c = ones(N,1);
disp('r');
r
disp(r);
1 1 1
disp(' ');
disp('c');
c
disp(c);
1 1 1
disp(' ');
 
% get length of c
N = length(c);
disp('N=length(c)');
N=length(c)
disp(N);
3
disp(' ');
 
% initialze vectors to ones
a = r(2);
b = c(3);
disp('a=r(2)');
a=r(2)
disp(a);
1
disp(' ');
disp('b=c(3)');
b=c(3)
disp(b);
1
disp(' ');
 
% define A
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
disp('A');
A
disp(A);
1 2 3 4 5 6 7 8 9
disp(' ');
 
% define A, alternate version
A = [1, 2, 3; 4, 5, 6; 7, 8, 9];
disp('A');
A
disp(A);
1 2 3 4 5 6 7 8 9
disp(' ');
 
% initialize matrices with zeros() and ones()
N=3; M=3;
A = zeros(N,M);
B = ones(N,M);
disp('A');
A
disp(A);
0 0 0 0 0 0 0 0 0
disp(' ');
disp('B');
B
disp(B);
1 1 1 1 1 1 1 1 1
disp(' ');
 
[N,M]=size(A);
disp('N = rows in A');
N = rows in A
disp(N);
3
disp(' ');
disp('M = cols in A');
M = cols in A
disp(M);
3
disp(' ');
 
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
disp('A');
A
disp(A);
1 2 3 4 5 6 7 8 9
disp(' ');
 
B = [ [0, 3, 7]', [5, 5, 9]', [1, 2, 7]'];
disp('B');
B
disp(B);
0 5 1 3 5 2 7 9 7
disp(' ');
 
s = A(1,3);
disp('s = M(1,3)');
s = M(1,3)
disp(s);
3
disp(' ');
 
S = A+B;
disp('S=A+B');
S=A+B
disp(S);
1 7 4 7 10 8 14 17 16
disp(' ');
 
D = A-B;
disp('D=A-B');
D=A-B
disp(D);
1 -3 2 1 0 4 0 -1 2
disp(' ');
 
P = A*B;
disp('P=AB');
P=AB
disp(P);
27 42 26 57 99 56 87 156 86
disp(' ');
 
% gdama01_07
% examples of matrix products
 
% define column-vector a
a = [1, 2, 3]';
disp('a');
a
disp(a);
1 2 3
disp(' ');
 
% define column-vector b
b = [2, 4, 6]';
disp('b');
b
disp(b);
2 4 6
disp(' ');
 
% define matrix A
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
disp('A');
A
disp(A);
1 2 3 4 5 6 7 8 9
disp(' ');
 
% define matrix N
B = [ [0, 3, 7]', [5, 5, 9]', [1, 2, 7]'];
disp('B');
B
disp(B);
0 5 1 3 5 2 7 9 7
disp(' ');
 
% inner product
s = a'*b;
disp('s=aT*b');
s=aT*b
disp(s);
28
disp(' ');
 
% outer product
T = a*b';
disp('T = a*bT');
T = a*bT
disp(T);
2 4 6 4 8 12 6 12 18
disp(' ');
 
% matrix times vector
c = A*a;
disp('c = A b');
c = A b
disp(c);
14 32 50
disp(' ');
 
% matrix times matrix
P = A*B;
disp('P = A B');
P = A B
disp(P);
27 42 26 57 99 56 87 156 86
disp(' ');
 
% identity matrix
I = eye(3);
disp('I');
I
disp(I);
1 0 0 0 1 0 0 0 1
disp(' ');
 
% elementwise product
d = a.*b;
disp('d = a .* b');
d = a .* b
disp(d);
2 8 18
disp(' ');
% gdama01_08
% sub-matrices using colon operator
 
a = [1, 2, 3]';
disp('a');
a
disp(a);
1 2 3
disp(' ');
 
% define matrix B
B = [ 1, 2, 3; 4, 5, 6; 7, 8, 9 ];
disp('B');
B
disp(B);
1 2 3 4 5 6 7 8 9
disp(' ');
 
% top two elements of vector a
x = a(1:2);
disp('x is top two elements of a');
x is top two elements of a
disp(x);
1 2
disp(' ');
 
% second columns of matrix B
y = B(:,2);
disp('y is second column of B');
y is second column of B
disp(y);
2 5 8
disp(' ');
 
% second row of matrix B
z = B(2,:);
disp('z is second row of B');
z is second row of B
disp(z);
4 5 6
disp(' ');
 
% lower-right 2x2 block of B
T = B(2:3,2:3);
disp('T is lower-right 2x2 block of B');
T is lower-right 2x2 block of B
disp(T);
5 6 8 9
disp(' ');
 
% counting by twos
r=[1:2:9];
disp('r=[1:2:9]');
r=[1:2:9]
disp(r);
1 3 5 7 9
disp(' ');
 
% counting down
r=[10:-1:1];
disp('r=[10:-1:1]');
r=[10:-1:1]
disp(r);
10 9 8 7 6 5 4 3 2 1
disp(' ');
 
% gdama01_09
% t-axis
 
N=11;
Dt = 2;
t = Dt*[0:N-1]';
 
disp('t-axis, transposed');
t-axis, transposed
disp(t');
Columns 1 through 10 0 2 4 6 8 10 12 14 16 18 Column 11 20
disp(' ');
% gdama01_10
% example of the matrix inverse, determinant, eigenvalues
 
% define matrix A
A = [ [1, 5, 13]', [2, 7, 17]', [3, 11, 19]'];
disp('A');
A
disp(A);
1 2 3 5 7 11 13 17 19
disp(' ');
 
% define vector b
b = [1, 2, 3]';
disp('b');
b
disp(b);
1 2 3
disp(' ');
 
% compute inverse
B = inv(A);
disp('B=inv(A)');
B=inv(A)
disp(B);
-2.2500 0.5417 0.0417 2.0000 -0.8333 0.1667 -0.2500 0.3750 -0.1250
disp(' ');
 
% check result
C=A*B;
disp('B=A*inv(A)');
B=A*inv(A)
disp(C);
1.0000 0 0.0000 0.0000 1.0000 0.0000 -0.0000 0 1.0000
disp(' ');
 
% check result
D=B*A;
disp('D=inv(A)*A');
D=inv(A)*A
disp(D);
1.0000 0.0000 -0.0000 0 1.0000 0.0000 0 0 1.0000
disp(' ');
 
% determinant
d = det(A);
disp('d=det(A)');
d=det(A)
disp(d);
24.0000
disp(' ');
 
% solve b = A*c
c = A\b;
disp('c=A\b');
c=A\b
disp(c);
-1.0417 0.8333 0.1250
disp(' ');
 
% check result
error=b-A*c;
disp('error');
error
disp(error);
1.0e-15 * 0.1110 -0.8882 0.4441
disp(' ');
 
% define matrix B
B = [ [1, 3, 4]', [2, 3, 2]', [0, 0, 4]'];
disp('B');
B
disp(B);
1 2 0 3 3 0 4 2 4
disp(' ');
 
% D=B*inv(A)
D=B/A;
disp('D=B/A');
D=B/A
disp(D);
1.7500 -1.1250 0.3750 -0.7500 -0.8750 0.6250 -6.0000 2.0000 0.0000
disp(' ');
 
% define symmetric matrix M
A = [ [1, 2, 0]', [2, 2, 0]', [0, 0, 4]'];
disp('A');
A
disp(A);
1 2 0 2 2 0 0 0 4
disp(' ');
 
% eigenvalues and eigenvectors
[V,LAMBDA] = eig(A);
disp('V');
V
disp(V);
-0.7882 0.6154 0 0.6154 0.7882 0 0 0 1.0000
disp(' ');
disp('LAMBDA');
LAMBDA
disp(LAMBDA);
-0.5616 0 0 0 3.5616 0 0 0 4.0000
disp(' ');
 
% check orthonormality
P=V'*V;
disp('P=V*V');
P=V*V
disp(P);
1 0 0 0 1 0 0 0 1
disp(' ');
 
% error E=M-V*LAMBDA*V';
E=M-V*LAMBDA*V'
E = 3×3
2.0000 1.0000 3.0000 1.0000 1.0000 3.0000 3.0000 3.0000 -1.0000
disp('E=M-V*LAMBDA*VT');
E=M-V*LAMBDA*VT
disp(E);
2.0000 1.0000 3.0000 1.0000 1.0000 3.0000 3.0000 3.0000 -1.0000
disp(' ');
% gdama01_11
% character string formatting amd lists (cellstrs) of character strings
 
myfilename='mydata.txt';
disp(myfilename);
mydata.txt
 
% format an integer
i=1;
myfilename=sprintf('myfile_%d.txt',i);
disp(myfilename);
myfile_1.txt
 
% format a fractioanl number
x=10.40;
mysentence = sprintf('position x is %.2f meters',x);
disp(mysentence);
position x is 10.40 meters
 
% alternative command
fprintf('position x is %.2f meters\n',x);
position x is 10.40 meters
 
% define list f character strings
mycolorlist = {'red', 'crimson', 'chartruse', 'teal'};
 
% access one element of the list
myfavoritecolor = char(mycolorlist(4));
disp('myfavoritecolor');
myfavoritecolor
disp(myfavoritecolor);
teal
% gdama01_12
% example of a for loop
 
% define a matrix A and a vector b
A = [ [1, 2, 3]', [4, 5, 6]', [7, 8, 9]'];
disp('A');
A
disp(A);
1 4 7 2 5 8 3 6 9
disp(' ');
 
% initialize b
b = [0, 0, 0]';
 
% use a for loop to extract the diagonal elements of
% the matrix M into the column-vector a
for i = (1:3)
b(i) = A(i,i);
end
 
disp('b');
b
disp(b);
1 5 9
disp(' ');
% gdama01_13
% example of two nested for loops
 
% define a matrices M and N
A = [ [1, 4, 7]', [2, 5, 8]', [3, 6, 9]'];
B = zeros(3,3);
 
% use two nested for loop to reverse the order
% of the elements in each row of the matrix M
for i = (1:3)
for j = (1:3)
B(i,4-j) = A(i,j);
end
end
 
disp('A');
A
disp(A);
1 2 3 4 5 6 7 8 9
disp(' ');
 
disp('B');
B
disp(B);
3 2 1 6 5 4 9 8 7
disp(' ');
% gdama01_14
% example of a for loop containing a consitional statement
 
% define vectors a and b
a = [ 1, 2, 1, 4, 3, 2, 6, 4, 9, 2, 1, 4 ]';
N = length(a);
b = zeros(N,1);
 
% use two nested for loop to reverse the order
% of the elements in each row of the matrix M
b = a;
for i = (1:N)
if ( a(i) >= 6 )
b(i) = 6;
else
b(i) = a(i);
end
end
 
% display vectors side by side
disp('a b');
a b
disp([a, b]);
1 1 2 2 1 1 4 4 3 3 2 2 6 6 4 4 9 6 2 2 1 1 4 4
% gdama01_15
% avoiding loops
 
% define a matrix A
A = [ [1, 2, 3]', [4, 5, 6]', [7, 8, 9]'];
disp('A');
A
disp(A);
1 4 7 2 5 8 3 6 9
disp(' ');
 
% extract diagonal
a = diag(A);
disp('a=diag(A)');
a=diag(A)
disp(a);
1 5 9
disp(' ');
 
% reverse order of rows
B = fliplr(A);
disp('B = fliplr(A)');
B = fliplr(A)
disp(B);
7 4 1 8 5 2 9 6 3
disp(' ');
 
% another way of reversing order order of rows
C = A(:,end:-1:1);
disp('A');
A
disp(A);
1 4 7 2 5 8 3 6 9
disp(' ');
 
% define vector a
a = [ 1, 2, 1, 4, 3, 2, 6, 4, 9, 2, 1, 4 ]';
disp('a');
a
disp(a);
1 2 1 4 3 2 6 4 9 2 1 4
disp(' ');
 
% clip to 6
b=a;
b(a>6)=6;
% display vectors side by side
disp('a b (via logical indexing)');
a b (via logical indexing)
disp([a, b]);
1 1 2 2 1 1 4 4 3 3 2 2 6 6 4 4 9 6 2 2 1 1 4 4
disp(' ');
 
% yet another way to clip to 6
b=a;
b(find(a>6))=6;
% display vectors side by side
disp('a b (via find)');
a b (via find)
disp([a, b]);
1 1 2 2 1 1 4 4 3 3 2 2 6 6 4 4 9 6 2 2 1 1 4 4
disp(' ');
 
% yet another way to clip to 6
b = a.*(a<6)+6.*(a>=6);
% display vectors side by side
disp('a b (via logical functions)');
a b (via logical functions)
disp([a, b]);
1 1 2 2 1 1 4 4 3 3 2 2 6 6 4 4 9 6 2 2 1 1 4 4
disp(' ');
 
 
% gdama01_16
% load and plot 1965-2017 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.
 
% load
D=load('../Data/global_temp.txt');
t=D(:,1);
d=D(:,2);
N=length(d);
 
% display first few lines
disp('Time Temp');
Time Temp
disp([t(1:5),d(1:5)]);
1.0e+03 * 1.9650 -0.0001 1.9660 -0.0001 1.9670 -0.0000 1.9680 -0.0001 1.9690 0.0001
disp(' ');
 
% plot data
figure();
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',12);
hold on;
axis( [1965, 2020, -0.5, 1.0] );
plot(t,d,'r-','LineWidth',3); % plot data as red lines
plot(t,d,'ko','LineWidth',3); % plot data as black circles
xlabel('calendar year');
ylabel('temperature anomaly, deg C');
title('global temperature data');
 
% add a caption
fprintf('Caption: Global temperature data for the time period 1965–2016\n');
Caption: Global temperature data for the time period 1965–2016
 
% gdama01_17
% write data to a text file
% this example writes global temperature data in deg F
 
cf = 1.8; % decC to degF conversion factor
 
D=load('../Data/global_temp.txt');
t=D(:,1);
d=D(:,2);
N=length(d);
 
% conversion
dF = cf*d;
 
% new file name
filename = '../TestFolder/global_temp_in_F.txt';
 
% concatenate output array
DF = [t,dF];
 
% write file with tab as delimiter
dlmwrite(filename,DF,'delimiter','\t');
 
% plot to check
figure(1);
clf;
set(gca,'LineWidth',3);
set(gca,'FontSize',14);
hold on;
axis( [1965, 2020, -1.0, 2.0] );
plot(t,dF,'r-','LineWidth',3);
plot(t,dF,'ko','LineWidth',3);
xlabel('calendar year');
ylabel('temperature anomaly, deg F');
% format a title string
titlestr=sprintf('temperature data from %d to %d',t(1),t(N));
% set the title to the title string
title(titlestr);
 
fprintf('Caption: Global temperature data for the time period 1965–2016\n');
Caption: Global temperature data for the time period 1965–2016
 
% gdama01_18
% Better quality plots
 
% Reference for 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.
 
% load
D=load('../Data/global_temp.txt');
t=D(:,1);
d=D(:,2);
N=length(d);
 
disp(sprintf('Number of data: %d\n', N));
Number of data: 52
 
% display first few lines
disp('Time Temp');
Time Temp
disp([t(1:5),d(1:5)]);
1.0e+03 * 1.9650 -0.0001 1.9660 -0.0001 1.9670 -0.0000 1.9680 -0.0001 1.9690 0.0001
disp(' ');
 
% plot data
figure(1);
clf; % clears the figure (just to be sure were starting with a blank figure)
% The following 'set' command specifies the line width of the figure axes.
% The figure does not look good when printed to paper when the line width is too thin
set(gca,'LineWidth',3);
% The following 'set' command specifies the font size of characters om the figure axes.
set(gca,'FontSize',14);
hold on; % dont let the above properties be changed by the plot() command
axis( [1965, 2020, -0.5, 1.0] );
plot(t,d,'r-','LineWidth',3);
plot(t,d,'ko','LineWidth',3);
xlabel('calendar year');
ylabel('temperature anomaly, deg C');
% format a title string
titlestr=sprintf('temperature data from %d to %d',t(1),t(N))
titlestr = 'temperature data from 1965 to 2016'
% set the title to the title string
title(titlestr);
 
fprintf('Caption: Global temperature data for the time period 1965–2016\n');
Caption: Global temperature data for the time period 1965–2016