-
Notifications
You must be signed in to change notification settings - Fork 5
/
ddTos.m
73 lines (66 loc) · 2.29 KB
/
ddTos.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
function varargout=ddTos(k,th,params,phi,xi,pxm,xver)
% [trWfax,dTinvdDdf2]=ddTos(k,th,params,phi,xi,pxm,xver)
%
% Calculates functions of the second derivatives of isotropic Tinv in
% Olhede & Simons (2013).
%
% k Wavenumber(s) at which this is to be evaluated [1/m]
% th The parameter vector with TWO elements (rest ignored)
% D Isotropic flexural rigidity [Nm]
% f2 The sub-surface to surface initial loading ratio
% params A structure with AT LEAST these constants that are known:
% DEL surface and subsurface density contrast [kg/m^3]
% g gravitational acceleration [m/s^2]
% phi Optionally precalculated phi, see PHIOS
% xi Optionally precalculated xi, see XIOS
% pxm Optionally precalculated (phi*xi-1)
%
% OUTPUT:
%
% trWfax The sum of the eigenvalues of dTinvdDdf2
% dTinvdDdf2 The second (cross) derivative of Tinv to D and f2
%
% Last modified by fjsimons-at-alum.mit.edu, 12/20/2010
defval('xver',0)
% Extract the parameters from the input
D=th(1);
f2=th(2);
DEL=params.DEL;
g=params.g;
if nargin<6
% First the auxiliary quantities
phi=phios(k,D,DEL,g);
xi = xios(k,D,DEL,g);
% Note that this has a zero at zero wavenumber
pxm=(phi.*xi-1);
end
% Produce the sum of the eigenvalues
warning off MATLAB:divideByZero
fax=2*k(:).^4/g/f2*dpos(DEL,-1,-1).*xi.^(-1)./pxm;
warning on MATLAB:divideByZero
trWfax=fax.*(dpos(DEL,0,1)*xi.^2-dpos(DEL,1,0)*pxm);
if nargout>=2 || xver==1
% First dTinvdDdf2, which is a symmetric matrix
warning off MATLAB:divideByZero
fax=2*k(:).^4/g/f2^2*dpos(DEL,-2,1).*xi./pxm.^3;
warning on MATLAB:divideByZero
dTinvdDdf2(:,1)=fax.*(1+dpos(DEL,-1,1)*xi.^2);
dTinvdDdf2(:,2)=fax.*(dpos(DEL,-1,1)*...
[dpos(DEL,-1,1)*xi.^3-1/2*phi.*xi.^2+3/2*xi]);
dTinvdDdf2(:,3)=fax.*(dpos(DEL,-2,2)*...
[dpos(DEL,-1,1)*xi.^4-phi.*xi.^3+2*xi.^2]);
else
dTinvdDdf2=NaN;
end
if xver==1
[invT,detT,L]=Tos(k,th,params,phi,xi,pxm);
randi=ceil(rand*length(dTinvdDdf2));
M=dTinvdDdf2;
Mrand=[M(randi,1) M(randi,2) ; M(randi,2) M(randi,3)];
Lrand=[L(randi,1) 0 ; L(randi,2) L(randi,3)];
% Check RELATIVE precision
difer((trWfax(randi)-trace(Lrand'*Mrand*Lrand))/trWfax(randi))
end
% And now for the output
varns={trWfax,dTinvdDdf2};
varargout=varns(1:nargout);