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