1 /*
2
3
4 !
5 ! Dalton, a molecular electronic structure program
6 ! Copyright (C) by the authors of Dalton.
7 !
8 ! This program is free software; you can redistribute it and/or
9 ! modify it under the terms of the GNU Lesser General Public
10 ! License version 2.1 as published by the Free Software Foundation.
11 !
12 ! This program is distributed in the hope that it will be useful,
13 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
14 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 ! Lesser General Public License for more details.
16 !
17 ! If a copy of the GNU LGPL v2.1 was not distributed with this
18 ! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
19 !
20
21 !
22
23 */
24 /*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
25 /* fun-lg93.c:
26
27 LG93 functional: gradient correction GGA exchange functional,
28 with the full LG93 exchange actually SLATER + LG93x.
29
30 Implemented by David Wilson (david.wilson@latrobe.edu.au) Jun 2005
31
32 Automatically generated code implementing LG93 functional and
33 its derivatives. It is generated by func-codegen.pl being a part of
34 a "Automatic code generation framework for analytical functional
35 derivative evaluation", Pawel Salek, 2005
36
37 This functional is connected by making following changes:
38 1. add "extern Functional lg93Functional;" to 'functionals.h'
39 2. add "&lg93Functional," to 'functionals.c'
40 3. add "fun-lg93.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.
41
42 This functional has been generated from following input:
43 ------ cut here -------
44 rho: rhoa + rhob;
45 grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
46 zeta: (rhoa-rhob)/(rhoa+rhob);
47
48 ad: 1.0*10^(-8);
49 a4: 29.790;
50 a6: 22.417;
51 a8: 12.119;
52 a10: 1570.1;
53 a12: 55.944;
54 b: 0.024974;
55 C0: (ad+0.1234)/b;
56
57 xa: grada/(rhoa^(4/3));
58 xb: gradb/(rhob^(4/3));
59
60 sa: 0.5*((3*%PI^2)^(-1/3))*xa;
61 sb: 0.5*((3*%PI^2)^(-1/3))*xb;
62
63 F1a: (1+C0*sa^2+a4*sa^4+a6*sa^6+a8*sa^8+a10*sa^10+a12*sa^12)^b;
64 F1b: (1+C0*sb^2+a4*sb^4+a6*sb^6+a8*sb^8+a10*sb^10+a12*sb^12)^b;
65 F2a: 1+ad*sa^2;
66 F2b: 1+ad*sb^2;
67 Exa: rhoa^(4/3)*F1a/F2a;
68 Exb: rhob^(4/3)*F1b/F2b;
69
70 K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);
71
72
73 ------ cut here -------
74 */
75
76
77 /* strictly conform to XOPEN ANSI C standard */
78 #if !defined(SYS_DEC)
79 /* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
80 * is not specified. */
81 #define _XOPEN_SOURCE 500
82 #define _XOPEN_SOURCE_EXTENDED 1
83 #endif
84 #include <math.h>
85 #include <stddef.h>
86 #include "general.h"
87
88 #define __CVERSION__
89
90 #include "functionals.h"
91
92 /* INTERFACE PART */
lg93_isgga(void)93 static integer lg93_isgga(void) { return 1; } /* FIXME: detect! */
94 static integer lg93_read(const char *conf_line);
95 static real lg93_energy(const FunDensProp* dp);
96 static void lg93_first(FunFirstFuncDrv *ds, real factor,
97 const FunDensProp* dp);
98 static void lg93_second(FunSecondFuncDrv *ds, real factor,
99 const FunDensProp* dp);
100 static void lg93_third(FunThirdFuncDrv *ds, real factor,
101 const FunDensProp* dp);
102 static void lg93_fourth(FunFourthFuncDrv *ds, real factor,
103 const FunDensProp* dp);
104
105 Functional LG93xFunctional = {
106 "LG93x", /* name */
107 lg93_isgga, /* gga-corrected */
108 1,
109 lg93_read,
110 NULL,
111 lg93_energy,
112 lg93_first,
113 lg93_second,
114 lg93_third,
115 lg93_fourth
116 };
117
118 /* IMPLEMENTATION PART */
119 static integer
lg93_read(const char * conf_line)120 lg93_read(const char *conf_line)
121 {
122 fun_set_hf_weight(0);
123 return 1;
124 }
125
126 static real
lg93_energy(const FunDensProp * dp)127 lg93_energy(const FunDensProp *dp)
128 {
129 real res;
130 real rhoa = dp->rhoa, rhob = dp->rhob;
131 real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;
132
133 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
134 real t11, t12;
135
136 t1 = 1/pow(3.0,0.666666666666667);
137 t2 = 1/pow(3.141592653589793,1.333333333333333);
138 t3 = pow(grada,2.0);
139 t4 = 1/pow(rhoa,2.666666666666667);
140 t5 = 1/pow(3.141592653589793,8.0);
141 t6 = 1/pow(3.0,0.333333333333333);
142 t7 = 1/pow(3.141592653589793,6.666666666666667);
143 t8 = 1/pow(3.141592653589793,5.333333333333333);
144 t9 = 1/pow(3.141592653589793,4.0);
145 t10 = 1/pow(3.141592653589793,2.666666666666667);
146 t11 = pow(gradb,2.0);
147 t12 = 1/pow(rhob,2.666666666666667);
148
149 /* code */
150 res = pow(rhob,1.333333333333333)*pow(1.6861979166666666E-4*
151 t5*pow(gradb,12.0)/pow(rhob,16.0)+0.056788917824074*t6*t7*
152 pow(gradb,10.0)/pow(rhob,13.33333333333333)+0.005259982638889*
153 t1*t8*pow(gradb,8.0)/pow(rhob,10.66666666666667)+0.038918402777778*
154 t9*pow(gradb,6.0)/pow(rhob,8.0)+0.620625*t10*t6*pow(gradb,
155 4.0)/pow(rhob,5.333333333333333)+1.235284796188036*t1*t2*t11*
156 t12+1.0,0.024974)/(2.5E-9*t1*t2*t11*t12+1.0)+pow(rhoa,1.333333333333333)*
157 pow(1.6861979166666666E-4*t5*pow(grada,12.0)/pow(rhoa,16.0)+
158 0.056788917824074*t6*t7*pow(grada,10.0)/pow(rhoa,13.33333333333333)+
159 0.005259982638889*t1*t8*pow(grada,8.0)/pow(rhoa,10.66666666666667)+
160 0.038918402777778*t9*pow(grada,6.0)/pow(rhoa,8.0)+0.620625*
161 t10*t6*pow(grada,4.0)/pow(rhoa,5.333333333333333)+1.235284796188036*
162 t1*t2*t3*t4+1.0,0.024974)/(2.5E-9*t1*t2*t3*t4+1.0);
163
164 return res;
165 }
166
167 static void
lg93_first_helper(real rhoa,real grada,real * res)168 lg93_first_helper(real rhoa, real grada, real *res)
169 { real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
170 real t11, t12, t13, t14, t15, t16, t17, t18;
171 real t19, t20, t21, t22, t23, t24, t25, t26;
172 real t27;
173
174 t1 = 1/pow(3.0,0.666666666666667);
175 t2 = 1/pow(3.141592653589793,1.333333333333333);
176 t3 = pow(grada,2.0);
177 t4 = 1/pow(rhoa,2.666666666666667);
178 t5 = 2.5E-9*t1*t2*t3*t4+1.0;
179 t6 = 1/pow(t5,2.0);
180 t7 = 1/pow(3.141592653589793,8.0);
181 t8 = pow(grada,12.0);
182 t9 = 1/pow(rhoa,16.0);
183 t10 = 1/pow(3.0,0.333333333333333);
184 t11 = 1/pow(3.141592653589793,6.666666666666667);
185 t12 = pow(grada,10.0);
186 t13 = 1/pow(rhoa,13.33333333333333);
187 t14 = 1/pow(3.141592653589793,5.333333333333333);
188 t15 = pow(grada,8.0);
189 t16 = 1/pow(rhoa,10.66666666666667);
190 t17 = 1/pow(3.141592653589793,4.0);
191 t18 = pow(grada,6.0);
192 t19 = 1/pow(rhoa,8.0);
193 t20 = 1/pow(3.141592653589793,2.666666666666667);
194 t21 = pow(grada,4.0);
195 t22 = 1/pow(rhoa,5.333333333333333);
196 t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
197 t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
198 t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
199 t10*t11*t12*t13+1.0;
200 t24 = pow(t23,0.024974);
201 t25 = 1/t5;
202 t26 = 1/pow(t23,0.975026);
203 t27 = pow(rhoa,1.333333333333333);
204
205 /* code */
206 res[0] = 0.024974*t25*t26*t27*(-0.002697916666667*t7*
207 t8/pow(rhoa,17.0)-0.757185570987654*t10*t11*t12/pow(rhoa,14.33333333333333)-
208 0.056106481481481*t1*t14*t15/pow(rhoa,11.66666666666667)-0.311347222222222*
209 t17*t18/pow(rhoa,9.0)-3.31*t10*t20*t21/pow(rhoa,6.333333333333333)-
210 3.294092789834761*t1*t2*t3/pow(rhoa,3.666666666666667))+6.666666666666667E-9*
211 t1*t2*t24*t3*t6/pow(rhoa,2.333333333333333)+1.333333333333333*
212 t24*t25*pow(rhoa,0.333333333333333);
213 res[1] = 0.024974*t25*t26*t27*(0.0020234375*t7*t9*pow(grada,
214 11.0)+0.567889178240741*t10*t11*t13*pow(grada,9.0)+0.042079861111111*
215 t1*t14*t16*pow(grada,7.0)+0.233510416666667*t17*t19*pow(grada,
216 5.0)+2.4825*t10*t20*t22*pow(grada,3.0)+2.470569592376071*t1*
217 t2*grada*t4)-5.0E-9*t1*t2*t24*t6*grada/t27;
218 }
219
220 static void
lg93_first(FunFirstFuncDrv * ds,real factor,const FunDensProp * dp)221 lg93_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
222 {
223 real res[2];
224
225 lg93_first_helper(dp->rhoa, dp->grada, res);
226 /* Final assignment */
227 ds->df1000 += factor*res[0];
228 ds->df0010 += factor*res[1];
229
230
231 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
232 fabs(dp->grada-dp->gradb)>1e-13)
233 lg93_first_helper(dp->rhob, dp->gradb, res);
234 ds->df0100 += factor*res[0];
235 ds->df0001 += factor*res[1];
236
237 }
238
239 static void
lg93_second_helper(real rhoa,real grada,real * res)240 lg93_second_helper(real rhoa, real grada, real *res)
241 {
242 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
243 real t11, t12, t13, t14, t15, t16, t17, t18;
244 real t19, t20, t21, t22, t23, t24, t25, t26;
245 real t27, t28, t29, t30, t31, t32, t33, t34;
246 real t35, t36, t37, t38, t39, t40, t41, t42;
247 real t43, t44, t45;
248
249 t1 = 1/pow(3.0,0.666666666666667);
250 t2 = 1/pow(3.141592653589793,1.333333333333333);
251 t3 = pow(grada,2.0);
252 t4 = 1/pow(rhoa,2.666666666666667);
253 t5 = 2.5E-9*t1*t2*t3*t4+1.0;
254 t6 = 1/pow(t5,2.0);
255 t7 = 1/pow(3.141592653589793,8.0);
256 t8 = pow(grada,12.0);
257 t9 = 1/pow(rhoa,16.0);
258 t10 = 1/pow(3.0,0.333333333333333);
259 t11 = 1/pow(3.141592653589793,6.666666666666667);
260 t12 = pow(grada,10.0);
261 t13 = 1/pow(rhoa,13.33333333333333);
262 t14 = 1/pow(3.141592653589793,5.333333333333333);
263 t15 = pow(grada,8.0);
264 t16 = 1/pow(rhoa,10.66666666666667);
265 t17 = 1/pow(3.141592653589793,4.0);
266 t18 = pow(grada,6.0);
267 t19 = 1/pow(rhoa,8.0);
268 t20 = 1/pow(3.141592653589793,2.666666666666667);
269 t21 = pow(grada,4.0);
270 t22 = 1/pow(rhoa,5.333333333333333);
271 t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
272 t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
273 t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
274 t10*t11*t12*t13+1.0;
275 t24 = pow(t23,0.024974);
276 t25 = 1/pow(rhoa,2.333333333333333);
277 t26 = 1/t5;
278 t27 = pow(rhoa,0.333333333333333);
279 t28 = 1/pow(rhoa,17.0);
280 t29 = 1/pow(rhoa,14.33333333333333);
281 t30 = 1/pow(rhoa,11.66666666666667);
282 t31 = 1/pow(rhoa,9.0);
283 t32 = 1/pow(rhoa,6.333333333333333);
284 t33 = 1/pow(rhoa,3.666666666666667);
285 t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
286 t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
287 t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
288 t7*t8*t28;
289 t35 = 1/pow(t23,0.975026);
290 t36 = pow(rhoa,1.333333333333333);
291 t37 = 1/t36;
292 t38 = pow(grada,11.0);
293 t39 = pow(grada,9.0);
294 t40 = pow(grada,7.0);
295 t41 = pow(grada,5.0);
296 t42 = pow(grada,3.0);
297 t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
298 t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
299 t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
300 t7*t38*t9;
301 t44 = 1/pow(t5,3.0);
302 t45 = 1/pow(t23,1.975026);
303
304 /* code */
305 res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
306 t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
307 res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
308 t6*t24*t37;
309 res[2] = 0.024974*t26*t35*t36*(0.045864583333333*t7*t8/
310 pow(rhoa,18.0)+10.85299318415638*t10*t11*t12/pow(rhoa,15.33333333333333)+
311 0.65457561728395*t1*t14*t15/pow(rhoa,12.66666666666667)+2.802125*
312 t17*t18/pow(rhoa,10.0)+20.96333333333333*t10*t20*t21/pow(rhoa,
313 7.333333333333333)+12.07834022939412*t1*t2*t3/pow(rhoa,4.666666666666667))+
314 2.96296296296296E-17*t10*t20*t21*t24*t44/pow(rhoa,6.0)-6.666666666666668E-9*
315 t1*t2*t24*t3*t6/pow(rhoa,3.333333333333333)+0.444444444444444*
316 t24*t26/pow(rhoa,0.666666666666667)-0.024350299324*t26*pow(t34,
317 2.0)*t36*t45+0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*
318 t1*t2*t3*t34*t6*t35*t25;
319 res[3] = -2.222222222222222E-17*t10*t20*t24*t42*t44/pow(rhoa,
320 5.0)-1.2487E-10*t1*t2*grada*t34*t6*t35*t37-0.024350299324*
321 t34*t43*t26*t45*t36+0.024974*(-6.588185579669522*t1*t2*grada*
322 t33-13.24*t10*t20*t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*
323 t1*t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*
324 t7*t38*t28)*t26*t35*t36+0.033298666666667*t43*t26*t35*t27+
325 1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25+6.666666666666667E-9*
326 t1*t2*grada*t6*t24*t25;
327 res[4] = 1.666666666666667E-17*t10*t20*t24*t3*t44/pow(rhoa,
328 4.0)-0.024350299324*t26*t36*pow(t43,2.0)*t45-2.4974E-10*t1*
329 t2*grada*t43*t6*t35*t37-5.0E-9*t1*t2*t6*t24*t37+0.024974*(2.470569592376071*
330 t1*t2*t4+7.4475*t10*t20*t3*t22+1.167552083333333*t17*t21*t19+
331 0.294559027777778*t1*t14*t18*t16+5.111002604166666*t10*t11*
332 t15*t13+0.0222578125*t7*t12*t9)*t26*t35*t36;
333
334 }
335
336 static void
lg93_second(FunSecondFuncDrv * ds,real factor,const FunDensProp * dp)337 lg93_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
338 {
339 real res[5];
340
341 lg93_second_helper(dp->rhoa, dp->grada, res);
342
343 ds->df1000 += factor*res[0];
344 ds->df0010 += factor*res[1];
345
346 ds->df2000 += factor*res[2];
347 ds->df1010 += factor*res[3];
348 ds->df0020 += factor*res[4];
349
350
351 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
352 fabs(dp->grada-dp->gradb)>1e-13)
353 lg93_second_helper(dp->rhob, dp->gradb, res);
354 ds->df0100 += factor*res[0];
355 ds->df0001 += factor*res[1];
356
357 ds->df0200 += factor*res[2];
358 ds->df0101 += factor*res[3];
359 ds->df0002 += factor*res[4];
360
361 }
362
363 static void
lg93_third_helper(real rhoa,real grada,real * res)364 lg93_third_helper(real rhoa, real grada, real *res)
365 {
366 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
367 real t11, t12, t13, t14, t15, t16, t17, t18;
368 real t19, t20, t21, t22, t23, t24, t25, t26;
369 real t27, t28, t29, t30, t31, t32, t33, t34;
370 real t35, t36, t37, t38, t39, t40, t41, t42;
371 real t43, t44, t45, t46, t47, t48, t49, t50;
372 real t51, t52, t53, t54, t55, t56, t57, t58;
373 real t59, t60, t61, t62, t63;
374
375 t1 = 1/pow(3.0,0.666666666666667);
376 t2 = 1/pow(3.141592653589793,1.333333333333333);
377 t3 = pow(grada,2.0);
378 t4 = 1/pow(rhoa,2.666666666666667);
379 t5 = 2.5E-9*t1*t2*t3*t4+1.0;
380 t6 = 1/pow(t5,2.0);
381 t7 = 1/pow(3.141592653589793,8.0);
382 t8 = pow(grada,12.0);
383 t9 = 1/pow(rhoa,16.0);
384 t10 = 1/pow(3.0,0.333333333333333);
385 t11 = 1/pow(3.141592653589793,6.666666666666667);
386 t12 = pow(grada,10.0);
387 t13 = 1/pow(rhoa,13.33333333333333);
388 t14 = 1/pow(3.141592653589793,5.333333333333333);
389 t15 = pow(grada,8.0);
390 t16 = 1/pow(rhoa,10.66666666666667);
391 t17 = 1/pow(3.141592653589793,4.0);
392 t18 = pow(grada,6.0);
393 t19 = 1/pow(rhoa,8.0);
394 t20 = 1/pow(3.141592653589793,2.666666666666667);
395 t21 = pow(grada,4.0);
396 t22 = 1/pow(rhoa,5.333333333333333);
397 t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
398 t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
399 t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
400 t10*t11*t12*t13+1.0;
401 t24 = pow(t23,0.024974);
402 t25 = 1/pow(rhoa,2.333333333333333);
403 t26 = 1/t5;
404 t27 = pow(rhoa,0.333333333333333);
405 t28 = 1/pow(rhoa,17.0);
406 t29 = 1/pow(rhoa,14.33333333333333);
407 t30 = 1/pow(rhoa,11.66666666666667);
408 t31 = 1/pow(rhoa,9.0);
409 t32 = 1/pow(rhoa,6.333333333333333);
410 t33 = 1/pow(rhoa,3.666666666666667);
411 t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
412 t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
413 t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
414 t7*t8*t28;
415 t35 = 1/pow(t23,0.975026);
416 t36 = pow(rhoa,1.333333333333333);
417 t37 = 1/t36;
418 t38 = pow(grada,11.0);
419 t39 = pow(grada,9.0);
420 t40 = pow(grada,7.0);
421 t41 = pow(grada,5.0);
422 t42 = pow(grada,3.0);
423 t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
424 t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
425 t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
426 t7*t38*t9;
427 t44 = 1/pow(t5,3.0);
428 t45 = 1/pow(rhoa,6.0);
429 t46 = 1/pow(rhoa,3.333333333333333);
430 t47 = 1/pow(rhoa,0.666666666666667);
431 t48 = pow(t34,2.0);
432 t49 = 1/pow(t23,1.975026);
433 t50 = 1/pow(rhoa,18.0);
434 t51 = 1/pow(rhoa,15.33333333333333);
435 t52 = 1/pow(rhoa,12.66666666666667);
436 t53 = 1/pow(rhoa,10.0);
437 t54 = 1/pow(rhoa,7.333333333333333);
438 t55 = 1/pow(rhoa,4.666666666666667);
439 t56 = 12.07834022939412*t1*t2*t3*t55+20.96333333333333*
440 t10*t20*t21*t54+2.802125*t17*t18*t53+0.65457561728395*t1*t14*
441 t15*t52+10.85299318415638*t10*t11*t12*t51+0.045864583333333*
442 t7*t8*t50;
443 t57 = 1/pow(rhoa,5.0);
444 t58 = -6.588185579669522*t1*t2*grada*t33-13.24*t10*t20*
445 t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*t1*
446 t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*t7*
447 t38*t28;
448 t59 = 1/pow(rhoa,4.0);
449 t60 = pow(t43,2.0);
450 t61 = 2.470569592376071*t1*t2*t4+7.4475*t10*t20*t3*t22+
451 1.167552083333333*t17*t21*t19+0.294559027777778*t1*t14*t18*
452 t16+5.111002604166666*t10*t11*t15*t13+0.0222578125*t7*t12*
453 t9;
454 t62 = 1/pow(t5,4.0);
455 t63 = 1/pow(t23,2.975026);
456
457 /* code */
458 res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
459 t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
460 res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
461 t6*t24*t37;
462 res[2] = 0.444444444444444*t24*t26*t47-6.666666666666668E-9*
463 t1*t2*t3*t6*t24*t46+2.96296296296296E-17*t10*t20*t21*t44*t24*
464 t45-0.024350299324*t48*t26*t49*t36+0.024974*t56*t26*t35*t36+
465 0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*t1*
466 t2*t3*t34*t6*t35*t25;
467 res[3] = 0.024974*t58*t26*t35*t36-0.024350299324*t34*
468 t43*t26*t49*t36+0.033298666666667*t43*t26*t35*t27-1.2487E-10*
469 t1*t2*grada*t34*t6*t35*t37+6.666666666666667E-9*t1*t2*grada*
470 t6*t24*t25+1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25-2.222222222222222E-17*
471 t10*t20*t42*t44*t24*t57;
472 res[4] = 0.024974*t61*t26*t35*t36-0.024350299324*t60*
473 t26*t49*t36-5.0E-9*t1*t2*t6*t24*t37-2.4974E-10*t1*t2*grada*
474 t43*t6*t35*t37+1.666666666666667E-17*t10*t20*t3*t44*t24*t59;
475 res[5] = 0.024974*t26*t35*t36*(-0.8255625*t7*t8/pow(rhoa,
476 19.0)-166.4125621570645*t10*t11*t12/pow(rhoa,16.33333333333333)-
477 8.291291152263373*t1*t14*t15/pow(rhoa,13.66666666666667)-28.02125*
478 t17*t18/pow(rhoa,11.0)-153.7311111111111*t10*t20*t21/pow(rhoa,
479 8.333333333333334)-56.36558773717258*t1*t2*t3/pow(rhoa,5.666666666666667))+
480 1.975308641975309E-25*t17*t18*t24*t62/pow(rhoa,9.666666666666666)-
481 2.074074074074074E-16*t10*t20*t21*t24*t44/pow(rhoa,7.0)+2.518518518518519E-8*
482 t1*t2*t24*t3*t6/pow(rhoa,4.333333333333333)-0.296296296296296*
483 t24*t26/pow(rhoa,1.666666666666667)+0.048092474272682*t26*
484 pow(t34,3.0)*t36*t63+0.033298666666667*t34*t26*t35*t47-4.9948E-10*
485 t1*t2*t3*t34*t6*t35*t46+2.21991111111111E-18*t10*t20*t21*t34*
486 t44*t35*t45-0.073050897972*t56*t34*t26*t49*t36-0.097401197296*
487 t48*t26*t49*t27+0.099896*t56*t26*t35*t27-4.8700598648E-10*
488 t1*t2*t3*t48*t6*t49*t25+4.9948E-10*t1*t2*t3*t56*t6*t35*t25;
489 res[6] = -1.481481481481482E-25*t17*t24*t41*t62/pow(rhoa,
490 8.666666666666666)-1.109955555555555E-18*t10*t20*t42*t34*t44*
491 t35*t57+0.011099555555556*t43*t26*t35*t47-1.664933333333334E-10*
492 t1*t2*t3*t43*t6*t35*t46-1.555555555555556E-8*t1*t2*grada*t6*
493 t24*t46+7.3997037037037E-19*t10*t20*t21*t43*t44*t35*t45+1.407407407407407E-16*
494 t10*t20*t42*t44*t24*t45+1.2175149662E-10*t1*t2*grada*t48*t6*
495 t49*t37-1.2487E-10*t1*t2*grada*t56*t6*t35*t37+0.048092474272682*
496 t48*t43*t26*t63*t36-0.024350299324*t56*t43*t26*t49*t36-0.048700598648*
497 t58*t34*t26*t49*t36+0.024974*(24.15668045878825*t1*t2*grada*
498 t55+83.85333333333332*t10*t20*t42*t54+16.81275*t17*t41*t53+
499 5.236604938271604*t1*t14*t40*t52+108.5299318415638*t10*t11*
500 t39*t51+0.550375*t7*t38*t50)*t26*t35*t36-0.064934131530667*
501 t34*t43*t26*t49*t27+0.066597333333333*t58*t26*t35*t27-3.24670657653333E-10*
502 t1*t2*t3*t34*t43*t6*t49*t25+3.32986666666667E-10*t1*t2*t3*
503 t58*t6*t35*t25+3.32986666666667E-10*t1*t2*grada*t34*t6*t35*
504 t25;
505 res[7] = 1.111111111111111E-25*t17*t21*t24*t62/pow(rhoa,
506 7.666666666666667)+4.16233333333333E-19*t10*t20*t3*t34*t44*
507 t35*t59-1.109955555555555E-18*t10*t20*t42*t43*t44*t35*t57-
508 8.88888888888889E-17*t10*t20*t3*t44*t24*t57+2.4350299324E-10*
509 t1*t2*grada*t34*t43*t6*t49*t37-2.4974E-10*t1*t2*grada*t58*
510 t6*t35*t37-1.2487E-10*t1*t2*t34*t6*t35*t37+0.048092474272682*
511 t34*t60*t26*t63*t36-0.024350299324*t34*t61*t26*t49*t36-0.048700598648*
512 t58*t43*t26*t49*t36+0.024974*(-6.588185579669522*t1*t2*t33-
513 39.72*t10*t20*t3*t32-9.340416666666666*t17*t21*t31-3.141962962962963*
514 t1*t14*t18*t30-68.14670138888889*t10*t11*t15*t29-0.356125*
515 t7*t12*t28)*t26*t35*t36-0.032467065765333*t60*t26*t49*t27+
516 0.033298666666667*t61*t26*t35*t27-1.623353288266666E-10*t1*
517 t2*t3*t60*t6*t49*t25+1.664933333333333E-10*t1*t2*t3*t61*t6*
518 t35*t25+3.32986666666667E-10*t1*t2*grada*t43*t6*t35*t25+6.666666666666667E-9*
519 t1*t2*t6*t24*t25;
520 res[8] = -8.33333333333333E-26*t17*t24*t42*t62/pow(rhoa,
521 6.666666666666667)+0.048092474272682*t26*t36*pow(t43,3.0)*
522 t63+1.2487E-18*t10*t20*t3*t43*t44*t35*t59+5.0E-17*t10*t20*
523 grada*t44*t24*t59+3.6525448986E-10*t1*t2*grada*t60*t6*t49*
524 t37-3.7461E-10*t1*t2*grada*t61*t6*t35*t37-3.7461E-10*t1*t2*
525 t43*t6*t35*t37-0.073050897972*t61*t43*t26*t49*t36+0.024974*
526 (14.895*t10*t20*grada*t22+4.670208333333333*t17*t42*t19+1.767354166666667*
527 t1*t14*t41*t16+40.88802083333333*t10*t11*t40*t13+0.222578125*
528 t7*t39*t9)*t26*t35*t36;
529
530 }
531
532 static void
lg93_third(FunThirdFuncDrv * ds,real factor,const FunDensProp * dp)533 lg93_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
534 {
535 real res[9];
536
537 lg93_third_helper(dp->rhoa, dp->grada, res);
538
539 ds->df1000 += factor*res[0];
540 ds->df0010 += factor*res[1];
541
542 ds->df2000 += factor*res[2];
543 ds->df1010 += factor*res[3];
544 ds->df0020 += factor*res[4];
545
546 ds->df3000 += factor*res[5];
547 ds->df2010 += factor*res[6];
548 ds->df1020 += factor*res[7];
549 ds->df0030 += factor*res[8];
550
551
552 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
553 fabs(dp->grada-dp->gradb)>1e-13)
554 lg93_third_helper(dp->rhob, dp->gradb, res);
555
556 ds->df0100 += factor*res[0];
557 ds->df0001 += factor*res[1];
558
559 ds->df0200 += factor*res[2];
560 ds->df0101 += factor*res[3];
561 ds->df0002 += factor*res[4];
562
563 ds->df0300 += factor*res[5];
564 ds->df0201 += factor*res[6];
565 ds->df0102 += factor*res[7];
566 ds->df0003 += factor*res[8];
567
568 }
569
570 static void
lg93_fourth_helper(real rhoa,real grada,real * res)571 lg93_fourth_helper(real rhoa, real grada, real *res)
572 {
573 real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
574 real t11, t12, t13, t14, t15, t16, t17, t18;
575 real t19, t20, t21, t22, t23, t24, t25, t26;
576 real t27, t28, t29, t30, t31, t32, t33, t34;
577 real t35, t36, t37, t38, t39, t40, t41, t42;
578 real t43, t44, t45, t46, t47, t48, t49, t50;
579 real t51, t52, t53, t54, t55, t56, t57, t58;
580 real t59, t60, t61, t62, t63, t64, t65, t66;
581 real t67, t68, t69, t70, t71, t72, t73, t74;
582 real t75, t76, t77, t78, t79, t80, t81, t82;
583 real t83, t84, t85;
584
585 t1 = 1/pow(3.0,0.666666666666667);
586 t2 = 1/pow(3.141592653589793,1.333333333333333);
587 t3 = pow(grada,2.0);
588 t4 = 1/pow(rhoa,2.666666666666667);
589 t5 = 2.5E-9*t1*t2*t3*t4+1.0;
590 t6 = 1/pow(t5,2.0);
591 t7 = 1/pow(3.141592653589793,8.0);
592 t8 = pow(grada,12.0);
593 t9 = 1/pow(rhoa,16.0);
594 t10 = 1/pow(3.0,0.333333333333333);
595 t11 = 1/pow(3.141592653589793,6.666666666666667);
596 t12 = pow(grada,10.0);
597 t13 = 1/pow(rhoa,13.33333333333333);
598 t14 = 1/pow(3.141592653589793,5.333333333333333);
599 t15 = pow(grada,8.0);
600 t16 = 1/pow(rhoa,10.66666666666667);
601 t17 = 1/pow(3.141592653589793,4.0);
602 t18 = pow(grada,6.0);
603 t19 = 1/pow(rhoa,8.0);
604 t20 = 1/pow(3.141592653589793,2.666666666666667);
605 t21 = pow(grada,4.0);
606 t22 = 1/pow(rhoa,5.333333333333333);
607 t23 = 1.6861979166666666E-4*t7*t8*t9+1.235284796188036*
608 t1*t2*t3*t4+0.620625*t10*t20*t21*t22+0.038918402777778*t17*
609 t18*t19+0.005259982638889*t1*t14*t15*t16+0.056788917824074*
610 t10*t11*t12*t13+1.0;
611 t24 = pow(t23,0.024974);
612 t25 = 1/pow(rhoa,2.333333333333333);
613 t26 = 1/t5;
614 t27 = pow(rhoa,0.333333333333333);
615 t28 = 1/pow(rhoa,17.0);
616 t29 = 1/pow(rhoa,14.33333333333333);
617 t30 = 1/pow(rhoa,11.66666666666667);
618 t31 = 1/pow(rhoa,9.0);
619 t32 = 1/pow(rhoa,6.333333333333333);
620 t33 = 1/pow(rhoa,3.666666666666667);
621 t34 = -3.294092789834761*t1*t2*t3*t33-3.31*t10*t20*t21*
622 t32-0.311347222222222*t17*t18*t31-0.056106481481481*t1*t14*
623 t15*t30-0.757185570987654*t10*t11*t12*t29-0.002697916666667*
624 t7*t8*t28;
625 t35 = 1/pow(t23,0.975026);
626 t36 = pow(rhoa,1.333333333333333);
627 t37 = 1/t36;
628 t38 = pow(grada,11.0);
629 t39 = pow(grada,9.0);
630 t40 = pow(grada,7.0);
631 t41 = pow(grada,5.0);
632 t42 = pow(grada,3.0);
633 t43 = 2.470569592376071*t1*t2*grada*t4+2.4825*t10*t20*
634 t42*t22+0.233510416666667*t17*t41*t19+0.042079861111111*t1*
635 t14*t40*t16+0.567889178240741*t10*t11*t39*t13+0.0020234375*
636 t7*t38*t9;
637 t44 = 1/pow(t5,3.0);
638 t45 = 1/pow(rhoa,6.0);
639 t46 = 1/pow(rhoa,3.333333333333333);
640 t47 = 1/pow(rhoa,0.666666666666667);
641 t48 = pow(t34,2.0);
642 t49 = 1/pow(t23,1.975026);
643 t50 = 1/pow(rhoa,18.0);
644 t51 = 1/pow(rhoa,15.33333333333333);
645 t52 = 1/pow(rhoa,12.66666666666667);
646 t53 = 1/pow(rhoa,10.0);
647 t54 = 1/pow(rhoa,7.333333333333333);
648 t55 = 1/pow(rhoa,4.666666666666667);
649 t56 = 12.07834022939412*t1*t2*t3*t55+20.96333333333333*
650 t10*t20*t21*t54+2.802125*t17*t18*t53+0.65457561728395*t1*t14*
651 t15*t52+10.85299318415638*t10*t11*t12*t51+0.045864583333333*
652 t7*t8*t50;
653 t57 = 1/pow(rhoa,5.0);
654 t58 = -6.588185579669522*t1*t2*grada*t33-13.24*t10*t20*
655 t42*t32-1.868083333333333*t17*t41*t31-0.448851851851852*t1*
656 t14*t40*t30-7.571855709876543*t10*t11*t39*t29-0.032375*t7*
657 t38*t28;
658 t59 = 1/pow(rhoa,4.0);
659 t60 = pow(t43,2.0);
660 t61 = 2.470569592376071*t1*t2*t4+7.4475*t10*t20*t3*t22+
661 1.167552083333333*t17*t21*t19+0.294559027777778*t1*t14*t18*
662 t16+5.111002604166666*t10*t11*t15*t13+0.0222578125*t7*t12*
663 t9;
664 t62 = 1/pow(t5,4.0);
665 t63 = 1/pow(rhoa,9.666666666666666);
666 t64 = 1/pow(rhoa,7.0);
667 t65 = 1/pow(rhoa,4.333333333333333);
668 t66 = 1/pow(rhoa,1.666666666666667);
669 t67 = pow(t34,3.0);
670 t68 = 1/pow(t23,2.975026);
671 t69 = 1/pow(rhoa,19.0);
672 t70 = 1/pow(rhoa,16.33333333333333);
673 t71 = 1/pow(rhoa,13.66666666666667);
674 t72 = 1/pow(rhoa,11.0);
675 t73 = 1/pow(rhoa,8.333333333333334);
676 t74 = 1/pow(rhoa,5.666666666666667);
677 t75 = -56.36558773717258*t1*t2*t3*t74-153.7311111111111*
678 t10*t20*t21*t73-28.02125*t17*t18*t72-8.291291152263373*t1*
679 t14*t15*t71-166.4125621570645*t10*t11*t12*t70-0.8255625*t7*
680 t8*t69;
681 t76 = 1/pow(rhoa,8.666666666666666);
682 t77 = 24.15668045878825*t1*t2*grada*t55+83.85333333333332*
683 t10*t20*t42*t54+16.81275*t17*t41*t53+5.236604938271604*t1*
684 t14*t40*t52+108.5299318415638*t10*t11*t39*t51+0.550375*t7*
685 t38*t50;
686 t78 = 1/pow(rhoa,7.666666666666667);
687 t79 = -6.588185579669522*t1*t2*t33-39.72*t10*t20*t3*t32-
688 9.340416666666666*t17*t21*t31-3.141962962962963*t1*t14*t18*
689 t30-68.14670138888889*t10*t11*t15*t29-0.356125*t7*t12*t28;
690 t80 = 1/
691 pow(rhoa,6.666666666666667);
692 t81 = pow(t43,3.0);
693 t82 = 14.895*t10*t20*grada*t22+4.670208333333333*t17*
694 t42*t19+1.767354166666667*t1*t14*t41*t16+40.88802083333333*
695 t10*t11*t40*t13+0.222578125*t7*t39*t9;
696 t83 = 1/pow(t5,5.0);
697 t84 = 1/pow(t23,3.975026);
698 t85 = 1/pow(rhoa,9.333333333333334);
699
700 /* code */
701 res[0] = 0.024974*t34*t26*t35*t36+1.333333333333333*t24*
702 t26*t27+6.666666666666667E-9*t1*t2*t3*t6*t24*t25;
703 res[1] = 0.024974*t43*t26*t35*t36-5.0E-9*t1*t2*grada*
704 t6*t24*t37;
705 res[2] = 0.444444444444444*t24*t26*t47-6.666666666666668E-9*
706 t1*t2*t3*t6*t24*t46+2.96296296296296E-17*t10*t20*t21*t44*t24*
707 t45-0.024350299324*t48*t26*t49*t36+0.024974*t56*t26*t35*t36+
708 0.066597333333333*t34*t26*t35*t27+3.32986666666667E-10*t1*
709 t2*t3*t34*t6*t35*t25;
710 res[3] = 0.024974*t58*t26*t35*t36-0.024350299324*t34*
711 t43*t26*t49*t36+0.033298666666667*t43*t26*t35*t27-1.2487E-10*
712 t1*t2*grada*t34*t6*t35*t37+6.666666666666667E-9*t1*t2*grada*
713 t6*t24*t25+1.664933333333333E-10*t1*t2*t3*t43*t6*t35*t25-2.222222222222222E-17*
714 t10*t20*t42*t44*t24*t57;
715 res[4] = 0.024974*t61*t26*t35*t36-0.024350299324*t60*
716 t26*t49*t36-5.0E-9*t1*t2*t6*t24*t37-2.4974E-10*t1*t2*grada*
717 t43*t6*t35*t37+1.666666666666667E-17*t10*t20*t3*t44*t24*t59;
718 res[5] = -0.296296296296296*t24*t26*t66+2.518518518518519E-8*
719 t1*t2*t3*t6*t24*t65-2.074074074074074E-16*t10*t20*t21*t44*
720 t24*t64+1.975308641975309E-25*t17*t18*t62*t24*t63+0.033298666666667*
721 t34*t26*t35*t47-4.9948E-10*t1*t2*t3*t34*t6*t35*t46+2.21991111111111E-18*
722 t10*t20*t21*t34*t44*t35*t45+0.048092474272682*t67*t26*t68*
723 t36-0.073050897972*t56*t34*t26*t49*t36+0.024974*t75*t26*t35*
724 t36-0.097401197296*t48*t26*t49*t27+0.099896*t56*t26*t35*t27-
725 4.8700598648E-10*t1*t2*t3*t48*t6*t49*t25+4.9948E-10*t1*t2*
726 t3*t56*t6*t35*t25;
727 res[6] = 0.024974*t77*t26*t35*t36-0.024350299324*t56*
728 t43*t26*t49*t36-0.048700598648*t58*t34*t26*t49*t36+0.048092474272682*
729 t48*t43*t26*t68*t36+0.066597333333333*t58*t26*t35*t27-0.064934131530667*
730 t34*t43*t26*t49*t27+0.011099555555556*t43*t26*t35*t47-1.2487E-10*
731 t1*t2*grada*t56*t6*t35*t37+1.2175149662E-10*t1*t2*grada*t48*
732 t6*t49*t37+3.32986666666667E-10*t1*t2*grada*t34*t6*t35*t25+
733 3.32986666666667E-10*t1*t2*t3*t58*t6*t35*t25-3.24670657653333E-10*
734 t1*t2*t3*t34*t43*t6*t49*t25-1.555555555555556E-8*t1*t2*grada*
735 t6*t24*t46-1.664933333333334E-10*t1*t2*t3*t43*t6*t35*t46-1.109955555555555E-18*
736 t10*t20*t42*t34*t44*t35*t57+1.407407407407407E-16*t10*t20*
737 t42*t44*t24*t45+7.3997037037037E-19*t10*t20*t21*t43*t44*t35*
738 t45-1.481481481481482E-25*t17*t41*t62*t24*t76;
739 res[7] = 0.024974*t79*t26*t35*t36-0.048700598648*t58*
740 t43*t26*t49*t36-0.024350299324*t34*t61*t26*t49*t36+0.048092474272682*
741 t34*t60*t26*t68*t36+0.033298666666667*t61*t26*t35*t27-0.032467065765333*
742 t60*t26*t49*t27-1.2487E-10*t1*t2*t34*t6*t35*t37-2.4974E-10*
743 t1*t2*grada*t58*t6*t35*t37+2.4350299324E-10*t1*t2*grada*t34*
744 t43*t6*t49*t37+6.666666666666667E-9*t1*t2*t6*t24*t25+3.32986666666667E-10*
745 t1*t2*grada*t43*t6*t35*t25+1.664933333333333E-10*t1*t2*t3*
746 t61*t6*t35*t25-1.623353288266666E-10*t1*t2*t3*t60*t6*t49*t25+
747 4.16233333333333E-19*t10*t20*t3*t34*t44*t35*t59-8.88888888888889E-17*
748 t10*t20*t3*t44*t24*t57-1.109955555555555E-18*t10*t20*t42*t43*
749 t44*t35*t57+1.111111111111111E-25*t17*t21*t62*t24*t78;
750 res[8] = 0.024974*t82*t26*t35*t36-0.073050897972*t61*
751 t43*t26*t49*t36+0.048092474272682*t81*t26*t68*t36-3.7461E-10*
752 t1*t2*t43*t6*t35*t37-3.7461E-10*t1*t2*grada*t61*t6*t35*t37+
753 3.6525448986E-10*t1*t2*grada*t60*t6*t49*t37+5.0E-17*t10*t20*
754 grada*t44*t24*t59+1.2487E-18*t10*t20*t3*t43*t44*t35*t59-8.33333333333333E-26*
755 t17*t42*t62*t24*t80;
756 res[9] = 0.024974*t26*t35*t36*(15.6856875*t7*t8/pow(rhoa,
757 20.0)+2718.071848565387*t10*t11*t12/pow(rhoa,17.33333333333333)+
758 113.3143124142661*t1*t14*t15/pow(rhoa,14.66666666666667)+308.23375*
759 t17*t18/pow(rhoa,12.0)+1281.092592592593*t10*t20*t21*t85+319.4049971773113*
760 t1*t2*t3*t80)-0.143076361365561*t26*pow(t34,4.0)*t36*t84-0.029598814814815*
761 t34*t26*t35*t66+2.51589925925926E-9*t1*t2*t3*t34*t6*t35*t65-
762 2.071917037037037E-17*t10*t20*t21*t34*t44*t35*t64+1.973254320987654E-26*
763 t17*t18*t34*t62*t35*t63-0.073050897972*t26*t36*t49*pow(t56,
764 2.0)-0.064934131530667*t48*t26*t49*t47+0.066597333333333*t56*
765 t26*t35*t47+9.7401197296E-10*t1*t2*t3*t48*t6*t49*t46-9.9896E-10*
766 t1*t2*t3*t56*t6*t35*t46-4.32894210204444E-18*t10*t20*t21*t48*
767 t44*t49*t45+4.43982222222222E-18*t10*t20*t21*t56*t44*t35*t45+
768 0.493827160493827*t24*t26*t4+0.288554845636095*t56*t48*t26*
769 t68*t36-0.097401197296*t75*t34*t26*t49*t36+0.256493196120973*
770 t67*t26*t68*t27-0.389604789184*t56*t34*t26*t49*t27+0.133194666666667*
771 t75*t26*t35*t27+1.282465980604865E-9*t1*t2*t3*t67*t6*t68*t25-
772 1.94802394592E-9*t1*t2*t3*t56*t34*t6*t49*t25+6.65973333333333E-10*
773 t1*t2*t3*t75*t6*t35*t25-1.1111111111111112E-7*t1*t2*t3*t6*
774 t24*t22+1.563786008230453E-15*t10*t20*t21*t44*t24*t19-3.29218106995885E-24*
775 t17*t18*t62*t24*t16+5.26748971193416E-33*t1*t14*t15*t83*t24*
776 t13;
777 res[10] = -3.95061728395062E-33*t1*t14*t24*t40*t83/pow(rhoa,
778 12.33333333333333)-1.109955555555556E-26*t17*t41*t34*t62*t35*
779 t76-0.007399703703704*t43*t26*t35*t66+6.28974814814815E-10*
780 t1*t2*t3*t43*t6*t35*t65+5.185185185185187E-8*t1*t2*grada*t6*
781 t24*t65-5.17979259259259E-18*t10*t20*t21*t43*t44*t35*t64-9.1358024691358E-16*
782 t10*t20*t42*t44*t24*t64+4.93313580246914E-27*t17*t18*t43*t62*
783 t35*t63+2.22222222222222E-24*t17*t41*t62*t24*t63+1.623353288266667E-18*
784 t10*t20*t42*t48*t44*t49*t57-1.664933333333333E-18*t10*t20*
785 t42*t56*t44*t35*t57-0.032467065765333*t34*t43*t26*t49*t47+
786 0.033298666666667*t58*t26*t35*t47+4.8700598648E-10*t1*t2*t3*
787 t34*t43*t6*t49*t46-4.9948E-10*t1*t2*t3*t58*t6*t35*t46-1.165453333333334E-9*
788 t1*t2*grada*t34*t6*t35*t46-2.16447105102222E-18*t10*t20*t21*
789 t34*t43*t44*t49*t45+2.21991111111111E-18*t10*t20*t21*t58*t44*
790 t35*t45+1.054457777777778E-17*t10*t20*t42*t34*t44*t35*t45-
791 2.40462371363412E-10*t1*t2*grada*t67*t6*t68*t37+3.6525448986E-10*
792 t1*t2*grada*t56*t34*t6*t49*t37-1.2487E-10*t1*t2*grada*t75*
793 t6*t35*t37-0.143076361365561*t67*t43*t26*t84*t36+0.144277422818047*
794 t58*t48*t26*t68*t36+0.144277422818047*t56*t34*t43*t26*t68*
795 t36-0.073050897972*t56*t58*t26*t49*t36-0.024350299324*t75*
796 t43*t26*t49*t36-0.073050897972*t77*t34*t26*t49*t36+0.024974*
797 (-112.7311754743452*t1*t2*grada*t74-614.9244444444444*t10*
798 t20*t42*t73-168.1275*t17*t41*t72-66.33032921810698*t1*t14*
799 t40*t71-1664.125621570645*t10*t11*t39*t70-9.906749999999999*
800 t7*t38*t69)*t26*t35*t36+0.19236989709073*t48*t43*t26*t68*t27-
801 0.097401197296*t56*t43*t26*t49*t27-0.194802394592*t58*t34*
802 t26*t49*t27+0.099896*t77*t26*t35*t27+9.61849485453648E-10*
803 t1*t2*t3*t48*t43*t6*t68*t25-4.8700598648E-10*t1*t2*grada*t48*
804 t6*t49*t25-4.8700598648E-10*t1*t2*t3*t56*t43*t6*t49*t25-9.7401197296E-10*
805 t1*t2*t3*t58*t34*t6*t49*t25+4.9948E-10*t1*t2*t3*t77*t6*t35*
806 t25+4.9948E-10*t1*t2*grada*t56*t6*t35*t25;
807 res[11] = 2.962962962962963E-33*t1*t14*t18*t24*t83/pow(rhoa,
808 11.33333333333333)+5.54977777777778E-27*t17*t21*t34*t62*t35*
809 t78-7.3997037037037E-27*t17*t41*t43*t62*t35*t76-1.444444444444445E-24*
810 t17*t21*t62*t24*t76-4.05838322066667E-19*t10*t20*t3*t48*t44*
811 t49*t59+4.16233333333333E-19*t10*t20*t3*t56*t44*t35*t59-0.048700598648*
812 t26*t36*t49*pow(t58,2.0)+2.16447105102222E-18*t10*t20*t42*
813 t34*t43*t44*t49*t57-2.21991111111111E-18*t10*t20*t42*t58*t44*
814 t35*t57-4.43982222222222E-18*t10*t20*t3*t34*t44*t35*t57-0.010822355255111*
815 t60*t26*t49*t47+0.011099555555556*t61*t26*t35*t47+1.623353288266667E-10*
816 t1*t2*t3*t60*t6*t49*t46-1.664933333333334E-10*t1*t2*t3*t61*
817 t6*t35*t46-7.76968888888889E-10*t1*t2*grada*t43*t6*t35*t46-
818 1.555555555555556E-8*t1*t2*t6*t24*t46-7.21490350340741E-19*
819 t10*t20*t21*t60*t44*t49*t45+7.3997037037037E-19*t10*t20*t21*
820 t61*t44*t35*t45+7.02971851851852E-18*t10*t20*t42*t43*t44*t35*
821 t45+4.74074074074074E-16*t10*t20*t3*t44*t24*t45-4.80924742726824E-10*
822 t1*t2*grada*t48*t43*t6*t68*t37+1.2175149662E-10*t1*t2*t48*
823 t6*t49*t37+2.4350299324E-10*t1*t2*grada*t56*t43*t6*t49*t37+
824 4.8700598648E-10*t1*t2*grada*t58*t34*t6*t49*t37-2.4974E-10*
825 t1*t2*grada*t77*t6*t35*t37-1.2487E-10*t1*t2*t56*t6*t35*t37-
826 0.143076361365561*t48*t60*t26*t84*t36+0.048092474272682*t48*
827 t61*t26*t68*t36+0.048092474272682*t56*t60*t26*t68*t36+0.19236989709073*
828 t58*t34*t43*t26*t68*t36-0.024350299324*t56*t61*t26*t49*t36-
829 0.048700598648*t77*t43*t26*t49*t36-0.048700598648*t79*t34*
830 t26*t49*t36+0.024974*(24.15668045878825*t1*t2*t55+251.56*t10*
831 t20*t3*t54+84.06375*t17*t21*t53+36.65623456790123*t1*t14*t18*
832 t52+976.769386574074*t10*t11*t15*t51+6.054125*t7*t12*t50)*
833 t26*t35*t36+0.128246598060486*t34*t60*t26*t68*t27-0.064934131530667*
834 t34*t61*t26*t49*t27-0.129868263061333*t58*t43*t26*t49*t27+
835 0.066597333333333*t79*t26*t35*t27+6.41232990302432E-10*t1*
836 t2*t3*t34*t60*t6*t68*t25-3.24670657653333E-10*t1*t2*t3*t34*
837 t61*t6*t49*t25-6.49341315306667E-10*t1*t2*t3*t58*t43*t6*t49*
838 t25-6.49341315306667E-10*t1*t2*grada*t34*t43*t6*t49*t25+3.32986666666667E-10*
839 t1*t2*t3*t79*t6*t35*t25+6.65973333333333E-10*t1*t2*grada*t58*
840 t6*t35*t25+3.32986666666667E-10*t1*t2*t34*t6*t35*t25;
841 res[12] = -2.222222222222222E-33*t1*t14*t24*t41*t83/pow(rhoa,
842 10.33333333333333)-2.08116666666667E-27*t17*t42*t34*t62*t35*
843 t80+8.32466666666667E-27*t17*t21*t43*t62*t35*t78+8.88888888888889E-25*
844 t17*t42*t62*t24*t78-1.2175149662E-18*t10*t20*t3*t34*t43*t44*
845 t49*t59+1.2487E-18*t10*t20*t3*t58*t44*t35*t59+1.2487E-18*t10*
846 t20*grada*t34*t44*t35*t59+1.623353288266667E-18*t10*t20*t42*
847 t60*t44*t49*t57-1.664933333333333E-18*t10*t20*t42*t61*t44*
848 t35*t57-6.65973333333333E-18*t10*t20*t3*t43*t44*t35*t57-2.0E-16*
849 t10*t20*grada*t44*t24*t57-7.21387114090236E-10*t1*t2*grada*
850 t34*t60*t6*t68*t37+3.6525448986E-10*t1*t2*grada*t34*t61*t6*
851 t49*t37+7.3050897972E-10*t1*t2*grada*t58*t43*t6*t49*t37+3.6525448986E-10*
852 t1*t2*t34*t43*t6*t49*t37-3.7461E-10*t1*t2*grada*t79*t6*t35*
853 t37-3.7461E-10*t1*t2*t58*t6*t35*t37-0.143076361365561*t34*
854 t81*t26*t84*t36+0.144277422818047*t58*t60*t26*t68*t36+0.144277422818047*
855 t34*t61*t43*t26*t68*t36-0.073050897972*t58*t61*t26*t49*t36-
856 0.073050897972*t79*t43*t26*t49*t36-0.024350299324*t82*t34*
857 t26*t49*t36+0.024974*(-79.44*t10*t20*grada*t32-37.36166666666666*
858 t17*t42*t31-18.85177777777778*t1*t14*t41*t30-545.1736111111111*
859 t10*t11*t40*t29-3.56125*t7*t39*t28)*t26*t35*t36+0.064123299030243*
860 t81*t26*t68*t27-0.097401197296*t61*t43*t26*t49*t27+0.033298666666667*
861 t82*t26*t35*t27+3.20616495151216E-10*t1*t2*t3*t81*t6*t68*t25-
862 4.8700598648E-10*t1*t2*grada*t60*t6*t49*t25-4.8700598648E-10*
863 t1*t2*t3*t61*t43*t6*t49*t25+1.664933333333333E-10*t1*t2*t3*
864 t82*t6*t35*t25+4.9948E-10*t1*t2*grada*t61*t6*t35*t25+4.9948E-10*
865 t1*t2*t43*t6*t35*t25;
866 res[13] = 1.666666666666667E-33*t1*t14*t21*t83*t24*t85-
867 0.143076361365561*t26*t36*pow(t43,4.0)*t84-8.32466666666667E-27*
868 t17*t42*t43*t62*t35*t80-5.0E-25*t17*t3*t62*t24*t80-0.073050897972*
869 t26*t36*t49*pow(t61,2.0)-2.4350299324E-18*t10*t20*t3*t60*t44*
870 t49*t59+2.4974E-18*t10*t20*t3*t61*t44*t35*t59+4.9948E-18*t10*
871 t20*grada*t43*t44*t35*t59+5.0E-17*t10*t20*t44*t24*t59-9.61849485453648E-10*
872 t1*t2*grada*t81*t6*t68*t37+7.3050897972E-10*t1*t2*t60*t6*t49*
873 t37+1.46101795944E-9*t1*t2*grada*t61*t43*t6*t49*t37-4.9948E-10*
874 t1*t2*grada*t82*t6*t35*t37-7.4922E-10*t1*t2*t61*t6*t35*t37+
875 0.288554845636095*t61*t60*t26*t68*t36-0.097401197296*t82*t43*
876 t26*t49*t36+0.024974*(14.895*t10*t20*t22+14.010625*t17*t3*
877 t19+8.836770833333333*t1*t14*t21*t16+286.2161458333333*t10*
878 t11*t18*t13+2.003203125*t7*t15*t9)*t26*t35*t36;
879
880 }
881
882 static void
lg93_fourth(FunFourthFuncDrv * ds,real factor,const FunDensProp * dp)883 lg93_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
884 {
885 real res[14];
886
887 lg93_fourth_helper(dp->rhoa, dp->grada, res);
888
889 ds->df1000 += factor*res[0];
890 ds->df0010 += factor*res[1];
891
892 ds->df2000 += factor*res[2];
893 ds->df1010 += factor*res[3];
894 ds->df0020 += factor*res[4];
895
896 ds->df3000 += factor*res[5];
897 ds->df2010 += factor*res[6];
898 ds->df1020 += factor*res[7];
899 ds->df0030 += factor*res[8];
900
901 ds->df4000 += factor*res[9];
902 ds->df3010 += factor*res[10];
903 ds->df2020 += factor*res[11];
904 ds->df1030 += factor*res[12];
905 ds->df0040 += factor*res[13];
906
907
908 if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
909 fabs(dp->grada-dp->gradb)>1e-13)
910 lg93_fourth_helper(dp->rhob, dp->gradb, res);
911
912 ds->df0100 += factor*res[0];
913 ds->df0001 += factor*res[1];
914
915 ds->df0200 += factor*res[2];
916 ds->df0101 += factor*res[3];
917 ds->df0002 += factor*res[4];
918
919 ds->df0300 += factor*res[5];
920 ds->df0201 += factor*res[6];
921 ds->df0102 += factor*res[7];
922 ds->df0003 += factor*res[8];
923
924 ds->df0400 += factor*res[9];
925 ds->df0301 += factor*res[10];
926 ds->df0202 += factor*res[11];
927 ds->df0103 += factor*res[12];
928 ds->df0004 += factor*res[13];
929
930 }
931