56 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
57 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
58 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
59 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
60 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
61 lza <= tp->
lmax[0] - lxa - lya; lza++) {
62 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
63 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
64 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
65 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
66 const double pab =
idx2(tp->
pab[0], jco, ico);
71 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
72 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
73 idx2(tp->
pab_prep[0], jco_l, ico_l) += 0.5 * lxa * lxb * pab;
75 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
76 jco_l = coset((lxb + 1), lyb, lzb);
79 ico_l = coset((lxa + 1), lya, lza);
80 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
83 ico_l = coset((lxa + 1), lya, lza);
84 jco_l = coset((lxb + 1), lyb, lzb);
86 2.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
90 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
91 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
92 idx2(tp->
pab_prep[0], jco_l, ico_l) += 0.5 * lya * lyb * pab;
94 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
95 jco_l = coset(lxb, (lyb + 1), lzb);
98 ico_l = coset(lxa, (lya + 1), lza);
99 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
102 ico_l = coset(lxa, (lya + 1), lza);
103 jco_l = coset(lxb, (lyb + 1), lzb);
105 2.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
109 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
110 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
111 idx2(tp->
pab_prep[0], jco_l, ico_l) += 0.5 * lza * lzb * pab;
113 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
114 jco_l = coset(lxb, lyb, (lzb + 1));
117 ico_l = coset(lxa, lya, (lza + 1));
118 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
121 ico_l = coset(lxa, lya, (lza + 1));
122 jco_l = coset(lxb, lyb, (lzb + 1));
124 2.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
142 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
143 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
144 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
145 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
146 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
147 lza <= tp->
lmax[0] - lxa - lya; lza++) {
148 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
149 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
150 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
151 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
152 const double pab =
idx2(tp->
pab[0], jco, ico);
158 ico_l = coset(lxa, lya, lza);
159 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
163 jco_l = coset((lxb + 1), lyb, lzb);
166 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
167 jco_l = coset(lxb, lyb, lzb);
170 ico_l = coset(lxa + 1, lya, lza);
175 ico_l = coset(lxa, lya, lza);
176 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
180 jco_l = coset(lxb, lyb + 1, lzb);
183 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
184 jco_l = coset(lxb, lyb, lzb);
187 ico_l = coset(lxa, lya + 1, lza);
192 ico_l = coset(lxa, lya, lza);
193 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
197 jco_l = coset(lxb, lyb, lzb + 1);
200 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
201 jco_l = coset(lxb, lyb, lzb);
204 ico_l = coset(lxa, lya, lza + 1);
230 assert(1 <= tp->
dir1 && tp->
dir1 <= 3);
231 assert(1 <= tp->
dir2 && tp->
dir2 <= 3);
233 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
234 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
235 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
236 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
237 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
238 lza <= tp->
lmax[0] - lxa - lya; lza++) {
239 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
240 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
241 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
242 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
243 const double pab =
idx2(tp->
pab[0], jco, ico);
252 ico_l = coset(lxa, lya, lza);
253 jco_l = coset(lxb, lyb, lzb);
256 ico_l = coset(lxa, lya, lza);
257 jco_l = coset((lxb + 2), lyb, lzb);
259 2.0 * tp->
zeta[1] * pab;
261 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
262 jco_l = coset((lxb + 1), lyb, lzb);
265 ico_l = coset((lxa + 1), lya, lza);
266 jco_l = coset((lxb + 1), lyb, lzb);
268 2.0 * tp->
zeta[0] * pab;
271 ico_l = coset(lxa, lya, lza);
272 jco_l = coset(
imax(lxb - 1, 0), (lyb + 1), lzb);
275 ico_l = coset(lxa, lya, lza);
276 jco_l = coset((lxb + 1), (lyb + 1), lzb);
278 2.0 * tp->
zeta[1] * pab;
280 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
281 jco_l = coset(lxb, (lyb + 1), lzb);
284 ico_l = coset((lxa + 1), lya, lza);
285 jco_l = coset(lxb, (lyb + 1), lzb);
287 2.0 * tp->
zeta[0] * pab;
290 ico_l = coset(lxa, lya, lza);
291 jco_l = coset(
imax(lxb - 1, 0), lyb, (lzb + 1));
294 ico_l = coset(lxa, lya, lza);
295 jco_l = coset((lxb + 1), lyb, (lzb + 1));
297 2.0 * tp->
zeta[1] * pab;
299 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
300 jco_l = coset(lxb, lyb, (lzb + 1));
303 ico_l = coset((lxa + 1), lya, lza);
304 jco_l = coset(lxb, lyb, (lzb + 1));
306 2.0 * tp->
zeta[0] * pab;
315 ico_l = coset(lxa, lya, lza);
316 jco_l = coset((lxb + 1),
imax(lyb - 1, 0), lzb);
319 ico_l = coset(lxa, lya, lza);
320 jco_l = coset((lxb + 1), (lyb + 1), lzb);
322 2.0 * tp->
zeta[1] * pab;
324 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
325 jco_l = coset((lxb + 1), lyb, lzb);
328 ico_l = coset(lxa, (lya + 1), lza);
329 jco_l = coset((lxb + 1), lyb, lzb);
331 2.0 * tp->
zeta[0] * pab;
334 ico_l = coset(lxa, lya, lza);
335 jco_l = coset(lxb, lyb, lzb);
338 ico_l = coset(lxa, lya, lza);
339 jco_l = coset(lxb, (lyb + 2), lzb);
341 2.0 * tp->
zeta[1] * pab;
343 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
344 jco_l = coset(lxb, (lyb + 1), lzb);
347 ico_l = coset(lxa, (lya + 1), lza);
348 jco_l = coset(lxb, (lyb + 1), lzb);
350 2.0 * tp->
zeta[0] * pab;
353 ico_l = coset(lxa, lya, lza);
354 jco_l = coset(lxb,
imax(lyb - 1, 0), (lzb + 1));
357 ico_l = coset(lxa, lya, lza);
358 jco_l = coset(lxb, (lyb + 1), (lzb + 1));
360 2.0 * tp->
zeta[1] * pab;
362 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
363 jco_l = coset(lxb, lyb, (lzb + 1));
366 ico_l = coset(lxa, (lya + 1), lza);
367 jco_l = coset(lxb, lyb, (lzb + 1));
369 2.0 * tp->
zeta[0] * pab;
378 ico_l = coset(lxa, lya, lza);
379 jco_l = coset((lxb + 1), lyb,
imax(lzb - 1, 0));
382 ico_l = coset(lxa, lya, lza);
383 jco_l = coset((lxb + 1), lyb, (lzb + 1));
385 2.0 * tp->
zeta[1] * pab;
387 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
388 jco_l = coset((lxb + 1), lyb, lzb);
391 ico_l = coset(lxa, lya, (lza + 1));
392 jco_l = coset((lxb + 1), lyb, lzb);
394 2.0 * tp->
zeta[0] * pab;
397 ico_l = coset(lxa, lya, lza);
398 jco_l = coset(lxb, (lyb + 1),
imax(lzb - 1, 0));
401 ico_l = coset(lxa, lya, lza);
402 jco_l = coset(lxb, (lyb + 1), (lzb + 1));
404 -2.0 * tp->
zeta[1] * pab;
406 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
407 jco_l = coset(lxb, (lyb + 1), lzb);
410 ico_l = coset(lxa, lya, (lza + 1));
411 jco_l = coset(lxb, (lyb + 1), lzb);
413 2.0 * tp->
zeta[0] * pab;
416 ico_l = coset(lxa, lya, lza);
417 jco_l = coset(lxb, lyb, lzb);
420 ico_l = coset(lxa, lya, lza);
421 jco_l = coset(lxb, lyb, (lzb + 2));
423 2.0 * tp->
zeta[1] * pab;
425 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
426 jco_l = coset(lxb, lyb, (lzb + 1));
429 ico_l = coset(lxa, lya, (lza + 1));
430 jco_l = coset(lxb, lyb, (lzb + 1));
432 2.0 * tp->
zeta[0] * pab;
456 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
457 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
458 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
459 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
460 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
461 lza <= tp->
lmax[0] - lxa - lya; lza++) {
462 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
463 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
464 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
465 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
466 const double pab =
idx2(tp->
pab[0], jco, ico);
474 ico_l = coset(lxa, lya, lza);
475 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
478 ico_l = coset(lxa, lya, lza);
479 jco_l = coset((lxb + 1), lyb, lzb);
482 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
483 jco_l = coset(lxb, lyb, lzb);
486 ico_l = coset((lxa + 1), lya, lza);
487 jco_l = coset(lxb, lyb, lzb);
491 ico_l = coset(lxa, lya, lza);
492 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
495 ico_l = coset(lxa, lya, lza);
496 jco_l = coset(lxb, (lyb + 1), lzb);
499 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
500 jco_l = coset(lxb, lyb, lzb);
503 ico_l = coset(lxa, (lya + 1), lza);
504 jco_l = coset(lxb, lyb, lzb);
508 ico_l = coset(lxa, lya, lza);
509 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
512 ico_l = coset(lxa, lya, lza);
513 jco_l = coset(lxb, lyb, (lzb + 1));
516 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
517 jco_l = coset(lxb, lyb, lzb);
520 ico_l = coset(lxa, lya, (lza + 1));
521 jco_l = coset(lxb, lyb, lzb);
545 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
546 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
547 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
548 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
549 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
550 lza <= tp->
lmax[0] - lxa - lya; lza++) {
551 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
552 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
553 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
554 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
555 const double pab =
idx2(tp->
pab[0], jco, ico);
564 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
565 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
568 ico_l = coset(
imax(lxa - 1, 0), lya, lza);
569 jco_l = coset((lxb + 1), lyb, lzb);
571 2.0 * lxa * tp->
zeta[1] * pab;
573 ico_l = coset((lxa + 1), lya, lza);
574 jco_l = coset(
imax(lxb - 1, 0), lyb, lzb);
576 2.0 * tp->
zeta[0] * lxb * pab;
578 ico_l = coset((lxa + 1), lya, lza);
579 jco_l = coset((lxb + 1), lyb, lzb);
581 4.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
585 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
586 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
589 ico_l = coset(lxa,
imax(lya - 1, 0), lza);
590 jco_l = coset(lxb, (lyb + 1), lzb);
592 2.0 * lya * tp->
zeta[1] * pab;
594 ico_l = coset(lxa, (lya + 1), lza);
595 jco_l = coset(lxb,
imax(lyb - 1, 0), lzb);
597 2.0 * tp->
zeta[0] * lyb * pab;
599 ico_l = coset(lxa, (lya + 1), lza);
600 jco_l = coset(lxb, (lyb + 1), lzb);
602 4.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
606 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
607 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
610 ico_l = coset(lxa, lya,
imax(lza - 1, 0));
611 jco_l = coset(lxb, lyb, (lzb + 1));
613 2.0 * lza * tp->
zeta[1] * pab;
615 ico_l = coset(lxa, lya, (lza + 1));
616 jco_l = coset(lxb, lyb,
imax(lzb - 1, 0));
618 2.0 * tp->
zeta[0] * lzb * pab;
620 ico_l = coset(lxa, lya, (lza + 1));
621 jco_l = coset(lxb, lyb, (lzb + 1));
623 4.0 * tp->
zeta[0] * tp->
zeta[1] * pab;
638 const int lx,
const int ly,
const int lz,
639 const double zet,
tensor *
const pab_prep) {
646 jco_l = coset(
imax(lx - 1, 0),
imax(ly - 1, 0), lz);
647 idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
649 jco_l = coset(lx + 1,
imax(ly - 1, 0), lz);
650 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
652 jco_l = coset(
imax(lx - 1, 0), ly + 1, lz);
653 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
655 jco_l = coset(lx + 1, ly + 1, lz);
656 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
661 jco_l = coset(lx,
imax(ly - 1, 0),
imax(lz - 1, 0));
662 idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
664 jco_l = coset(lx, ly + 1,
imax(lz - 1, 0));
665 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
667 jco_l = coset(lx,
imax(ly - 1, 0), lz + 1);
668 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
670 jco_l = coset(lx, ly + 1, lz + 1);
671 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
676 jco_l = coset(
imax(lx - 1, 0), ly,
imax(lz - 1, 0));
677 idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
679 jco_l = coset(
imax(lx - 1, 0), ly, lz + 1);
680 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
682 jco_l = coset(lx + 1, ly,
imax(lz - 1, 0));
683 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
685 jco_l = coset(lx + 1, ly, lz + 1);
686 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
698 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
699 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
700 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
701 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
702 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
703 lza <= tp->
lmax[0] - lxa - lya; lza++) {
704 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
705 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
706 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
707 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
708 const double pab =
idx2(tp->
pab[0], jco, ico);
715 if ((tp->
dir1 ==
'X' && tp->
dir2 ==
'Y') ||
716 (tp->
dir1 ==
'Y' && tp->
dir2 ==
'X')) {
718 ico_l = coset(
imax(lxa - 1, 0),
imax(lya - 1, 0), lza);
719 func_a = lxa * lya * pab;
723 ico_l = coset(lxa + 1,
imax(lya - 1, 0), lza);
724 func_a = -2.0 * tp->
zeta[0] * lya * pab;
728 ico_l = coset(
imax(lxa - 1, 0), lya + 1, lza);
729 func_a = -2.0 * tp->
zeta[0] * lxa * pab;
733 ico_l = coset(lxa + 1, lya + 1, lza);
734 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
737 }
else if ((tp->
dir1 ==
'Y' && tp->
dir2 ==
'Z') ||
738 (tp->
dir1 ==
'Z' && tp->
dir2 ==
'Y')) {
740 ico_l = coset(lxa,
imax(lya - 1, 0),
imax(lza - 1, 0));
741 func_a = lya * lza * pab;
745 ico_l = coset(lxa, lya + 1,
imax(lza - 1, 0));
746 func_a = -2.0 * tp->
zeta[0] * lza * pab;
750 ico_l = coset(lxa,
imax(lya - 1, 0), lza + 1);
751 func_a = -2.0 * tp->
zeta[0] * lya * pab;
755 ico_l = coset(lxa, lya + 1, lza + 1);
756 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
759 }
else if ((tp->
dir1 ==
'Z' && tp->
dir2 ==
'X') ||
760 (tp->
dir1 ==
'X' && tp->
dir2 ==
'Z')) {
762 ico_l = coset(
imax(lxa - 1, 0), lya,
imax(lza - 1, 0));
763 func_a = lza * lxa * pab;
767 ico_l = coset(
imax(lxa - 1, 0), lya, lza + 1);
768 func_a = -2.0 * tp->
zeta[0] * lxa * pab;
772 ico_l = coset(lxa + 1, lya,
imax(lza - 1, 0));
773 func_a = -2.0 * tp->
zeta[0] * lza * pab;
777 ico_l = coset(lxa + 1, lya, lza + 1);
778 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
792 const int lx,
const int ly,
const int lz,
793 const double zet,
tensor *
const pab_prep) {
799 jco_l = coset(
imax(lx - 2, 0), ly, lz);
800 idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
802 jco_l = coset(lx, ly, lz);
803 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
805 jco_l = coset(lx + 2, ly, lz);
806 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
810 jco_l = coset(lx,
imax(ly - 2, 0), lz);
811 idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
813 jco_l = coset(lx, ly, lz);
814 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
816 jco_l = coset(lx, ly + 2, lz);
817 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
821 jco_l = coset(lx, ly,
imax(lz - 2, 0));
822 idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
824 jco_l = coset(lx, ly, lz);
825 idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
827 jco_l = coset(lx, ly, lz + 2);
828 idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
831 printf(
"Wrong value for ider: should be 1, 2, or 3\n");
843 for (
int lxa = 0; lxa <= tp->
lmax[0]; lxa++) {
844 for (
int lxb = 0; lxb <= tp->
lmax[1]; lxb++) {
845 for (
int lya = 0; lya <= tp->
lmax[0] - lxa; lya++) {
846 for (
int lyb = 0; lyb <= tp->
lmax[1] - lxb; lyb++) {
847 for (
int lza =
imax(tp->
lmin[0] - lxa - lya, 0);
848 lza <= tp->
lmax[0] - lxa - lya; lza++) {
849 for (
int lzb =
imax(tp->
lmin[1] - lxb - lyb, 0);
850 lzb <= tp->
lmax[1] - lxb - lyb; lzb++) {
851 const int ico = tp->
offset[0] + coset(lxa, lya, lza);
852 const int jco = tp->
offset[1] + coset(lxb, lyb, lzb);
853 const double pab =
idx2(tp->
pab[0], jco, ico);
862 ico_l = coset(
imax(lxa - 2, 0), lya, lza);
863 func_a = lxa * (lxa - 1) * pab;
867 ico_l = coset(lxa, lya, lza);
868 func_a = -2.0 * tp->
zeta[0] * (2 * lxa + 1) * pab;
872 ico_l = coset(lxa + 2, lya, lza);
873 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
879 ico_l = coset(lxa,
imax(lya - 2, 0), lza);
880 func_a = lya * (lya - 1) * pab;
884 ico_l = coset(lxa, lya, lza);
885 func_a = -2.0 * tp->
zeta[0] * (2 * lya + 1) * pab;
889 ico_l = coset(lxa, lya + 2, lza);
890 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
896 ico_l = coset(lxa, lya,
imax(lza - 2, 0));
897 func_a = lza * (lza - 1) * pab;
901 ico_l = coset(lxa, lya, lza);
902 func_a = -2.0 * tp->
zeta[0] * (2 * lza + 1) * pab;
906 ico_l = coset(lxa, lya, lza + 2);
907 func_a = 4.0 * tp->
zeta[0] * tp->
zeta[0] * pab;
912 printf(
"Wrong value for ider: should be 'X', 'Y' or 'Z'.\n");