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