1 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
2 /* fun-mpwx.c:
3
4 Automatically generated code implementing MPWX functional and
5 its derivatives. It is generated by func-codegen.pl being a part of
6 a "Automatic code generation framework for analytical functional
7 derivative evaluation", Pawel Salek, 2005
8
9 This functional is connected by making following changes:
10 1. add "extern Functional mpwxFunctional;" to 'functionals.h'
11 2. add "&mpwxFunctional," to 'functionals.c'
12 3. add "fun-mpwx.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.
13
14 mPWx exchange functional (actually gradient correction, for full
15 exchange functional, mPWx + Slater)
16
17 Implemented by David Wilson (david.wilson@latrobe.edu.au).
18 Sept 2007
19
20 NOTE: The original Adamo/Barone JCP reference had a typo [1] as published in [2].
21 (i) b is 0.00426, not 0.0046 as originally published
22 (ii) exponent (d) is 3.72, not 3.73.
23
24 References
25 [1] C. Adamo and V. Barone, J. Chem. Phys., 108, 664-675 (1998).
26 [2] Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 109, 5656-5667 (2005).
27
28
29 This functional has been generated from following input:
30 ------ cut here -------
31 rho: rhoa + rhob;
32 grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
33 zeta: (rhoa-rhob)/(rhoa+rhob);
34
35 Ax: -3.0/4.0*((6.0/%PI)^(1/3));
36 Axinv: 1.0/Ax;
37 bb: 0.00426;
38 beta: 5*((36*%PI)^(-5/3));
39 cc: 1.6455;
40 dd: 3.72;
41 em6: 10^(-6);
42 bmb: bb-beta;
43
44 xa: grada*rhoa^(-4/3);
45 xb: gradb*rhob^(-4/3);
46 xa2: xa^2;
47 xb2: xb^2;
48 xad: xa^(dd);
49 xbd: xb^(dd);
50
51 F0a: bb*xa^2-bmb*xa^2*exp(-cc*xa^2)-em6*xad;
52 F0b: bb*xb^2-bmb*xb^2*exp(-cc*xb^2)-em6*xbd;
53 F1a: 1+6*bb*xa*asinh(xa)-em6*xad*Axinv;
54 F1b: 1+6*bb*xb*asinh(xb)-em6*xbd*Axinv;
55 Exa: rhoa^(4.0/3.0)*F0a/F1a;
56 Exb: rhob^(4.0/3.0)*F0b/F1b;
57
58 K(rhoa,rhob,grada,gradb,gradab):=-(Exa+Exb);
59
60
61
62
63 ------ cut here -------
64 */
65
66
67 /* strictly conform to XOPEN ANSI C standard */
68 #if !defined(SYS_DEC)
69 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
70 * is not specified. */
71 #define _XOPEN_SOURCE 500
72 #define _XOPEN_SOURCE_EXTENDED 1
73 #endif
74 #include <math.h>
75 #include <stddef.h>
76 #include "general.h"
77
78 #define __CVERSION__
79
80 #include "functionals.h"
81
82 /* INTERFACE PART */
mpwx_isgga(void)83 static integer mpwx_isgga(void) { return 1; } /* FIXME: detect! */
84 static integer mpwx_read(const char *conf_line);
85 static real mpwx_energy(const FunDensProp* dp);
86 static void mpwx_first(FunFirstFuncDrv *ds, real factor,
87 const FunDensProp* dp);
88 static void mpwx_second(FunSecondFuncDrv *ds, real factor,
89 const FunDensProp* dp);
90 static void mpwx_third(FunThirdFuncDrv *ds, real factor,
91 const FunDensProp* dp);
92 static void mpwx_fourth(FunFourthFuncDrv *ds, real factor,
93 const FunDensProp* dp);
94
95 Functional mPWxFunctional = {
96 "mPWx", /* name */
97 mpwx_isgga, /* gga-corrected */
98 1,
99 mpwx_read,
100 NULL,
101 mpwx_energy,
102 mpwx_first,
103 mpwx_second,
104 mpwx_third,
105 mpwx_fourth
106 };
107
108 /* IMPLEMENTATION PART */
109 static integer
mpwx_read(const char * conf_line)110 mpwx_read(const char *conf_line)
111 {
112 fun_set_hf_weight(0);
113 return 1;
114 }
115
116 static real
mpwx_energy(const FunDensProp * dp)117 mpwx_energy(const FunDensProp *dp)
118 {
119 real res;
120 real rhoa = dp->rhoa, rhob = dp->rhob;
121 real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
122
123 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
124 real t11, t12;
125
126 t1 = pow(3.141592653589793,0.333333333333333);
127 t2 = 1/pow(rhoa,1.333333333333333);
128 t3 = grada*t2;
129 t4 = pow(t3,3.72);
130 t5 = pow(grada,2.0);
131 t6 = 1/pow(rhoa,2.666666666666667);
132 t7 = 0.00426-0.138888888888889/(pow(3.141592653589793,
133 1.666666666666667)*pow(36.0,0.666666666666667));
134 t8 = 1/pow(rhob,1.333333333333333);
135 t9 = gradb*t8;
136 t10 = pow(t9,3.72);
137 t11 = pow(gradb,2.0);
138 t12 = 1/pow(rhob,2.666666666666667);
139
140 /* code */
141 res = -1.0*(-1.0*t11*t12*t7/pow(2.718281828459045,1.6455*
142 t11*t12)+0.00426*t11*t12-1.0E-6*t10)*pow(rhob,1.333333333333333)/
143 (0.02556*gradb*asinh(t9)*t8+7.337616108654726E-7*t1*t10+1.0)-
144 1.0*(-1.0*t5*t6*t7/pow(2.718281828459045,1.6455*t5*t6)+0.00426*
145 t5*t6-1.0E-6*t4)*pow(rhoa,1.333333333333333)/(7.337616108654726E-7*
146 t1*t4+0.02556*grada*asinh(t3)*t2+1.0);
147
148 return res;
149 }
150
151 static void
mpwx_first_helper(real rhoa,real grada,real * res)152 mpwx_first_helper(real rhoa, real grada, real *res)
153 { real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
154 real t11, t12, t13, t14, t15, t16, t17, t18;
155 real t19;
156
157 t1 = pow(3.141592653589793,0.333333333333333);
158 t2 = 1/pow(rhoa,1.333333333333333);
159 t3 = grada*t2;
160 t4 = pow(t3,3.72);
161 t5 = asinh(t3);
162 t6 = 7.337616108654726E-7*t1*t4+0.02556*grada*t5*t2+1.0;
163 t7 = 1/
164 t6;
165 t8 = pow(grada,2.0);
166 t9 = 1/pow(rhoa,2.666666666666667);
167 t10 = 0.00426-0.138888888888889/(pow(3.141592653589793,
168 1.666666666666667)*pow(36.0,0.666666666666667));
169 t11 = 1/pow(2.718281828459045,1.6455*t8*t9);
170 t12 = -1.0*t10*t11*t8*t9+0.00426*t8*t9-1.0E-6*t4;
171 t13 = 1/sqrt(t8*t9+1.0);
172 t14 = 1/pow(rhoa,3.666666666666667);
173 t15 = 1/pow(rhoa,2.333333333333333);
174 t16 = 1/pow(t6,2.0);
175 t17 = pow(rhoa,1.333333333333333);
176 t18 = 1/pow(fabs(rhoa),3.626666666666667);
177 t19 = pow(grada,2.72);
178
179 /* code */
180 res[0] = -1.0*t17*t7*(-4.388*t10*t11*pow(grada,4.0)/pow(rhoa,
181 6.333333333333333)+4.96E-6*t15*t18*pow(grada,3.72)+2.666666666666667*
182 t10*t11*t14*t8-0.01136*t8*t14)-1.333333333333333*t12*t7*pow(rhoa,
183 0.333333333333333)+t12*t16*t17*(-3.639457589892744E-6*t1*t15*
184 pow(t3,2.72)*grada-0.03408*grada*t5*t15-0.03408*t8*t13*t14);
185 res[1] = t16*t17*t12*(2.729593192419558E-6*t1*t19*t2*
186 t18+0.02556*t5*t2+0.02556*grada*t13*t9)-1.0*t17*t7*(3.291*
187 t10*t11*pow(grada,3.0)/pow(rhoa,5.333333333333333)-2.0*t10*
188 t11*t9*grada+0.00852*grada*t9-3.72E-6*t19*t2*t18);
189 }
190
191 static void
mpwx_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)192 mpwx_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
193 {
194 real res[2];
195
196 mpwx_first_helper(dp->rhoa, dp->grada, res);
197 /* Final assignment */
198 ds->df1000 += factor*res[0];
199 ds->df0010 += factor*res[1];
200
201
202 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
203 fabs(dp->grada-dp->gradb)>1e-13)
204 mpwx_first_helper(dp->rhob, dp->gradb, res);
205 ds->df0100 += factor*res[0];
206 ds->df0001 += factor*res[1];
207
208 }
209
210 static void
mpwx_second_helper(real rhoa,real grada,real * res)211 mpwx_second_helper(real rhoa, real grada, real *res)
212 {
213 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
214 real t11, t12, t13, t14, t15, t16, t17, t18;
215 real t19, t20, t21, t22, t23, t24, t25, t26;
216 real t27, t28, t29, t30, t31, t32, t33, t34;
217 real t35, t36, t37, t38, t39, t40;
218
219 t1 = pow(3.141592653589793,0.333333333333333);
220 t2 = 1/pow(rhoa,1.333333333333333);
221 t3 = grada*t2;
222 t4 = pow(t3,3.72);
223 t5 = asinh(t3);
224 t6 = 7.337616108654726E-7*t1*t4+0.02556*grada*t5*t2+1.0;
225 t7 = 1/
226 t6;
227 t8 = pow(rhoa,0.333333333333333);
228 t9 = pow(grada,2.0);
229 t10 = 1/pow(rhoa,2.666666666666667);
230 t11 = 0.00426-0.138888888888889/(pow(3.141592653589793,
231 1.666666666666667)*pow(36.0,0.666666666666667));
232 t12 = 1/pow(2.718281828459045,1.6455*t9*t10);
233 t13 = -1.0*t10*t11*t12*t9-1.0E-6*t4+0.00426*t9*t10;
234 t14 = sqrt(t9*
235 t10+1.0);
236 t15 = 1/t14;
237 t16 = 1/pow(rhoa,3.666666666666667);
238 t17 = -0.03408*t9*t15*t16;
239 t18 = 1/pow(rhoa,2.333333333333333);
240 t19 = -0.03408*grada*t5*t18;
241 t20 = 1/pow(t6,2.0);
242 t21 = pow(rhoa,1.333333333333333);
243 t22 = pow(grada,4.0);
244 t23 = 1/pow(rhoa,6.333333333333333);
245 t24 = pow(grada,3.72);
246 t25 = fabs(rhoa);
247 t26 = 1/pow(t25,3.626666666666667);
248 t27 = 2.666666666666667*t11*t12*t16*t9+4.96E-6*t24*t18*
249 t26-0.01136*t9*t16-4.388*t11*t22*t23*t12;
250 t28 = pow(grada,3.0);
251 t29 = 1/pow(rhoa,5.333333333333333);
252 t30 = pow(grada,2.72);
253 t31 = -2.0*t10*t11*t12*grada-3.72E-6*t30*t2*t26+3.291*
254 t11*t28*t29*t12+0.00852*grada*t10;
255 t32 = 2.729593192419558E-6*t1*t30*t2*t26+0.02556*t5*t2+
256 0.02556*grada*t15*t10;
257 t33 = 1/pow(rhoa,4.666666666666667);
258 t34 = 1/pow(rhoa,7.333333333333333);
259 t35 = 1/pow(t25,5.626666666666667);
260 t36 = 1/pow(rhoa,3.333333333333333);
261 t37 = 1/pow(t14,3.0);
262 t38 = -3.639457589892744E-6*t1*t24*t18*t26+t19+t17;
263 t39 = 1/
264 pow(t6,3.0);
265 t40 = pow(grada,1.72);
266
267 /* code */
268 res[0] = t13*t20*t21*(-3.639457589892744E-6*t1*t18*pow(t3,
269 2.72)*grada+t19+t17)-1.0*t21*t27*t7-1.333333333333333*t7*t8*
270 t13;
271 res[1] = t20*t21*t13*t32-1.0*t21*t31*t7;
272 res[2] = -1.0*t21*t7*(-19.254544*t11*t12*pow(grada,6.0)/
273 pow(rhoa,10.0)-9.777777777777779*t11*t12*t33*t9-1.798826666666667E-5*
274 t24*t2*t35+0.041653333333333*t9*t33-1.1573333333333333E-5*
275 t24*t36*t26+39.492*t11*t22*t34*t12)-0.444444444444444*t13*
276 t7/pow(rhoa,0.666666666666667)-2.0*t13*t21*pow(t38,2.0)*t39+
277 2.0*t20*t21*t27*t38+2.666666666666667*t20*t8*t13*t38-2.666666666666667*
278 t7*t8*t27+t20*t21*t13*(8.492067709749736E-6*t1*t24*t36*t26+
279 1.3199099526011019E-5*t1*t24*t2*t35+0.07952*grada*t5*t36+0.1704*
280 t9*t15*t33-0.04544*t22*t37*t34);
281 res[3] = -1.0*t21*t7*(14.440908*t11*t12*pow(grada,5.0)/
282 pow(rhoa,9.0)+5.333333333333333*t11*t12*t16*grada+1.84512E-5*
283 t30*t18*t26-0.02272*grada*t16-26.328*t11*t28*t23*t12)-2.0*
284 t13*t21*t32*t38*t39+t20*t21*t27*t32+1.333333333333333*t20*
285 t8*t13*t32-1.333333333333333*t7*t8*t31+t20*t21*t38*t31+t20*
286 t21*t13*(-1.3538782234401007E-5*t1*t30*t18*t26-0.03408*t5*
287 t18-0.10224*grada*t15*t16+0.03408*t28*t37*t23);
288 res[4] = -1.0*t21*t7*(-10.830681*t11*t12*t22/pow(rhoa,
289 8.0)-1.0118400000000001E-5*t40*t2*t26+16.455*t11*t9*t29*t12-
290 2.0*t10*t11*t12+0.00852*t10)-2.0*t13*t21*pow(t32,2.0)*t39+
291 2.0*t20*t21*t31*t32+t20*t21*t13*(7.424493483381199E-6*t1*t40*
292 t2*t26+0.05112*t15*t10-0.02556*t9*t37*t29);
293
294 }
295
296 static void
mpwx_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)297 mpwx_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
298 {
299 real res[5];
300
301 mpwx_second_helper(dp->rhoa, dp->grada, res);
302
303 ds->df1000 += factor*res[0];
304 ds->df0010 += factor*res[1];
305
306 ds->df2000 += factor*res[2];
307 ds->df1010 += factor*res[3];
308 ds->df0020 += factor*res[4];
309
310
311 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
312 fabs(dp->grada-dp->gradb)>1e-13)
313 mpwx_second_helper(dp->rhob, dp->gradb, res);
314 ds->df0100 += factor*res[0];
315 ds->df0001 += factor*res[1];
316
317 ds->df0200 += factor*res[2];
318 ds->df0101 += factor*res[3];
319 ds->df0002 += factor*res[4];
320
321 }
322
323 static void
mpwx_third_helper(real rhoa,real grada,real * res)324 mpwx_third_helper(real rhoa, real grada, real *res)
325 {
326 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
327 real t11, t12, t13, t14, t15, t16, t17, t18;
328 real t19, t20, t21, t22, t23, t24, t25, t26;
329 real t27, t28, t29, t30, t31, t32, t33, t34;
330 real t35, t36, t37, t38, t39, t40, t41, t42;
331 real t43, t44, t45, t46, t47, t48, t49, t50;
332 real t51, t52, t53, t54, t55, t56, t57, t58;
333 real t59, t60, t61, t62, t63;
334
335 t1 = pow(3.141592653589793,0.333333333333333);
336 t2 = 1/pow(rhoa,1.333333333333333);
337 t3 = grada*t2;
338 t4 = pow(t3,3.72);
339 t5 = asinh(t3);
340 t6 = 7.337616108654726E-7*t1*t4+0.02556*grada*t5*t2+1.0;
341 t7 = 1/
342 t6;
343 t8 = pow(rhoa,0.333333333333333);
344 t9 = pow(grada,2.0);
345 t10 = 1/pow(rhoa,2.666666666666667);
346 t11 = 0.00426-0.138888888888889/(pow(3.141592653589793,
347 1.666666666666667)*pow(36.0,0.666666666666667));
348 t12 = 1/pow(2.718281828459045,1.6455*t9*t10);
349 t13 = -1.0*t10*t11*t12*t9-1.0E-6*t4+0.00426*t9*t10;
350 t14 = sqrt(t9*
351 t10+1.0);
352 t15 = 1/t14;
353 t16 = 1/pow(rhoa,3.666666666666667);
354 t17 = -0.03408*t9*t15*t16;
355 t18 = 1/pow(rhoa,2.333333333333333);
356 t19 = -0.03408*grada*t5*t18;
357 t20 = 1/pow(t6,2.0);
358 t21 = pow(rhoa,1.333333333333333);
359 t22 = pow(grada,4.0);
360 t23 = 1/pow(rhoa,6.333333333333333);
361 t24 = pow(grada,3.72);
362 t25 = fabs(rhoa);
363 t26 = 1/pow(t25,3.626666666666667);
364 t27 = 2.666666666666667*t11*t12*t16*t9+4.96E-6*t24*t18*
365 t26-0.01136*t9*t16-4.388*t11*t22*t23*t12;
366 t28 = pow(grada,3.0);
367 t29 = 1/pow(rhoa,5.333333333333333);
368 t30 = pow(grada,2.72);
369 t31 = -2.0*t10*t11*t12*grada-3.72E-6*t30*t2*t26+3.291*
370 t11*t28*t29*t12+0.00852*grada*t10;
371 t32 = 2.729593192419558E-6*t1*t30*t2*t26+0.02556*t5*t2+
372 0.02556*grada*t15*t10;
373 t33 = 1/pow(rhoa,0.666666666666667);
374 t34 = 1/pow(rhoa,4.666666666666667);
375 t35 = pow(grada,6.0);
376 t36 = 1/pow(rhoa,10.0);
377 t37 = 1/pow(rhoa,7.333333333333333);
378 t38 = 1/pow(t25,5.626666666666667);
379 t39 = 1/pow(rhoa,3.333333333333333);
380 t40 = -9.777777777777779*t11*t12*t34*t9-1.798826666666667E-5*
381 t24*t2*t38+0.041653333333333*t9*t34-1.1573333333333333E-5*
382 t24*t39*t26+39.492*t11*t22*t37*t12-19.254544*t11*t35*t36*t12;
383 t41 = 1/
384 pow(t14,3.0);
385 t42 = 8.492067709749736E-6*t1*t24*t39*t26+1.3199099526011019E-5*
386 t1*t24*t2*t38+0.07952*grada*t5*t39+0.1704*t9*t15*t34-0.04544*
387 t22*t41*t37;
388 t43 = -3.639457589892744E-6*t1*t24*t18*t26+t19+t17;
389 t44 = 1/
390 pow(t6,3.0);
391 t45 = pow(t43,2.0);
392 t46 = pow(grada,5.0);
393 t47 = 1/pow(rhoa,9.0);
394 t48 = 5.333333333333333*t11*t12*t16*grada+1.84512E-5*
395 t30*t18*t26-0.02272*grada*t16+14.440908*t11*t46*t47*t12-26.328*
396 t11*t28*t23*t12;
397 t49 = -1.3538782234401007E-5*t1*t30*t18*t26-0.03408*t5*
398 t18-0.10224*grada*t15*t16+0.03408*t28*t41*t23;
399 t50 = 1/pow(rhoa,8.0);
400 t51 = pow(grada,1.72);
401 t52 = -1.0118400000000001E-5*t51*t2*t26-10.830681*t11*
402 t22*t50*t12+16.455*t11*t9*t29*t12-2.0*t10*t11*t12+0.00852*
403 t10;
404 t53 = 7.424493483381199E-6*t1*t51*t2*t26+0.05112*t15*
405 t10-0.02556*t9*t41*t29;
406 t54 = pow(t32,2.0);
407 t55 = 1/pow(rhoa,5.666666666666667);
408 t56 = 1/pow(rhoa,11.0);
409 t57 = 1/pow(rhoa,8.333333333333334);
410 t58 = 1/pow(rhoa,0.333333333333333);
411 t59 = 1/pow(t25,7.626666666666667);
412 t60 = 1/pow(rhoa,4.333333333333333);
413 t61 = 1/pow(t14,5.0);
414 t62 = 1/pow(t6,4.0);
415 t63 = pow(grada,0.72);
416
417 /* code */
418 res[0] = t13*t20*t21*(-3.639457589892744E-6*t1*t18*pow(t3,
419 2.72)*grada+t19+t17)-1.0*t21*t27*t7-1.333333333333333*t7*t8*
420 t13;
421 res[1] = t20*t21*t13*t32-1.0*t21*t31*t7;
422 res[2] = -1.0*t21*t40*t7-2.0*t13*t21*t44*t45+2.0*t20*
423 t21*t27*t43+2.666666666666667*t20*t8*t13*t43+t20*t21*t13*t42-
424 2.666666666666667*t7*t8*t27-0.444444444444444*t7*t33*t13;
425 res[3] = -1.0*t21*t48*t7+t20*t21*t13*t49-2.0*t13*t21*
426 t32*t43*t44+t20*t21*t27*t32+1.333333333333333*t20*t8*t13*t32-
427 1.333333333333333*t7*t8*t31+t20*t21*t43*t31;
428 res[4] = -1.0*t21*t52*t7-2.0*t13*t21*t44*t54+t20*t21*
429 t13*t53+2.0*t20*t21*t31*t32;
430 res[5] = -1.0*t21*t7*(-84.488939072*t11*t12*pow(grada,
431 8.0)/pow(rhoa,13.66666666666667)+45.62962962962963*t11*t12*
432 t55*t9+1.0121398044444446E-4*t24*t58*t59-0.194382222222222*
433 t9*t55+6.595697777777779E-5*t24*t18*t38+3.857777777777778E-5*
434 t24*t60*t26-332.5128888888888*t11*t22*t57*t12+365.836336*t11*
435 t35*t56*t12)+0.296296296296296*t13*t7/pow(rhoa,1.666666666666667)+
436 6.0*t13*t21*pow(t43,3.0)*t62-6.0*t21*t27*t44*t45-8.0*t44*t8*
437 t13*t45-6.0*t13*t21*t42*t43*t44+3.0*t20*t21*t40*t43+8.0*t20*
438 t8*t27*t43+1.333333333333333*t20*t33*t13*t43+3.0*t20*t21*t27*
439 t42+4.0*t20*t8*t13*t42-4.0*t7*t8*t40-1.333333333333333*t7*
440 t33*t27+t20*t21*t13*(-2.8306892365832456E-5*t1*t24*t60*t26-
441 4.83966982620404E-5*t1*t24*t18*t38-7.4266933333022E-5*t1*t24*
442 t58*t59-0.265066666666667*grada*t5*t60-0.901226666666667*t9*
443 t15*t55+0.560426666666667*t22*t41*t57-0.18176*t35*t61*t56);
444 res[6] = -1.0*t21*t7*(63.366704304*t11*t12*pow(grada,
445 7.0)/pow(rhoa,12.66666666666667)-19.55555555555556*t11*t12*
446 t34*grada-6.691635200000001E-5*t30*t2*t38+0.083306666666667*
447 grada*t34-4.30528E-5*t30*t39*t26+190.1466666666666*t11*t28*
448 t37*t12-245.495436*t11*t46*t36*t12)+6.0*t13*t21*t32*t45*t62-
449 4.0*t13*t21*t43*t44*t49+2.0*t20*t21*t27*t49+2.666666666666667*
450 t20*t8*t13*t49-2.666666666666667*t7*t8*t48+2.0*t20*t21*t43*
451 t48-2.0*t21*t31*t44*t45-4.0*t21*t27*t32*t43*t44-2.0*t13*t21*
452 t32*t42*t44-5.333333333333333*t44*t8*t13*t43*t32+t20*t21*t40*
453 t32+2.666666666666667*t20*t8*t27*t32+0.444444444444444*t20*
454 t33*t13*t32+2.666666666666667*t20*t8*t43*t31+t20*t21*t42*t31-
455 0.444444444444444*t7*t33*t31+t20*t21*t13*(3.159049188026902E-5*
456 t1*t30*t39*t26+4.910065023676099E-5*t1*t30*t2*t38+0.07952*
457 t5*t39+0.42032*grada*t15*t34-0.35216*t28*t41*t37+0.13632*t46*
458 t61*t36);
459 res[7] = -1.0*t21*t7*(-47.525028228*t11*t12*t35/pow(rhoa,
460 11.66666666666667)+5.018726400000001E-5*t51*t18*t26+5.333333333333333*
461 t11*t12*t16-0.02272*t16+158.849988*t11*t22*t47*t12-96.536*
462 t11*t9*t23*t12)+6.0*t13*t21*t43*t54*t62-2.0*t21*t27*t44*t54-
463 2.666666666666667*t44*t8*t13*t54-2.0*t13*t21*t43*t44*t53+t20*
464 t21*t27*t53+1.333333333333333*t20*t8*t13*t53-1.333333333333333*
465 t7*t8*t52+t20*t21*t43*t52-4.0*t13*t21*t32*t44*t49+2.0*t20*
466 t21*t31*t49+2.0*t20*t21*t32*t48-4.0*t21*t31*t32*t43*t44+2.666666666666667*
467 t20*t8*t31*t32+t20*t21*t13*(-3.682548767757074E-5*t1*t51*t18*
468 t26-0.13632*t15*t16+0.20448*t9*t41*t23-0.10224*t22*t61*t47);
469 res[8] = -1.0*t21*t7*(35.643771171*t11*t12*t46/pow(rhoa,
470 10.66666666666667)-1.7403648000000005E-5*t63*t2*t26-97.47612899999999*
471 t11*t28*t50*t12+39.492*t11*grada*t29*t12)+6.0*t13*t21*pow(t32,
472 3.0)*t62-6.0*t21*t31*t44*t54-6.0*t13*t21*t32*t44*t53+3.0*t20*
473 t21*t31*t53+3.0*t20*t21*t32*t52+t20*t21*t13*(1.2770128791415663E-5*
474 t1*t63*t2*t26-0.10224*grada*t41*t29+0.07668*t28*t61*t50);
475
476 }
477
478 static void
mpwx_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)479 mpwx_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
480 {
481 real res[9];
482
483 mpwx_third_helper(dp->rhoa, dp->grada, res);
484
485 ds->df1000 += factor*res[0];
486 ds->df0010 += factor*res[1];
487
488 ds->df2000 += factor*res[2];
489 ds->df1010 += factor*res[3];
490 ds->df0020 += factor*res[4];
491
492 ds->df3000 += factor*res[5];
493 ds->df2010 += factor*res[6];
494 ds->df1020 += factor*res[7];
495 ds->df0030 += factor*res[8];
496
497
498 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
499 fabs(dp->grada-dp->gradb)>1e-13)
500 mpwx_third_helper(dp->rhob, dp->gradb, res);
501
502 ds->df0100 += factor*res[0];
503 ds->df0001 += factor*res[1];
504
505 ds->df0200 += factor*res[2];
506 ds->df0101 += factor*res[3];
507 ds->df0002 += factor*res[4];
508
509 ds->df0300 += factor*res[5];
510 ds->df0201 += factor*res[6];
511 ds->df0102 += factor*res[7];
512 ds->df0003 += factor*res[8];
513
514 }
515
516 static void
mpwx_fourth_helper(real rhoa,real grada,real * res)517 mpwx_fourth_helper(real rhoa, real grada, real *res)
518 {
519 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
520 real t11, t12, t13, t14, t15, t16, t17, t18;
521 real t19, t20, t21, t22, t23, t24, t25, t26;
522 real t27, t28, t29, t30, t31, t32, t33, t34;
523 real t35, t36, t37, t38, t39, t40, t41, t42;
524 real t43, t44, t45, t46, t47, t48, t49, t50;
525 real t51, t52, t53, t54, t55, t56, t57, t58;
526 real t59, t60, t61, t62, t63, t64, t65, t66;
527 real t67, t68, t69, t70, t71, t72, t73, t74;
528 real t75, t76, t77, t78, t79, t80, t81, t82;
529 real t83, t84, t85, t86, t87, t88, t89;
530
531 t1 = pow(3.141592653589793,0.333333333333333);
532 t2 = 1/pow(rhoa,1.333333333333333);
533 t3 = grada*t2;
534 t4 = pow(t3,3.72);
535 t5 = asinh(t3);
536 t6 = 7.337616108654726E-7*t1*t4+0.02556*grada*t5*t2+1.0;
537 t7 = 1/
538 t6;
539 t8 = pow(rhoa,0.333333333333333);
540 t9 = pow(grada,2.0);
541 t10 = 1/pow(rhoa,2.666666666666667);
542 t11 = 0.00426-0.138888888888889/(pow(3.141592653589793,
543 1.666666666666667)*pow(36.0,0.666666666666667));
544 t12 = 1/pow(2.718281828459045,1.6455*t9*t10);
545 t13 = -1.0*t10*t11*t12*t9-1.0E-6*t4+0.00426*t9*t10;
546 t14 = sqrt(t9*
547 t10+1.0);
548 t15 = 1/t14;
549 t16 = 1/pow(rhoa,3.666666666666667);
550 t17 = -0.03408*t9*t15*t16;
551 t18 = 1/pow(rhoa,2.333333333333333);
552 t19 = -0.03408*grada*t5*t18;
553 t20 = 1/pow(t6,2.0);
554 t21 = pow(rhoa,1.333333333333333);
555 t22 = pow(grada,4.0);
556 t23 = 1/pow(rhoa,6.333333333333333);
557 t24 = pow(grada,3.72);
558 t25 = fabs(rhoa);
559 t26 = 1/pow(t25,3.626666666666667);
560 t27 = 2.666666666666667*t11*t12*t16*t9+4.96E-6*t24*t18*
561 t26-0.01136*t9*t16-4.388*t11*t22*t23*t12;
562 t28 = pow(grada,3.0);
563 t29 = 1/pow(rhoa,5.333333333333333);
564 t30 = pow(grada,2.72);
565 t31 = -2.0*t10*t11*t12*grada-3.72E-6*t30*t2*t26+3.291*
566 t11*t28*t29*t12+0.00852*grada*t10;
567 t32 = 2.729593192419558E-6*t1*t30*t2*t26+0.02556*t5*t2+
568 0.02556*grada*t15*t10;
569 t33 = 1/pow(rhoa,0.666666666666667);
570 t34 = 1/pow(rhoa,4.666666666666667);
571 t35 = pow(grada,6.0);
572 t36 = 1/pow(rhoa,10.0);
573 t37 = 1/pow(rhoa,7.333333333333333);
574 t38 = 1/pow(t25,5.626666666666667);
575 t39 = 1/pow(rhoa,3.333333333333333);
576 t40 = -9.777777777777779*t11*t12*t34*t9-1.798826666666667E-5*
577 t24*t2*t38+0.041653333333333*t9*t34-1.1573333333333333E-5*
578 t24*t39*t26+39.492*t11*t22*t37*t12-19.254544*t11*t35*t36*t12;
579 t41 = 1/
580 pow(t14,3.0);
581 t42 = 8.492067709749736E-6*t1*t24*t39*t26+1.3199099526011019E-5*
582 t1*t24*t2*t38+0.07952*grada*t5*t39+0.1704*t9*t15*t34-0.04544*
583 t22*t41*t37;
584 t43 = -3.639457589892744E-6*t1*t24*t18*t26+t19+t17;
585 t44 = 1/
586 pow(t6,3.0);
587 t45 = pow(t43,2.0);
588 t46 = pow(grada,5.0);
589 t47 = 1/pow(rhoa,9.0);
590 t48 = 5.333333333333333*t11*t12*t16*grada+1.84512E-5*
591 t30*t18*t26-0.02272*grada*t16+14.440908*t11*t46*t47*t12-26.328*
592 t11*t28*t23*t12;
593 t49 = -1.3538782234401007E-5*t1*t30*t18*t26-0.03408*t5*
594 t18-0.10224*grada*t15*t16+0.03408*t28*t41*t23;
595 t50 = 1/pow(rhoa,8.0);
596 t51 = pow(grada,1.72);
597 t52 = -1.0118400000000001E-5*t51*t2*t26-10.830681*t11*
598 t22*t50*t12+16.455*t11*t9*t29*t12-2.0*t10*t11*t12+0.00852*
599 t10;
600 t53 = 7.424493483381199E-6*t1*t51*t2*t26+0.05112*t15*
601 t10-0.02556*t9*t41*t29;
602 t54 = pow(t32,2.0);
603 t55 = 1/pow(rhoa,1.666666666666667);
604 t56 = 1/pow(rhoa,5.666666666666667);
605 t57 = pow(grada,8.0);
606 t58 = 1/pow(rhoa,13.66666666666667);
607 t59 = 1/pow(rhoa,11.0);
608 t60 = 1/pow(rhoa,8.333333333333334);
609 t61 = 1/pow(rhoa,0.333333333333333);
610 t62 = 1/pow(t25,7.626666666666667);
611 t63 = 1/pow(rhoa,4.333333333333333);
612 t64 = 45.62962962962963*t11*t12*t56*t9+1.0121398044444446E-4*
613 t24*t61*t62-0.194382222222222*t9*t56+6.595697777777779E-5*
614 t24*t18*t38+3.857777777777778E-5*t24*t63*t26-332.5128888888888*
615 t11*t22*t60*t12+365.836336*t11*t35*t59*t12-84.488939072*t11*
616 t57*t58*t12;
617 t65 = 1/pow(t14,5.0);
618 t66 = -2.8306892365832456E-5*t1*t24*t63*t26-4.83966982620404E-5*
619 t1*t24*t18*t38-7.4266933333022E-5*t1*t24*t61*t62-0.265066666666667*
620 grada*t5*t63-0.901226666666667*t9*t15*t56+0.560426666666667*
621 t22*t41*t60-0.18176*t35*t65*t59;
622 t67 = 1/pow(t6,4.0);
623 t68 = pow(t43,3.0);
624 t69 = pow(grada,7.0);
625 t70 = 1/pow(rhoa,12.66666666666667);
626 t71 = -19.55555555555556*t11*t12*t34*grada-6.691635200000001E-5*
627 t30*t2*t38+0.083306666666667*grada*t34-4.30528E-5*t30*t39*
628 t26+63.366704304*t11*t69*t70*t12+190.1466666666666*t11*t28*
629 t37*t12-245.495436*t11*t46*t36*t12;
630 t72 = 3.159049188026902E-5*t1*t30*t39*t26+4.910065023676099E-5*
631 t1*t30*t2*t38+0.07952*t5*t39+0.42032*grada*t15*t34-0.35216*
632 t28*t41*t37+0.13632*t46*t65*t36;
633 t73 = 1/pow(rhoa,11.66666666666667);
634 t74 = 5.018726400000001E-5*t51*t18*t26+5.333333333333333*
635 t11*t12*t16-0.02272*t16-47.525028228*t11*t35*t73*t12+158.849988*
636 t11*t22*t47*t12-96.536*t11*t9*t23*t12;
637 t75 = -3.682548767757074E-5*t1*t51*t18*t26-0.13632*t15*
638 t16+0.20448*t9*t41*t23-0.10224*t22*t65*t47;
639 t76 = 1/pow(rhoa,10.66666666666667);
640 t77 = pow(grada,0.72);
641 t78 = -1.7403648000000005E-5*t77*t2*t26+39.492*t11*grada*
642 t29*t12-97.47612899999999*t11*t28*t50*t12+35.643771171*t11*
643 t46*t76*t12;
644 t79 = 1.2770128791415663E-5*t1*t77*t2*t26-0.10224*grada*
645 t41*t29+0.07668*t28*t65*t50;
646 t80 = pow(t32,3.0);
647 t81 = 1/pow(rhoa,6.666666666666667);
648 t82 = 1/pow(rhoa,14.66666666666667);
649 t83 = 1/pow(rhoa,12.0);
650 t84 = 1/pow(rhoa,9.333333333333334);
651 t85 = pow(rhoa,0.666666666666667);
652 t86 = 1/pow(t25,9.626666666666667);
653 t87 = 1/pow(t14,7.0);
654 t88 = 1/pow(t6,5.0);
655 t89 = 1/pow(grada,0.28);
656
657 /* code */
658 res[0] = t13*t20*t21*(-3.639457589892744E-6*t1*t18*pow(t3,
659 2.72)*grada+t19+t17)-1.0*t21*t27*t7-1.333333333333333*t7*t8*
660 t13;
661 res[1] = t20*t21*t13*t32-1.0*t21*t31*t7;
662 res[2] = -1.0*t21*t40*t7-2.0*t13*t21*t44*t45+2.0*t20*
663 t21*t27*t43+2.666666666666667*t20*t8*t13*t43+t20*t21*t13*t42-
664 2.666666666666667*t7*t8*t27-0.444444444444444*t7*t33*t13;
665 res[3] = -1.0*t21*t48*t7+t20*t21*t13*t49-2.0*t13*t21*
666 t32*t43*t44+t20*t21*t27*t32+1.333333333333333*t20*t8*t13*t32-
667 1.333333333333333*t7*t8*t31+t20*t21*t43*t31;
668 res[4] = -1.0*t21*t52*t7-2.0*t13*t21*t44*t54+t20*t21*
669 t13*t53+2.0*t20*t21*t31*t32;
670 res[5] = -1.0*t21*t64*t7+6.0*t13*t21*t67*t68+t20*t21*
671 t13*t66-6.0*t21*t27*t44*t45-8.0*t44*t8*t13*t45-6.0*t13*t21*
672 t42*t43*t44+3.0*t20*t21*t40*t43+8.0*t20*t8*t27*t43+1.333333333333333*
673 t20*t33*t13*t43+3.0*t20*t21*t27*t42+4.0*t20*t8*t13*t42-4.0*
674 t7*t8*t40-1.333333333333333*t7*t33*t27+0.296296296296296*t7*
675 t55*t13;
676 res[6] = t20*t21*t13*t72-1.0*t21*t7*t71+6.0*t13*t21*t32*
677 t45*t67-4.0*t13*t21*t43*t44*t49+2.0*t20*t21*t27*t49+2.666666666666667*
678 t20*t8*t13*t49-2.666666666666667*t7*t8*t48+2.0*t20*t21*t43*
679 t48-2.0*t21*t31*t44*t45-4.0*t21*t27*t32*t43*t44-2.0*t13*t21*
680 t32*t42*t44-5.333333333333333*t44*t8*t13*t43*t32+t20*t21*t40*
681 t32+2.666666666666667*t20*t8*t27*t32+0.444444444444444*t20*
682 t33*t13*t32+2.666666666666667*t20*t8*t43*t31+t20*t21*t42*t31-
683 0.444444444444444*t7*t33*t31;
684 res[7] = t20*t21*t13*t75-1.0*t21*t7*t74+6.0*t13*t21*t43*
685 t54*t67-2.0*t21*t27*t44*t54-2.666666666666667*t44*t8*t13*t54-
686 2.0*t13*t21*t43*t44*t53+t20*t21*t27*t53+1.333333333333333*
687 t20*t8*t13*t53-1.333333333333333*t7*t8*t52+t20*t21*t43*t52-
688 4.0*t13*t21*t32*t44*t49+2.0*t20*t21*t31*t49+2.0*t20*t21*t32*
689 t48-4.0*t21*t31*t32*t43*t44+2.666666666666667*t20*t8*t31*t32;
690 res[8] = 6.0*t13*t21*t67*t80+t20*t21*t13*t79-1.0*t21*
691 t7*t78-6.0*t21*t31*t44*t54-6.0*t13*t21*t32*t44*t53+3.0*t20*
692 t21*t31*t53+3.0*t20*t21*t32*t52;
693 res[9] = -1.0*t21*t7*(-370.737464647936*t11*t12*pow(grada,
694 10.0)/pow(rhoa,17.33333333333333)-258.5679012345679*t11*t12*
695 t81*t9-7.719252908562964E-4*t24*t85*t86+1.101499259259259*
696 t9*t81-4.048559217777778E-4*t24*t2*t62-2.938083555555556E-4*
697 t24*t39*t38-1.6717037037037035E-4*t24*t29*t26+2971.163555555555*
698 t11*t22*t84*t12-5483.266252444444*t11*t35*t83*t12+2759.972009685333*
699 t11*t57*t82*t12)-0.493827160493827*t13*t7/pow(rhoa,2.666666666666667)-
700 24.0*t13*t21*pow(t43,4.0)*t88+24.0*t21*t27*t67*t68+32.0*t67*
701 t8*t13*t68+36.0*t13*t21*t42*t45*t67-8.0*t13*t21*t43*t44*t66+
702 4.0*t20*t21*t27*t66+5.333333333333333*t20*t8*t13*t66-5.333333333333333*
703 t7*t8*t64+4.0*t20*t21*t43*t64-12.0*t21*t40*t44*t45-32.0*t44*
704 t8*t27*t45-5.333333333333332*t44*t33*t13*t45-24.0*t21*t27*
705 t42*t43*t44-6.0*t13*t21*pow(t42,2.0)*t44-32.0*t44*t8*t13*t42*
706 t43+16.0*t20*t8*t40*t43+5.333333333333332*t20*t33*t27*t43-
707 1.185185185185185*t20*t55*t13*t43+6.0*t20*t21*t40*t42+2.666666666666666*
708 t20*t33*t13*t42-2.666666666666666*t7*t33*t40+1.185185185185185*
709 t7*t55*t27+16.0*t20*t8*t42*t27+t20*t21*t13*(1.2266320025194064E-4*
710 t1*t24*t29*t26+2.1558529225818E-4*t1*t24*t39*t38+2.97067733332088E-4*
711 t1*t24*t2*t62+5.664091448865145E-4*t1*t24*t85*t86+1.148622222222222*
712 grada*t5*t29+5.460373333333334*t9*t15*t81-5.871857777777778*
713 t22*t41*t84+4.241066666666667*t35*t65*t83-1.211733333333333*
714 t57*t87*t82);
715 res[10] = -1.0*t21*t7*(278.053098485952*t11*t12*pow(grada,
716 9.0)/pow(rhoa,16.33333333333333)+91.25925925925925*t11*t12*
717 t56*grada+3.765160072533334E-4*t30*t61*t62-0.388764444444444*
718 grada*t56+2.453599573333334E-4*t30*t18*t38+1.4350933333333335E-4*
719 t30*t63*t26-1480.218666666666*t11*t28*t60*t12+3289.317933333333*
720 t11*t46*t59*t12-1879.878894352*t11*t69*t58*t12)-24.0*t13*t21*
721 t32*t68*t88-6.0*t13*t21*t43*t44*t72+3.0*t20*t21*t27*t72+4.0*
722 t20*t8*t13*t72-4.0*t7*t8*t71+3.0*t20*t21*t43*t71+6.0*t21*t31*
723 t67*t68+18.0*t13*t21*t45*t49*t67+18.0*t21*t27*t32*t45*t67+
724 18.0*t13*t21*t32*t42*t43*t67-2.0*t13*t21*t32*t44*t66-12.0*
725 t21*t27*t43*t44*t49-6.0*t13*t21*t42*t44*t49+3.0*t20*t21*t40*
726 t49+1.333333333333333*t20*t33*t13*t49-6.0*t21*t44*t45*t48+
727 3.0*t20*t21*t42*t48-1.333333333333333*t7*t33*t48-6.0*t21*t31*
728 t42*t43*t44-6.0*t21*t32*t40*t43*t44-6.0*t21*t27*t32*t42*t44-
729 16.0*t44*t8*t13*t49*t43+8.0*t20*t8*t48*t43+t20*t21*t64*t32+
730 24.0*t67*t8*t13*t45*t32-16.0*t44*t8*t27*t43*t32-2.666666666666666*
731 t44*t33*t13*t43*t32-8.0*t44*t8*t13*t42*t32+4.0*t20*t8*t40*
732 t32+1.333333333333333*t20*t33*t27*t32-0.296296296296296*t20*
733 t55*t13*t32+t20*t21*t66*t31+0.296296296296296*t7*t55*t31-8.0*
734 t44*t8*t45*t31+1.333333333333333*t20*t33*t43*t31+4.0*t20*t8*
735 t42*t31+8.0*t20*t8*t49*t27+t20*t21*t13*(-1.0530163960089674E-4*
736 t1*t30*t63*t26-1.800357175347903E-4*t1*t30*t18*t38-2.762729919988419E-4*
737 t1*t30*t61*t62-0.265066666666667*t5*t63-2.06752*grada*t15*
738 t56+3.142933333333333*t28*t41*t60-2.77184*t46*t65*t59+0.9088*
739 t69*t87*t58);
740 res[11] = -1.0*t21*t7*(-208.539823864464*t11*t12*t57/
741 pow(rhoa,15.33333333333333)-1.8201247744000003E-4*t51*t2*t38-
742 19.55555555555556*t11*t12*t34+0.083306666666667*t34-1.1710361600000001E-4*
743 t51*t39*t26+1251.492410004*t11*t35*t70*t12+634.7973333333333*
744 t11*t9*t37*t12-1853.24986*t11*t22*t36*t12)-24.0*t13*t21*t45*
745 t54*t88-4.0*t13*t21*t43*t44*t75+2.0*t20*t21*t27*t75+2.666666666666667*
746 t20*t8*t13*t75-2.666666666666667*t7*t8*t74+2.0*t20*t21*t43*
747 t74-4.0*t13*t21*t32*t44*t72+2.0*t20*t21*t31*t72+2.0*t20*t21*
748 t32*t71+12.0*t21*t27*t43*t54*t67+6.0*t13*t21*t42*t54*t67+6.0*
749 t13*t21*t45*t53*t67+24.0*t13*t21*t32*t43*t49*t67+12.0*t21*
750 t31*t32*t45*t67-2.0*t21*t40*t44*t54+16.0*t67*t8*t13*t43*t54-
751 5.333333333333333*t44*t8*t27*t54-0.888888888888889*t44*t33*
752 t13*t54-4.0*t21*t27*t43*t44*t53-2.0*t13*t21*t42*t44*t53-5.333333333333333*
753 t44*t8*t13*t43*t53+t20*t21*t40*t53+2.666666666666667*t20*t8*
754 t27*t53+0.444444444444444*t20*t33*t13*t53-2.0*t21*t44*t45*
755 t52+2.666666666666667*t20*t8*t43*t52+t20*t21*t42*t52-0.444444444444444*
756 t7*t33*t52-4.0*t13*t21*t44*pow(t49,2.0)+4.0*t20*t21*t48*t49-
757 8.0*t21*t31*t43*t44*t49-8.0*t21*t27*t32*t44*t49-8.0*t21*t32*
758 t43*t44*t48-4.0*t21*t31*t32*t42*t44-10.66666666666667*t44*
759 t8*t13*t49*t32+5.333333333333333*t20*t8*t48*t32-10.66666666666667*
760 t44*t8*t43*t31*t32+0.888888888888889*t20*t33*t31*t32+5.333333333333333*
761 t20*t8*t49*t31+t20*t21*t13*(8.592613791433175E-5*t1*t51*t39*
762 t26+1.3355376864398991E-4*t1*t51*t2*t38+0.49984*t15*t34-1.4768*
763 t9*t41*t37+1.73808*t22*t65*t36-0.6816*t35*t87*t70);
764 res[12] = -1.0*t21*t7*(156.404867898348*t11*t12*t69/pow(rhoa,
765 14.33333333333333)+8.632209408000002E-5*t77*t18*t26-807.9254798759999*
766 t11*t46*t73*t12+953.099928*t11*t28*t47*t12-210.624*t11*grada*
767 t23*t12)-24.0*t13*t21*t43*t80*t88+6.0*t21*t27*t67*t80+8.0*
768 t67*t8*t13*t80-2.0*t13*t21*t43*t44*t79+t20*t21*t27*t79+1.333333333333333*
769 t20*t8*t13*t79-1.333333333333333*t7*t8*t78+t20*t21*t43*t78-
770 6.0*t13*t21*t32*t44*t75+3.0*t20*t21*t31*t75+3.0*t20*t21*t32*
771 t74+18.0*t13*t21*t49*t54*t67+18.0*t21*t31*t43*t54*t67+18.0*
772 t13*t21*t32*t43*t53*t67-6.0*t21*t44*t48*t54-8.0*t44*t8*t31*
773 t54-6.0*t13*t21*t44*t49*t53+3.0*t20*t21*t48*t53-6.0*t21*t31*
774 t43*t44*t53-6.0*t21*t27*t32*t44*t53+3.0*t20*t21*t49*t52-6.0*
775 t21*t32*t43*t44*t52-12.0*t21*t31*t32*t44*t49-8.0*t44*t8*t13*
776 t53*t32+4.0*t20*t8*t52*t32+4.0*t20*t8*t53*t31+t20*t21*t13*
777 (-6.333983880542168E-5*t1*t77*t18*t26+0.54528*grada*t41*t23-
778 1.0224*t28*t65*t47+0.5112*t46*t87*t73);
779 res[13] = -1.0*t21*t7*(-117.303650923761*t11*t12*t35/
780 pow(rhoa,13.33333333333333)-1.2530626560000006E-5*t89*t2*t26+
781 499.0127963939999*t11*t22*t76*t12-422.3965589999999*t11*t9*
782 t50*t12+39.492*t11*t29*t12)-24.0*t13*t21*pow(t32,4.0)*t88+
783 24.0*t21*t31*t67*t80-8.0*t13*t21*t32*t44*t79+4.0*t20*t21*t31*
784 t79+4.0*t20*t21*t32*t78+36.0*t13*t21*t53*t54*t67-12.0*t21*
785 t44*t52*t54-6.0*t13*t21*t44*pow(t53,2.0)+6.0*t20*t21*t52*t53-
786 24.0*t21*t31*t32*t44*t53+t20*t21*t13*(9.194492729819279E-6*
787 t1*t89*t2*t26-0.10224*t41*t29+0.53676*t9*t65*t50-0.3834*t22*
788 t87*t76);
789
790 }
791
792 static void
mpwx_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)793 mpwx_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
794 {
795 real res[14];
796
797 mpwx_fourth_helper(dp->rhoa, dp->grada, res);
798
799 ds->df1000 += factor*res[0];
800 ds->df0010 += factor*res[1];
801
802 ds->df2000 += factor*res[2];
803 ds->df1010 += factor*res[3];
804 ds->df0020 += factor*res[4];
805
806 ds->df3000 += factor*res[5];
807 ds->df2010 += factor*res[6];
808 ds->df1020 += factor*res[7];
809 ds->df0030 += factor*res[8];
810
811 ds->df4000 += factor*res[9];
812 ds->df3010 += factor*res[10];
813 ds->df2020 += factor*res[11];
814 ds->df1030 += factor*res[12];
815 ds->df0040 += factor*res[13];
816
817
818 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
819 fabs(dp->grada-dp->gradb)>1e-13)
820 mpwx_fourth_helper(dp->rhob, dp->gradb, res);
821
822 ds->df0100 += factor*res[0];
823 ds->df0001 += factor*res[1];
824
825 ds->df0200 += factor*res[2];
826 ds->df0101 += factor*res[3];
827 ds->df0002 += factor*res[4];
828
829 ds->df0300 += factor*res[5];
830 ds->df0201 += factor*res[6];
831 ds->df0102 += factor*res[7];
832 ds->df0003 += factor*res[8];
833
834 ds->df0400 += factor*res[9];
835 ds->df0301 += factor*res[10];
836 ds->df0202 += factor*res[11];
837 ds->df0103 += factor*res[12];
838 ds->df0004 += factor*res[13];
839
840 }
841