-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgradHO_3D.m
189 lines (184 loc) · 5.32 KB
/
gradHO_3D.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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
function [DTDX,DTDY,DTDZ] = gradHO_3D(X,Y,z,T)
% Andrew Wheeler 2016
D16 = zeros(7,10);
D16(1,1)=-4192853913297.768d0/1.0d12;
D16(1,2)=5339040076397.708d0/1.0d12;
D16(1,3)=-995327078300.349d0/1.0d12;
D16(1,4)=-150859084799.5902d0/1.0d12;
D16(2,1)=-403794985785.2408d0/1.0d12;
D16(2,3)=431827299959.2449d0/1.0d12;
D16(2,4)=-44851702678.4067d0/1.0d12;
D16(2,5)=16819388504.4025d0/1.0d12;
D16(3,1)=157062810031.9159d0/1.0d12;
D16(3,2)=-900989848376.5991d0/1.0d12;
D16(3,4)=959478115560.413d0/1.0d12;
D16(3,5)=-250644037526.018d0/1.0d12;
D16(3,6)=35092960310.2883d0/1.0d12;
D16(4,1)=16968084745.7541d0/1.0d12;
D16(4,2)=66702567274.9842d0/1.0d12;
D16(4,3)=-683894144860.1648d0/1.0d12;
D16(4,5)=753863390649.9535d0/1.0d12;
D16(4,6)=-178653360538.646d0/1.0d12;
D16(4,7)=25013462728.1191d0/1.0d12;
D16(5,2)=-26519952061.4978d0/1.0d12;
D16(5,3)=189413141579.3246d0/1.0d12;
D16(5,4)=-799266426974.1559d0/1.0d12;
D16(5,6)=799266426974.1559d0/1.0d12;
D16(5,7)=-189413141579.3246d0/1.0d12;
D16(5,8)=26519952061.4978d0/1.0d12;
%
[ni,nj,nk] = size(T);
DTDX(1:ni,1:nj,1:nk) = 0.0;
DTDY(1:ni,1:nj,1:nk) = 0.0;
DTDZ(1:ni,1:nj,1:nk) = 0.0;
for k=1:nk
% body
I = 5:ni-4;
J = 1:nj;
DTI(I,J) = D16(5,2)*T(I-3,J) ...
+ D16(5,3)*T(I-2,J) ...
+ D16(5,4)*T(I-1,J) ...
- D16(5,2)*T(I+3,J) ...
- D16(5,3)*T(I+2,J) ...
- D16(5,4)*T(I+1,J) ;
DXI(I,J) = D16(5,2)*X(I-3,J) ...
+ D16(5,3)*X(I-2,J) ...
+ D16(5,4)*X(I-1,J) ...
- D16(5,2)*X(I+3,J) ...
- D16(5,3)*X(I+2,J) ...
- D16(5,4)*X(I+1,J) ;
DYI(I,J) = D16(5,2)*Y(I-3,J) ...
+ D16(5,3)*Y(I-2,J) ...
+ D16(5,4)*Y(I-1,J) ...
- D16(5,2)*Y(I+3,J) ...
- D16(5,3)*Y(I+2,J) ...
- D16(5,4)*Y(I+1,J) ;
%
for II=1:4
DTI(II,J) = D16(II,1)*T(1,J) ...
+ D16(II,2)*T(2,J) ...
+ D16(II,3)*T(3,J) ...
+ D16(II,4)*T(4,J) ...
+ D16(II,5)*T(5,J) ...
+ D16(II,6)*T(6,J) ...
+ D16(II,7)*T(7,J);
DXI(II,J) = D16(II,1)*X(1,J) ...
+ D16(II,2)*X(2,J) ...
+ D16(II,3)*X(3,J) ...
+ D16(II,4)*X(4,J) ...
+ D16(II,5)*X(5,J) ...
+ D16(II,6)*X(6,J) ...
+ D16(II,7)*X(7,J);
DYI(II,J) = D16(II,1)*Y(1,J) ...
+ D16(II,2)*Y(2,J) ...
+ D16(II,3)*Y(3,J) ...
+ D16(II,4)*Y(4,J) ...
+ D16(II,5)*Y(5,J) ...
+ D16(II,6)*Y(6,J) ...
+ D16(II,7)*Y(7,J);
%DTDX(1:ni,1:nj,1:nk) = 0.0;
end
%
for II=ni-3:ni
DTI(II,J) = -D16(ni-II+1,1)*T(ni,J) ...
- D16(ni-II+1,2)*T(ni-1,J) ...
- D16(ni-II+1,3)*T(ni-2,J) ...
- D16(ni-II+1,4)*T(ni-3,J) ...
- D16(ni-II+1,5)*T(ni-4,J) ...
- D16(ni-II+1,6)*T(ni-5,J) ...
- D16(ni-II+1,7)*T(ni-6,J);
DXI(II,J) = -D16(ni-II+1,1)*X(ni,J) ...
- D16(ni-II+1,2)*X(ni-1,J) ...
- D16(ni-II+1,3)*X(ni-2,J) ...
- D16(ni-II+1,4)*X(ni-3,J) ...
- D16(ni-II+1,5)*X(ni-4,J) ...
- D16(ni-II+1,6)*X(ni-5,J) ...
- D16(ni-II+1,7)*X(ni-6,J);
DYI(II,J) = -D16(ni-II+1,1)*Y(ni,J) ...
- D16(ni-II+1,2)*Y(ni-1,J) ...
- D16(ni-II+1,3)*Y(ni-2,J) ...
- D16(ni-II+1,4)*Y(ni-3,J) ...
- D16(ni-II+1,5)*Y(ni-4,J) ...
- D16(ni-II+1,6)*Y(ni-5,J) ...
- D16(ni-II+1,7)*Y(ni-6,J);
end
%
%
I = 1:ni;
J = 5:nj-4;
DTJ(I,J) = D16(5,2)*T(I,J-3) ...
+ D16(5,3)*T(I,J-2) ...
+ D16(5,4)*T(I,J-1) ...
- D16(5,2)*T(I,J+3) ...
- D16(5,3)*T(I,J+2) ...
- D16(5,4)*T(I,J+1) ;
DXJ(I,J) = D16(5,2)*X(I,J-3) ...
+ D16(5,3)*X(I,J-2) ...
+ D16(5,4)*X(I,J-1) ...
- D16(5,2)*X(I,J+3) ...
- D16(5,3)*X(I,J+2) ...
- D16(5,4)*X(I,J+1) ;
DYJ(I,J) = D16(5,2)*Y(I,J-3) ...
+ D16(5,3)*Y(I,J-2) ...
+ D16(5,4)*Y(I,J-1) ...
- D16(5,2)*Y(I,J+3) ...
- D16(5,3)*Y(I,J+2) ...
- D16(5,4)*Y(I,J+1) ;
%
for JJ=1:4
DTJ(I,JJ) = D16(JJ,1)*T(I,1) ...
+ D16(JJ,2)*T(I,2) ...
+ D16(JJ,3)*T(I,3) ...
+ D16(JJ,4)*T(I,4) ...
+ D16(JJ,5)*T(I,5) ...
+ D16(JJ,6)*T(I,6) ...
+ D16(JJ,7)*T(I,7);
DXJ(I,JJ) = D16(JJ,1)*X(I,1) ...
+ D16(JJ,2)*X(I,2) ...
+ D16(JJ,3)*X(I,3) ...
+ D16(JJ,4)*X(I,4) ...
+ D16(JJ,5)*X(I,5) ...
+ D16(JJ,6)*X(I,6) ...
+ D16(JJ,7)*X(I,7);
DYJ(I,JJ) = D16(JJ,1)*Y(I,1) ...
+ D16(JJ,2)*Y(I,2) ...
+ D16(JJ,3)*Y(I,3) ...
+ D16(JJ,4)*Y(I,4) ...
+ D16(JJ,5)*Y(I,5) ...
+ D16(JJ,6)*Y(I,6) ...
+ D16(JJ,7)*Y(I,7);
end
%
for JJ=nj-3:nj
DTJ(I,JJ) = -D16(nj-JJ+1,1)*T(I,nj) ...
- D16(nj-JJ+1,2)*T(I,nj-1) ...
- D16(nj-JJ+1,3)*T(I,nj-2) ...
- D16(nj-JJ+1,4)*T(I,nj-3) ...
- D16(nj-JJ+1,5)*T(I,nj-4) ...
- D16(nj-JJ+1,6)*T(I,nj-5) ...
- D16(nj-JJ+1,7)*T(I,nj-6);
DXJ(I,JJ) = -D16(nj-JJ+1,1)*X(I,nj) ...
- D16(nj-JJ+1,2)*X(I,nj-1) ...
- D16(nj-JJ+1,3)*X(I,nj-2) ...
- D16(nj-JJ+1,4)*X(I,nj-3) ...
- D16(nj-JJ+1,5)*X(I,nj-4) ...
- D16(nj-JJ+1,6)*X(I,nj-5) ...
- D16(nj-JJ+1,7)*X(I,nj-6);
DYJ(I,JJ) = -D16(nj-JJ+1,1)*Y(I,nj) ...
- D16(nj-JJ+1,2)*Y(I,nj-1) ...
- D16(nj-JJ+1,3)*Y(I,nj-2) ...
- D16(nj-JJ+1,4)*Y(I,nj-3) ...
- D16(nj-JJ+1,5)*Y(I,nj-4) ...
- D16(nj-JJ+1,6)*Y(I,nj-5) ...
- D16(nj-JJ+1,7)*Y(I,nj-6);
end
%
DTDX(:,:,k) = (DTI.*DYJ - DTJ.*DYI)./(DXI.*DYJ - DXJ.*DYI);
DTDY(:,:,k) = (DTI.*DXJ - DTJ.*DXI)./(DYI.*DXJ - DYJ.*DXI);
end
DTDZ(:,:,1) = (T(:,:,2) - T(:,:,1))/(z(2)-z(1));
DTDZ(:,:,end) = (T(:,:,end) - T(:,:,end-1))/(z(end)-z(end-1));
for k=2:length(z)-1
DTDZ(:,:,k) = (T(:,:,k+1)-T(:,:,k-1))/(z(k+1)-z(k-1));
end
return