forked from ABHModels/raytransfer
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmetric.cpp
More file actions
313 lines (297 loc) · 16.3 KB
/
metric.cpp
File metadata and controls
313 lines (297 loc) · 16.3 KB
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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
/* Metric and derivative of metric used throughout code */
/* If possible, use Maple or Mathematica to optimize code */
void metric(long double r, long double th, long double g[][4])
{
long double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19;
long double gtt, grr, gthth, gpp, gtp;
t1 = r * r;
t2 = pow(t1, 0.2e1);
t3 = r * t1;
t4 = t1 * t2;
t5 = a22 + t1;
t6 = sin(th);
t6 = pow(t6, 0.2e1);
t7 = spin * spin;
t8 = cos(th);
t8 = t7 * pow(t8, 0.2e1);
t9 = (t8 + t1) * r + epsi3;
t10 = t7 + t1;
t11 = a13 + t3;
t12 = -t7 * r * t5 * t6 + t10 * t11;
t13 = 0.1e1 / r;
t8 = epsi3 * t13 + t1 + t8;
t14 = 0.2e1 * r;
t15 = -t14 + t10;
t16 = pow(t13, 0.2e1);
t17 = a52 * t16 + 0.1e1;
t18 = a22 * t16 + 0.1e1;
t16 = t10 * (a13 * t13 * t16 + 0.1e1);
t19 = -t7 * t18 * t6 + t16;
t15 = 0.1e1 / t15;
t17 = 0.1e1 / t17;
t19 = pow(t19, -0.2e1);
t12 = pow(t12, -0.2e1);
gtt = (0.2e1 * r * t2 - t7 * (-t6 * pow(t5, 0.2e1) + t2) - t4) * r * t9 * t12;
grr = t8 * t15 * t17;
gthth = t8;
gpp = t6 * (t6 * (-t4 * t7 * t10 + 0.2e1 * t3 * t2 * t7) + pow(t10, 0.2e1) * pow(t11, 0.2e1)) * t9 * t12 * t13;
gtp = -spin * (t16 * t18 - t1 + t14 - t7) * t8 * t6 * t19;
g[0][0] = gtt;
g[0][3] = gtp;
g[1][1] = grr;
g[2][2] = gthth;
g[3][0] = g[0][3];
g[3][3] = gpp;
}
void metric_rderivatives(long double r, long double th, long double dg[][4])
{
long double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20, t21, t22, t23, t24, t25, t26, t27, t28, t29;
long double dgttdr, dgtpdr, dgppdr;
t1 = r * r;
t2 = pow(t1, 0.2e1);
t3 = r * t1;
t4 = t1 * t2;
t5 = a22 + t1;
t6 = sin(th);
t6 = pow(t6, 0.2e1);
t7 = spin * spin;
t8 = t5 * t6;
t9 = cos(th);
t9 = t7 * pow(t9, 0.2e1);
t10 = (t1 + t9) * r + epsi3;
t11 = t1 + t7;
t12 = a13 + t3;
t13 = t11 * t12;
t14 = t8 * t7 * r - t13;
t15 = 2;
t5 = (double) t15 * r * t2 - t7 * (-t6 * pow(t5, 0.2e1) + t2) - t4;
t16 = 0.3e1 * t1;
t17 = t16 + t9;
t18 = (0.5e1 / 0.2e1 * t1 + 0.3e1 / 0.2e1 * t7) * r + a13;
t16 = r * t18 - t7 * (a22 + t16) * t6 / 0.2e1;
t14 = 0.1e1 / t14;
t19 = pow(t14, 0.2e1);
t4 = t6 * ((double) t15 * t3 * t2 * t7 - t4 * t7 * t11) + pow(t11, 0.2e1) * pow(t12, 0.2e1);
t12 = 0.1e1 / r;
t20 = t10 * t12;
t21 = t7 * a22;
t22 = 0.3e1 / 0.5e1;
t23 = pow(t12, 0.2e1);
t24 = t12 * t23;
t25 = a13 * t24 + 0.1e1;
t26 = a22 * t23 + 0.1e1;
t27 = t11 * t25;
t28 = (double) t15 * r;
t29 = -t7 * t26 * t6 + t27;
t29 = 0.1e1 / t29;
dgttdr = 0.4e1 * t10 * r * t19 * (r * (t3 * (-0.3e1 / 0.2e1 * r + 0.5e1 / 0.2e1) - t7 * (t1 - t8)) + t16 * t5 * t14) + t5 * t19 * (r * t17 + t10);
dgppdr = 0.4e1 * t6 * t10 * t19 * (t13 * t18 + (0.7e1 / 0.2e1 * (-0.4e1 / 0.7e1 * r + 0.1e1) * r - 0.3e1 / 0.2e1 * t7) * t7 * t2 * t6 + t4 * t16 * t14 * t12) + t6 * t4 * t19 * t12 * (t17 - t20);
dgtpdr = t6 * spin * (0.5e1 * t20 * (a13 * (t1 * (t1 / 0.5e1 + (a22 + t7) * t22) + t21) - 0.2e1 / 0.5e1 * t3 * (t3 - t21)) * t19 + (t27 * t26 - t1 + t28 - t7) * pow(t29, 0.2e1) * (epsi3 * t23 - t28 + (double) t15 * (epsi3 * t12 + t1 + t9) * t29 * ((double) t15 * (t21 * t24 * t6 + r * t25) - 0.3e1 * t11 * a13 * pow(t23, 0.2e1))));
dg[0][0] = dgttdr;
dg[0][3] = dgtpdr;
dg[3][0] = dg[0][3];
dg[3][3] = dgppdr;
}
void metric_r2derivatives(long double r, long double th, long double dg2[][4])
{
long double t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20, t21, t22, t23, t24, t25, t26, t27, t28, t29, t30, t31, t32, t33, t34, t35, t36, t37, t38, t39, t40, t41, t42, t43, t44, t45, t46, t47, t48, t49, t50, t51, t52, t53, t54, t55, t56, t57, t58, t59, t60, t61, t62, t63, t64, t65, t66, t67, t68, t69, t70, t71, t72, t73, t74, t75, t76, t77, t78, t79, t80, t81, t82, t83, t84, t85, t86, t87, t88, t89, t90, t91, t92, t93, t94, t95, t96, t97, t98, t99, t100, t101, t102, t103, t104, t105, t106, t107, t108, t109, t110, t111;
long double dgttdr2, dgtpdr2, dgppdr2;
t1 = r * r;
t2 = pow(t1, 0.2e1);
t3 = pow(t2, 0.2e1);
t4 = r * t2;
t5 = r * t1;
t6 = t5 * t3;
t7 = r * t3;
t8 = t5 * t2;
t9 = t1 * t2;
t10 = a22 + t1;
t11 = epsi3 + t5;
t12 = sin(th);
t13 = 2;
t14 = (double) t13 * a22;
t15 = spin * spin;
t16 = pow(t15, 0.2e1);
t17 = t15 * t16;
t18 = 0.3e1 * t15;
t19 = a22 * a22;
t20 = a22 * t19;
t21 = 0.5e1 * t15;
t22 = 0.6e1 * a22;
t23 = a22 - t18;
t24 = t5 * a22;
t25 = a13 + t5;
t26 = 0.3e1 / 0.2e1;
t27 = 0.1e1 / 0.2e1;
t28 = t19 * t26;
t29 = a13 * a22;
t30 = 0.4e1 * a13;
t31 = cos(th);
t32 = t26 * t15;
t33 = t22 - t32;
t34 = a13 + epsi3;
t35 = 0.4e1 * t15;
t36 = 0.12e2 * a22;
t37 = a13 / 0.7e1;
t38 = 0.6e1 / 0.7e1;
t39 = 0.3e1 * a22;
t40 = 0.6e1 * epsi3;
t41 = a13 * epsi3;
t42 = (double) t13 * a13;
t43 = epsi3 + t42;
t44 = a22 * epsi3;
t45 = -0.21e2 / 0.2e1 * epsi3;
t46 = t29 * epsi3;
t47 = t15 * t20;
t31 = pow(t31, 0.2e1);
t48 = t15 * t31;
t49 = (double) t13 * t15;
t50 = 36;
t51 = -0.26e2 / 0.3e1;
t52 = 0.2e1 / 0.3e1;
t53 = t52 * a22;
t54 = 0.12e2 * a13;
t55 = 0.13e2 / 0.3e1;
t56 = 0.5e1 / 0.2e1 * t19;
t57 = t29 * t15;
t58 = 0.4e1 * a22;
t59 = a13 * a13;
t60 = a13 * t59;
t61 = t59 * t16;
t62 = t61 * t19;
t63 = 0.3e1 * t59;
t64 = 0.16e2 * a13;
t65 = t59 * a22;
t66 = t26 * epsi3;
t67 = 0.10e2 / 0.3e1;
t68 = 0.7e1 / 0.2e1;
t69 = t27 * t15;
t70 = t67 * a13;
t71 = t52 * epsi3;
t72 = 0.3e1 * epsi3;
t73 = t27 * t19;
t74 = 0.28e2 / 0.3e1 * a22;
t75 = (double) t13 * epsi3;
t76 = t27 * epsi3;
t77 = -a13 * t52 - t76;
t78 = 0.15e2 / 0.2e1;
t79 = 0.8e1 / 0.3e1;
t80 = t79 * t15;
t81 = epsi3 / 0.4e1;
t82 = a13 + t81;
t83 = t19 * epsi3;
t84 = a13 - 0.8e1 / 0.7e1 * epsi3;
t85 = 0.1e1 / 0.3e1;
t86 = t85 * t15;
t87 = 0.18e2 * a22;
t88 = a13 + t75;
t89 = a22 * t88;
t90 = a22 - t15;
t91 = t59 * epsi3;
t92 = t52 * t15;
t93 = t3 * ((r * t27 - t79) * r + t92) - t30 * (0.7e1 / 0.6e1 * t16 * t84 + t83 + t86 * a22 * (a13 - 0.13e2 / 0.2e1 * epsi3));
t32 = t48 * (t1 * (((t1 * (t5 * ((double) t13 * t4 + 0.6e1 * t15 * (a13 * t51 + t53) + 0.6e1 * t29) + t35 * (a13 * (a13 - 0.14e2 * a22) + t56 * t15)) + t63 * (t15 * (-t58 + t21) + t19)) * r - t64 * t19 * t16) * r - 0.8e1 * t65 * (a22 - t32) * t15) + 0.12e2 * (a13 * (t15 * (-a22 * t55 - t49) + t19) + t2 * (-a13 + 0.16e2 / 0.3e1 * a22 + t1)) * t8 + t4 * (t5 * (t1 * (t1 * (t36 - t49) + a22 * (a22 * (double) t50 - t35)) + a13 * (a13 - 0.20e2 * a22) + t15 * (0.34e2 * t19 + t54)) + t57 * (-0.28e2 * a22 - 0.18e2 * t15)) + t62) / 0.6e1;
t50 = (int) (0.2e1 / 0.7e1);
t94 = 0.3e1 / 0.7e1;
t95 = (double) t50 * t15;
t96 = 0.1e1 / 0.42e2;
t97 = -0.23e2 / 0.42e2;
t98 = t15 * a13;
t99 = 0.4e1 / 0.3e1 * t15;
t100 = 0.1e1 / 0.5e1;
t101 = epsi3 / 0.14e2 + a13;
t102 = 0.4e1 / 0.5e1;
t103 = a13 * (a13 - 0.8e1 / 0.5e1 * epsi3);
t104 = a13 - 0.26e2 / 0.17e2 * epsi3;
t105 = t85 * epsi3;
t106 = a13 - t105;
t107 = 0.4e1 / 0.15e2;
t108 = a13 * t1;
t31 = r * (0.7e1 / 0.5e1 * t31 * (t1 * (t2 * ((t5 * ((-r / 0.14e2 + (double) t50) * r - t15 / 0.21e2) + 0.4e1 / 0.3e1 * t98) * r - t37 * (a13 + 0.16e2 * t15)) - 0.16e2 / 0.21e2 * t61) + (t1 * (((t5 * ((a13 * t94 - t95) * r - a13 * t38 + t16 * t96) + a13 * (0.25e2 / 0.21e2 * t16 + t37)) * r + t98 * (a13 * t97 - t95)) * r + t98 * (t16 * (double) t50 + 0.4e1 / 0.7e1 * a13)) + t61) * r - 0.5e1 / 0.14e2 * t17 * t59) + t108 * ((t16 * (-a13 * t68 + t75) + t41) * t107 - 0.9e1 / 0.5e1 * t2 * (-0.4e1 / 0.9e1 * t15 + a13 - 0.10e2 / 0.9e1 * epsi3))) + t15 * (-0.19e2 / 0.30e2 * t91 * t1 - t91 * t86);
t50 = 60;
t12 = pow(t12, 0.2e1);
t97 = pow(t12, 0.2e1);
t107 = pow(t10, 0.2e1);
t109 = (double) t13 * r;
t20 = (-0.4e1 * (double) t13 * t1 * ((t3 * (t34 - t22) + t47 * t34) * r + t41 * t19 * t23) + 0.12e2 * t4 * t3 - 0.28e2 * t2 * (t4 * (a22 * (epsi3 * t38 + a13 - t39) + t15 * (0.3e1 / 0.14e2 * epsi3 + t37)) + t46 * (a22 - 0.13e2 / 0.7e1 * t15)) + 0.4e1 * (t2 * (t19 * (a22 * t34 + t15 * (-0.9e1 * a13 + t45)) + t5 * (t1 * (t1 * (t33 + t1) + a22 * (t36 - t35)) + ((0.13e2 * a22 - 0.9e1 / 0.2e1 * t15) * a22 - t40) * a22 + t41)) + t48 * (t1 * (t1 * (t5 * (t2 * t27 + t28) - 0.3e1 * t29 * t23) + t30 * t19 * t15) + (double) t13 * t9 * (-a22 * (a13 - 0.5e1 * a22) + t18 * a13 + t24) + t4 * (t5 * (a13 - t14) + t19 * (t22 - t21)) + t25 * t15 * t20)) * r + 0.4e1 * t47 * t41 - 0.16e2 * t9 * (t15 * (-t42 * epsi3 - t20) + t44 * (-a22 * t26 + a13) + (a22 * (a13 + 0.9e1 / 0.4e1 * epsi3) + t15 * t43) * a22 * r)) * t16 * t97;
t23 = (0.12e2 * (t4 * t93 - t85 * t2 * (t5 * (t15 * (0.11e2 * a13 * (a13 - 0.54e2 / 0.11e2 * a22 - 0.30e2 / 0.11e2 * epsi3) + 0.6e1 * t44 + t28 * t15) + t29 * (a13 + 0.7e1 * epsi3)) + t41 * (-t15 * (a13 + t87) + t29) - 0.9e1 * t19 * (a13 - 0.5e1 / 0.6e1 * epsi3) * t16) - 0.4e1 / 0.3e1 * t9 * (t2 * (a22 * (a13 + t72) + t15 * (-a13 * t26 + epsi3 - t39)) - t15 * (0.6e1 * (-t76 + t19) * a13 - 0.6e1 * t83 - t14 * t82 * t15) - t46) + t1 * (t9 * ((t1 * ((a22 * t51 + t49 - t66) * r + t15 * (t14 - t69) - t19 * t68 + t70 - t71) - t15 * (t53 * t15 - 0.4e1 * a13 + 0.4e1 * t73) - a13 * (a13 - t75 - t74) - 0.32e2 / 0.3e1 * t44) * r + (-t44 * t78 + t80 * (a13 - 0.3e1 / 0.4e1 * epsi3)) * a22 + (-t71 + t19) * a13 + t16 * t77) - t91 * (a22 - 0.5e1 / 0.3e1 * t15) * t90) + t57 * (t5 * (t15 * (epsi3 * t55 - 0.5e1 * a13) + t89) + t41 * t90) - a13 * t16 * t19 * (a13 - t75) * r) * r - 0.12e2 * t32) * t15 * t12;
t28 = t15 + t1;
t32 = t15 * r;
t10 = t32 * t10 * t12 - t25 * t28;
t34 = 0.3e1 / 0.10e2 * t15;
t36 = -0.7e1 / 0.3e1;
t37 = t15 * a22;
t38 = 0.20e2 / 0.3e1 * a22;
t46 = 0.3e1 / 0.4e1 * t15;
t47 = 0.21e2;
t51 = 0.14e2 / 0.3e1;
t53 = -0.1e1 / 0.8e1;
t55 = 0.13e2 / 0.6e1;
t93 = epsi3 / 0.16e2;
t105 = a13 + t105;
t110 = t35 * a22;
t72 = -a13 + t72;
t111 = t68 * t15;
t33 = t48 * t27 * (t2 * (t1 * (t5 * (t1 * (-t5 * t52 + 0.4e1 / 0.3e1 * a13 + 0.64e2 / 0.3e1 * a22) + (t95 * (a22 + a13) + t29) * t51) - t79 * (t15 * (-t46 * t19 - 0.5e1 / 0.2e1 * a13 * (a13 - 0.14e2 / 0.5e1 * a22)) + t65)) + t59 * (t15 * (t15 * t47 - t58) + t19) * t85) + 0.4e1 * t1 * (t61 * a22 + t6) + t4 * (t1 * ((t1 * (-0.4e1 / 0.9e1 * t1 * t33 + a22 * (-0.20e2 / 0.3e1 * t15 + t74)) + a13 * (a13 - t38) - 0.4e1 * t15 * (-a13 - t56 + t37)) * r + t29 * (t35 + t38)) + t57 * (t58 + t49)) + t62);
t33 = t5 * (t9 * (((t1 * (0.4e1 * a22 * t55 + 0.4e1 * epsi3 * t53 + t1 * (r * t26 + t79) + 0.4e1 * a13 - 0.4e1 * t69) + 0.32e2 / 0.3e1 * a22 * (a13 - 0.5e1 / 0.16e2 * epsi3) + 0.32e2 / 0.3e1 * t15 * (0.7e1 / 0.16e2 * a13 - 0.3e1 / 0.8e1 * a22 + t93)) * r - 0.28e2 / 0.3e1 * a22 * t84 - 0.28e2 / 0.3e1 * t15 * (t94 * (t19 * t36 + a13) - t37 / 0.7e1)) * r + 0.7e1 * (t44 * t27 + 0.40e2 / 0.21e2 * t15 * (a13 - 0.9e1 / 0.20e2 * epsi3)) * a22 + 0.7e1 * (0.2e1 / 0.21e2 * epsi3 + t19) * a13 - t16 * t77) + t99 * t91 * (a22 + t111)) + t33;
t38 = -0.77e2 / 0.6e1;
t53 = t65 * t17;
t55 = 0.15e2 / 0.4e1;
t58 = 0.3e1 * a13;
t74 = 0.17e2 / 0.4e1;
t77 = t15 * t59;
t35 = t3 * (a13 * (t17 * (double) t13 * t43 + t91 + t98 * (t37 * t47 - t75)) + (t4 * (t1 * (((-a22 * t27 + t15 * t55 + t1) * r + t35 + t58 - t76) * r + t15 * (t14 + t111) + 0.8e1 * a13 - 0.5e1 * epsi3) + t15 * (t15 * (a22 * t68 + t46) + 0.14e2 * t101) + (0.33e2 / 0.4e1 * a13 - 0.3e1 * epsi3) * a13) + t15 * (t15 * (-0.18e2 * a13 * t106 + 0.6e1 * t37 * t106) - t63 * (a13 - 0.8e1 * epsi3)) - t65 * t72) * r) + t44 * t17 * t60;
t6 = t48 * t26 * (t1 * (((-t27 * t6 - a13 * (t16 * (a13 * t47 - t18 * a22) + t59 * t90) * t85) * r + 0.7e1 * t61 * (a22 + t69)) * r + t92 * t60 * t90) - (double) t13 * (t6 * (a13 - t15 + t1) - t53) * r + 0.4e1 * a13 * t4 * (-t4 * (a22 + t86) + t98 * (t15 * t67 + a22)) + t9 * (((-a13 * (t15 * (t39 - t92) + a13) + (t1 * (t1 * (-a22 * t67 - t86) + t15 * (-0.19e2 / 0.3e1 * a22 + t15 / 0.6e1) + 0.6e1 * a13) + t15 * (-t110 + t64) + t63) * r) * r - a22 * (t59 + t17) - t98 * (a13 * t38 - t49)) * r + t98 * (t14 * t15 - t30)) + t16 * t60 * t90);
t6 = t7 * ((-0.5e1 * t16 * (-0.11e2 / 0.5e1 * a13 - 0.2e1 / 0.5e1 * epsi3) - 0.5e1 * t103 + 0.21e2 * t37 * (a13 - 0.13e2 / 0.21e2 * epsi3) + (0.9e1 / 0.2e1 * a22 * (a13 - 0.13e2 / 0.9e1 * epsi3) + 0.9e1 / 0.2e1 * t15 * (t79 * (a13 + t93) - t86)) * t1) * r + 0.6e1 * t15 * (t59 * t74 + a22 * t16 / 0.6e1) - 0.6e1 * a13 * (-a22 * (a13 - t66) + t16)) + t6;
t18 = pow(t28, 0.2e1);
t25 = pow(t25, 0.2e1);
t14 = (0.12e2 * ((t1 * (t68 * (a22 * (a22 + epsi3 / 0.21e2) - t96 * t15 * epsi3) + t14 * t1) - 0.5e1 / 0.6e1 * t44 * (-t102 * t15 + a22)) * r + 0.5e1 / 0.3e1 * t48 * (-t100 * t24 * (r + 0.1e1) - t9 / 0.20e2 + r * (r * a22 * (-0.3e1 / 0.4e1 * a22 + t34) + t19) - t34 * t19)) * r + 0.12e2 * t27 * t3 * (-r + 0.1e1) + 0.12e2 * (t1 * (t1 * (t1 * (-t15 / 0.12e2 - 0.11e2 / 0.6e1 * a22) + a22 * (a22 * t36 - t86)) + a22 * (-epsi3 - 0.5e1 / 0.4e1 * t37)) + t83) * r - 0.12e2 * t81 * t15 * t19) * t15 * t9 * t12;
t24 = t32 * t12;
t14 = t24 * (0.6e1 * r * t33 + 0.6e1 * t85 * t57 * t1 * (t5 * (t15 * t72 + t89) + t41 * (a22 + 0.11e2 * t15)) + 0.6e1 * t9 * ((a13 * (t16 * (0.5e1 / 0.3e1 * a13 + t71) + t19 * t88 - t80 * t82 * a22) + ((t2 * (t1 * (a22 * t51 + t80) + t15 * (t22 + t69) + t19 * t78 - t70 + t71) + t15 * (a13 * (a13 + t71 - t87) + t44 * (double) t13) - t29 * t105 + t56 * t16) * r + t15 * (0.20e2 / 0.3e1 * (0.3e1 / 0.5e1 * epsi3 + t19) * a13 + 0.20e2 / 0.3e1 * t73 * epsi3 + t110 * (a13 - t71)) + t41 * (-0.4e1 / 0.3e1 * a22 + a13)) * r) * r - a22 * (-a22 * t27 * t43 * t16 + t91) + t21 * t41 * (a13 - 0.6e1 / 0.5e1 * a22)) + 0.6e1 * t62 * epsi3 + t14);
t19 = (t48 + t1) * r + epsi3;
t21 = a22 + t15;
t22 = t5 * (-t37 + t5);
t26 = a13 * (t1 * (t1 * t100 + 0.3e1 / 0.5e1 * t21) + t37) - 0.2e1 / 0.5e1 * t22;
t33 = t37 * t12;
t34 = 0.1e1 / r;
t36 = pow(t34, 0.2e1);
t38 = t34 * t36;
t39 = a13 * t38 + 0.1e1;
t44 = a22 * t36 + 0.1e1;
t46 = t28 * t39;
t47 = t46 * t44 - t1 + t109 - t15;
t44 = -t15 * t44 * t12 + t46;
t39 = (double) t13 * (r * t39 + t33 * t38) - t58 * t28 * pow(t36, 0.2e1);
t46 = epsi3 * t34 + t1 + t48;
t10 = 0.1e1 / t10;
t51 = pow(t10, 0.2e1);
t55 = pow(t51, 0.2e1);
t10 = t10 * t51;
t44 = 0.1e1 / t44;
t56 = pow(t44, 0.2e1);
t57 = t12 * spin;
t10 = -(double) t13 * t57 * (-((0.5e1 * t37 + t2) * a13 + 0.3e1 * t108 * t21 - (double) t13 * t22) * t19 * ((0.5e1 * t5 + t42) * r + 0.3e1 * t15 * t1 * (0.1e1 - t12) - t33) * t10 * t34 + t47 * t56 * (epsi3 * t38 - t46 * t44 * (t38 * (t34 * (t54 * t28 * t34 - 0.6e1 * t33) - 0.10e2 * a13) + (double) t13) + 0.1e1)) + 0.5e1 * t51 * t34 * t26 * t12 * spin * (-t34 * (-(double) t13 * t5 + epsi3 + t19) + 0.3e1 * t1 + t48) + 0.6e1 * t57 * (-t47 * t46 * pow(t56, 0.2e1) * pow(t39, 0.2e1) + (-(double) t13 * t2 + a13 * (t1 * t52 + t15) + (a13 + t32) * a22) * t19 * t51) + t57 * (-0.30e2 * t26 * ((t1 * t85 + t15) * a13 - t52 * r * (t33 + t2)) * t19 * t36 * t10 + 0.4e1 * t47 * (-epsi3 * t36 + t109) * t44 * t56 * t39);
t1 = (double) t13 * t12 * (pow(t18, 0.2e1) * pow(t25, 0.2e1) * t11 + t24 * (-0.4e1 * (double) t13 * t77 * t2 * (t41 * (a22 + t49) + t32 * (a22 * (a13 + t40) + t15 * (epsi3 * t74 - t42))) + 0.4e1 * t27 * t5 * (t5 * (a13 * (-0.9e1 * t15 * (t105 * a22 * t16 + t91) + t41 * (0.15e2 * t16 + t29)) + t2 * (a13 * (t16 * (-0.85e2 / 0.2e1 * a13 - 0.10e2 * epsi3) + t41 - 0.45e2 * t37 * (a13 - 0.7e1 / 0.15e2 * epsi3)) + (t3 + t15 * (0.34e2 * a13 * t104 - 0.45e2 * t37 * (a13 - 0.17e2 / 0.45e2 * epsi3)) + t17 * (-t66 - t30) + t59 * (a13 + t45)) * r)) + t53 * (a13 - 0.9e1 * epsi3)) - 0.4e1 * t5 * t6 - 0.4e1 * t68 * t77 * t1 * (t4 * (a22 * (a13 + 0.9e1 / 0.7e1 * epsi3) + t15 * (-0.13e2 / 0.7e1 * a13 + 0.109e3 / 0.14e2 * epsi3)) + t41 * t15 * (a22 + t15 / 0.7e1)) - 0.4e1 * t35 + t14)) * t55 * t38;
dgttdr2 = (t109 * t17 * pow(t107, 0.2e1) * t11 * t12 * t97 + t20 - t23 + (double) t50 * (t15 * t31 + t9 * (t15 * t8 - t17 * t85 * t43 - t91 + 0.68e2 / 0.3e1 * t98 * t104) / 0.10e2 - 0.2e1 / 0.5e1 * t2 * (-t3 * (a13 - t99 - t76) + t98 * (-0.6e1 * t15 * t106 + t41)) + r * (t3 * (((t16 / 0.15e2 - 0.16e2 / 0.15e2 * a13 + t71) * r + t15 * (t15 * t100 + a13 - 0.3e1 / 0.10e2 * epsi3)) * r + t15 * (-t16 / 0.30e2 - 0.28e2 / 0.15e2 * t101) - t27 * a13 * (a13 - 0.6e1 / 0.5e1 * epsi3)) + t61 * epsi3) + t4 * (a13 * (t16 * (-0.67e2 / 0.2e1 * a13 + 0.29e2 * epsi3) + t41) + t7) / 0.15e2 + t52 * (t16 * (a13 * t102 - epsi3 * t100) + t103) * t3) * t5) * t55;
dgppdr2 = t1;
dgtpdr2 = t10;
dg2[0][0] = dgttdr2;
dg2[0][3] = dgtpdr2;
dg2[3][0] = dg2[0][3];
dg2[3][3] = dgppdr2;
}
void uppermetric(long double r, long double th, long double rth[])
{
long double gurr, guthth;
long double t1, t2, t3, t4;
t1 = cos(th);
t2 = 0.1e1 / r;
t3 = r * r;
t4 = spin * spin;
t1 = t4 * pow(t1, 0.2e1) + epsi3 * t2 + t3;
t1 = 0.1e1 / t1;
gurr = t1 * (-0.2e1 * r + t3 + t4) * (a52 * pow(t2, 0.2e1) + 0.1e1);
guthth = t1;
rth[0] = gurr;
rth[1] = guthth;
}