35 |
|
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
36 |
|
|
37 |
|
typedef struct { |
38 |
< |
FVECT r_i, r_i1, e_i; |
39 |
< |
double nf, I1, I2, J2; |
38 |
> |
FVECT r_i, r_i1, e_i, rI2_eJ2; |
39 |
> |
double nf, I1, I2; |
40 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
41 |
|
|
42 |
|
|
138 |
|
static void |
139 |
|
comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop) |
140 |
|
{ |
141 |
< |
FVECT v1; |
142 |
< |
double dot_e, dot_er, dot_r, dot_r1; |
141 |
> |
FVECT vcp; |
142 |
> |
double dot_e, dot_er, dot_r, dot_r1, J2; |
143 |
> |
int i; |
144 |
|
|
145 |
|
VSUB(ftp->r_i, ap0, rop); |
146 |
|
VSUB(ftp->r_i1, ap1, rop); |
147 |
|
VSUB(ftp->e_i, ap1, ap0); |
148 |
< |
VCROSS(v1, ftp->e_i, ftp->r_i); |
149 |
< |
ftp->nf = 1.0/DOT(v1,v1); |
148 |
> |
VCROSS(vcp, ftp->e_i, ftp->r_i); |
149 |
> |
ftp->nf = 1.0/DOT(vcp,vcp); |
150 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
151 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
152 |
|
dot_r = DOT(ftp->r_i,ftp->r_i); |
155 |
|
sqrt( ftp->nf ); |
156 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + |
157 |
|
dot_e*ftp->I1 )*0.5*ftp->nf; |
158 |
< |
ftp->J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - |
159 |
< |
dot_er/dot_e*ftp->I2; |
158 |
> |
J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - dot_er/dot_e*ftp->I2; |
159 |
> |
for (i = 3; i--; ) |
160 |
> |
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
161 |
|
} |
162 |
|
|
163 |
|
|
178 |
|
static void |
179 |
|
comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) |
180 |
|
{ |
181 |
< |
FVECT v1, v2; |
181 |
> |
FVECT vcp; |
182 |
|
FVECT m1[3], m2[3], m3[3], m4[3]; |
183 |
|
double d1, d2, d3, d4; |
184 |
|
double I3, J3, K3; |
193 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
194 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
195 |
|
/* intermediate matrices */ |
196 |
< |
VCROSS(v1, nrm, ftp->e_i); |
197 |
< |
for (j = 3; j--; ) |
196 |
< |
v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; |
197 |
< |
compose_matrix(m1, v1, v2); |
196 |
> |
VCROSS(vcp, nrm, ftp->e_i); |
197 |
> |
compose_matrix(m1, vcp, ftp->rI2_eJ2); |
198 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
199 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
200 |
|
compose_matrix(m4, ftp->r_i, ftp->e_i); |
201 |
< |
VCROSS(v1, ftp->r_i, ftp->e_i); |
202 |
< |
d1 = DOT(nrm, v1); |
201 |
> |
VCROSS(vcp, ftp->r_i, ftp->e_i); |
202 |
> |
d1 = DOT(nrm, vcp); |
203 |
|
d2 = -d1*ftp->I2; |
204 |
|
d1 *= 2.0; |
205 |
|
for (i = 3; i--; ) /* final matrix sum */ |
251 |
|
f1 = 2.0*DOT(nrm, vcp); |
252 |
|
VCROSS(vcp, nrm, ftp->e_i); |
253 |
|
for (i = 3; i--; ) |
254 |
< |
grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + |
255 |
< |
f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) ); |
254 |
> |
grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
255 |
|
} |
256 |
|
|
257 |
|
|