1!
2!   compute the values of all (x-xp)**lp*exp(..)
3!
4! still requires the old trick:
5!  new trick to avoid to many exps (reuse the result from the previous gridpoint):
6!  exp( -a*(x+d)**2)=exp(-a*x**2)*(-2*a*x*d)*exp(-a*d**2)
7!  exp(-2*a*(x+d)*d)=exp(-2*a*x*d)*exp(-2*a*d**2)
8
9      iaxis = 3
10      t_exp_1 = EXP(-zetp*dr(iaxis)**2)
11      t_exp_2 = t_exp_1**2
12      t_exp_min_1 = EXP(-zetp*(+dr(iaxis) - roffset(iaxis))**2)
13      t_exp_min_2 = EXP(-2*zetp*(+dr(iaxis) - roffset(iaxis))*(-dr(iaxis)))
14      t_exp_plus_1 = EXP(-zetp*(-roffset(iaxis))**2)
15      t_exp_plus_2 = EXP(-2*zetp*(-roffset(iaxis))*(+dr(iaxis)))
16      DO ig = 0, lb_cube(iaxis), -1
17         rpg = REAL(ig, dp)*dr(iaxis) - roffset(iaxis)
18         t_exp_min_1 = t_exp_min_1*t_exp_min_2*t_exp_1
19         t_exp_min_2 = t_exp_min_2*t_exp_2
20         pg = t_exp_min_1
21         ! pg  = EXP(-zetp*rpg**2)
22         DO icoef = 0, lp
23            pol_z(1, icoef, ig) = pg
24            pg = pg*(rpg)
25         ENDDO
26
27         rpg = REAL(1 - ig, dp)*dr(iaxis) - roffset(iaxis)
28         t_exp_plus_1 = t_exp_plus_1*t_exp_plus_2*t_exp_1
29         t_exp_plus_2 = t_exp_plus_2*t_exp_2
30         pg = t_exp_plus_1
31         ! pg  = EXP(-zetp*rpg**2)
32         DO icoef = 0, lp
33            pol_z(2, icoef, ig) = pg
34            pg = pg*(rpg)
35         ENDDO
36      ENDDO
37
38      iaxis = 2
39      t_exp_1 = EXP(-zetp*dr(iaxis)**2)
40      t_exp_2 = t_exp_1**2
41      t_exp_min_1 = EXP(-zetp*(+dr(iaxis) - roffset(iaxis))**2)
42      t_exp_min_2 = EXP(-2*zetp*(+dr(iaxis) - roffset(iaxis))*(-dr(iaxis)))
43      t_exp_plus_1 = EXP(-zetp*(-roffset(iaxis))**2)
44      t_exp_plus_2 = EXP(-2*zetp*(-roffset(iaxis))*(+dr(iaxis)))
45      DO ig = 0, lb_cube(iaxis), -1
46         rpg = REAL(ig, dp)*dr(iaxis) - roffset(iaxis)
47         t_exp_min_1 = t_exp_min_1*t_exp_min_2*t_exp_1
48         t_exp_min_2 = t_exp_min_2*t_exp_2
49         pg = t_exp_min_1
50         ! pg  = EXP(-zetp*rpg**2)
51         DO icoef = 0, lp
52            pol_y(1, icoef, ig) = pg
53            pg = pg*(rpg)
54         ENDDO
55
56         rpg = REAL(1 - ig, dp)*dr(iaxis) - roffset(iaxis)
57         t_exp_plus_1 = t_exp_plus_1*t_exp_plus_2*t_exp_1
58         t_exp_plus_2 = t_exp_plus_2*t_exp_2
59         pg = t_exp_plus_1
60         ! pg  = EXP(-zetp*rpg**2)
61         DO icoef = 0, lp
62            pol_y(2, icoef, ig) = pg
63            pg = pg*(rpg)
64         ENDDO
65      ENDDO
66
67      iaxis = 1
68      t_exp_1 = EXP(-zetp*dr(iaxis)**2)
69      t_exp_2 = t_exp_1**2
70      t_exp_min_1 = EXP(-zetp*(+dr(iaxis) - roffset(iaxis))**2)
71      t_exp_min_2 = EXP(-2*zetp*(+dr(iaxis) - roffset(iaxis))*(-dr(iaxis)))
72      t_exp_plus_1 = EXP(-zetp*(-roffset(iaxis))**2)
73      t_exp_plus_2 = EXP(-2*zetp*(-roffset(iaxis))*(+dr(iaxis)))
74      DO ig = 0, lb_cube(1), -1
75
76         rpg = REAL(ig, dp)*dr(1) - roffset(1)
77         t_exp_min_1 = t_exp_min_1*t_exp_min_2*t_exp_1
78         t_exp_min_2 = t_exp_min_2*t_exp_2
79         pg = t_exp_min_1
80         !pg  = EXP(-zetp*rpg**2)
81         DO icoef = 0, lp
82            pol_x(icoef, ig) = pg
83            pg = pg*(rpg)
84         ENDDO
85
86         rpg = REAL(1 - ig, dp)*dr(1) - roffset(1)
87         t_exp_plus_1 = t_exp_plus_1*t_exp_plus_2*t_exp_1
88         t_exp_plus_2 = t_exp_plus_2*t_exp_2
89         pg = t_exp_plus_1
90         ! pg  = EXP(-zetp*rpg**2)
91         DO icoef = 0, lp
92            pol_x(icoef, 1 - ig) = pg
93            pg = pg*(rpg)
94         ENDDO
95      ENDDO
96
97