1      Subroutine hfmke(Axyz,Bxyz,alpha,ES,E,pf,nd,NPP,MXD,La,Lb)
2c $Id$
3
4      Implicit real*8 (a-h,o-z)
5      Implicit integer (i-n)
6
7c--> Cartesian Coordinates
8
9      Dimension Axyz(3),Bxyz(3)
10
11c--> Exponents & Scaling Factors
12
13      Dimension alpha(2,NPP),ES(3,NPP)
14
15c--> Hermite Linear Expansion Coefficients
16
17      Dimension E(3,NPP,0:MXD,0:(La+Lb),0:La,0:Lb)
18
19c--> Scratch Space
20
21      Dimension pf(2,NPP)
22c
23c Define the linear expansion coefficients for the product of two CGF in
24c terms of HGF.
25c
26c   Recursion Formulas:
27c
28c      Ia+1,Ib;n            Ia,Ib;n            Ia,Ib;n           Ia,Ib;n
29c     E           = (1/2p) E        - (b/p) R E        + (Ip+1) E
30c      Ip                   Ip-1               Ip                Ip+1
31c
32c      Ia,Ib+1;n            Ia,Ib;n            Ia,Ib;n           Ia,Ib;n
33c     E           = (1/2p) E        + (a/p) R E        + (Ip+1) E
34c      Ip                   Ip-1               Ip                Ip+1
35c
36c   Initial Values:
37c
38c       0,0;0            0,0;0            0,0;0
39c     Ex      = ESx,   Ey      = ESy,   Ez      = ESz,
40c       0                0                0
41c
42c                               /  PI   \ 1/2      /   a b     2 \
43c      where typically   ESx = | ------- |     EXP| - -----  Rx   |
44c                               \ a + b /          \  a + b      /
45c
46c N.B. The prefactors for the overlap distribution of the two CGF
47c are typically used to scale the expansion coefficients. Products of
48c contraction coefficients may also be incorporated.
49c
50c   Indices for E(k,m,n,Ip,Ia,Ib):
51c
52c     1 [  k  ] 3            Cartesian Components {X,Y,Z}
53c     1 [  m  ] NPP          CGF Product Pair Index
54c     0 [  nd ] MXD          Derivative Index
55c     0 [  Ip ] Ia+Ib        Angular Momentum Index for the HGF
56c     0 [  Ia ] La           Angular Momentum Index for the CGF on Center "A"
57c     0 [  Ib ] Lb           Angular Momentum Index for the CGF on Center "B"
58c
59c******************************************************************************
60
61c Compute the vector between the centers.
62
63      Rx = Axyz(1) - Bxyz(1)
64      Ry = Axyz(2) - Bxyz(2)
65      Rz = Axyz(3) - Bxyz(3)
66
67      if( nd.eq.0 )then
68
69c Define the expansion coefficients for the products of two CGF.
70
71c Define E(Ip,Ia,Ib) for Ip=0, Ia=0, Ib=0.
72
73       do 100 m = 1,NPP
74        E(1,m,nd,0,0,0) = ES(1,m)
75        E(2,m,nd,0,0,0) = ES(2,m)
76        E(3,m,nd,0,0,0) = ES(3,m)
77  100  continue
78
79c Compute the prefactor for the 1st term of the recursion formulas.
80
81       if( La.gt.0 .or. Lb.gt.0 )then
82        do 110 m = 1,NPP
83         pf(1,m) = 0.5D0/( alpha(1,m) + alpha(2,m) )
84  110   continue
85       end if
86
87c Define E(Ip,Ia,Ib) for Ip = 0,Ia+Ib; Ia = 1,La; Ib = 0.
88
89       if( La.gt.0 )then
90
91c Compute the prefactor for the 2nd term of the recursion formula.
92
93        do 200 m = 1,NPP
94         pf(2,m) = -2.D0*alpha(2,m)*pf(1,m)
95  200   continue
96
97c            ===>   Ip = 0,Ia; Ia = 1; Ib = 0   <===
98
99        do 210 m = 1,NPP
100
101         Ex2 = Rx*E(1,m,nd,0,0,0)
102         Ey2 = Ry*E(2,m,nd,0,0,0)
103         Ez2 = Rz*E(3,m,nd,0,0,0)
104
105         E(1,m,nd,0,1,0) = pf(2,m)*Ex2
106         E(2,m,nd,0,1,0) = pf(2,m)*Ey2
107         E(3,m,nd,0,1,0) = pf(2,m)*Ez2
108
109         Ex1 = E(1,m,nd,0,0,0)
110         Ey1 = E(2,m,nd,0,0,0)
111         Ez1 = E(3,m,nd,0,0,0)
112
113         E(1,m,nd,1,1,0) = pf(1,m)*Ex1
114         E(2,m,nd,1,1,0) = pf(1,m)*Ey1
115         E(3,m,nd,1,1,0) = pf(1,m)*Ez1
116
117  210   continue
118
119        do 260 Ia = 2,La
120
121c            ===>   Ip = 0; Ia = 2,La; Ib = 0   <===
122
123         do 220 m = 1,NPP
124
125          Ex2 = Rx*E(1,m,nd,0,Ia-1,0)
126          Ex3 =    E(1,m,nd,1,Ia-1,0)
127
128          Ey2 = Ry*E(2,m,nd,0,Ia-1,0)
129          Ey3 =    E(2,m,nd,1,Ia-1,0)
130
131          Ez2 = Rz*E(3,m,nd,0,Ia-1,0)
132          Ez3 =    E(3,m,nd,1,Ia-1,0)
133
134          E(1,m,nd,0,Ia,0) = pf(2,m)*Ex2 + Ex3
135          E(2,m,nd,0,Ia,0) = pf(2,m)*Ey2 + Ey3
136          E(3,m,nd,0,Ia,0) = pf(2,m)*Ez2 + Ez3
137
138  220    continue
139
140c            ===>   Ip = 1,Ia-2; Ia = 2,La; Ib = 0   <===
141
142          do 240 Ip = 1,Ia-2
143
144           do m = 1,NPP
145
146            Ex1 =    E(1,m,nd,Ip-1,Ia-1,0)
147            Ex2 = Rx*E(1,m,nd,Ip  ,Ia-1,0)
148            Ex3 =    E(1,m,nd,Ip+1,Ia-1,0)
149
150            Ey1 =    E(2,m,nd,Ip-1,Ia-1,0)
151            Ey2 = Ry*E(2,m,nd,Ip  ,Ia-1,0)
152            Ey3 =    E(2,m,nd,Ip+1,Ia-1,0)
153
154            Ez1 =    E(3,m,nd,Ip-1,Ia-1,0)
155            Ez2 = Rz*E(3,m,nd,Ip  ,Ia-1,0)
156            Ez3 =    E(3,m,nd,Ip+1,Ia-1,0)
157
158            E(1,m,nd,Ip,Ia,0) = pf(1,m)*Ex1 + pf(2,m)*Ex2 + (Ip+1)*Ex3
159            E(2,m,nd,Ip,Ia,0) = pf(1,m)*Ey1 + pf(2,m)*Ey2 + (Ip+1)*Ey3
160            E(3,m,nd,Ip,Ia,0) = pf(1,m)*Ez1 + pf(2,m)*Ez2 + (Ip+1)*Ez3
161
162         enddo
163
164  240    continue
165
166c            ===>   Ip = Ia-1,Ia; Ia = 2,La; Ib = 0   <===
167
168         Ip = Ia-1
169cbreakakakakaka
170         do  m = 1,NPP
171            Ex1=E(1,m,nd,Ip,Ia-1,0)
172            Ey1=E(2,m,nd,Ip,Ia-1,0)
173            Ez1=E(3,m,nd,Ip,Ia-1,0)
174            E(1,m,nd,Ip+1,Ia,0) = pf(1,m)*Ex1
175            E(2,m,nd,Ip+1,Ia,0) = pf(1,m)*Ey1
176            E(3,m,nd,Ip+1,Ia,0) = pf(1,m)*Ez1
177            E(1,m,nd,Ip,Ia,0) =   pf(2,m)*Rx*Ex1
178            E(2,m,nd,Ip,Ia,0) =   pf(2,m)*Ry*Ey1
179            E(3,m,nd,Ip,Ia,0) =   pf(2,m)*Rz*Ez1
180
181         enddo
182
183         do  m = 1,NPP
184          E(1,m,nd,Ip,Ia,0) = E(1,m,nd,Ip,Ia,0) +
185     +           pf(1,m)*E(1,m,nd,Ip-1,Ia-1,0)
186          E(2,m,nd,Ip,Ia,0) = E(2,m,nd,Ip,Ia,0) +
187     +         pf(1,m)*E(2,m,nd,Ip-1,Ia-1,0)
188          E(3,m,nd,Ip,Ia,0) = E(3,m,nd,Ip,Ia,0) +
189     +         pf(1,m)*E(3,m,nd,Ip-1,Ia-1,0)
190
191
192        enddo
193
194  260   continue
195
196       end if
197
198c Define E(Ip,Ia,Ib) for Ip=0,Ia+Ib, Ia=0,La, Ib=1,Lb.
199
200       if( Lb.gt.0 )then
201
202c Compute the prefactor for the 2nd term of the recursion formula.
203
204        do 300 m = 1,NPP
205         pf(2,m) = 2.D0*alpha(1,m)*pf(1,m)
206  300   continue
207
208c    ===>   Ip = 0,Ia+Ib; Ia = 0; Ib = 1   <===
209
210        do 310 m = 1,NPP
211
212         Ex2 = Rx*E(1,m,nd,0,0,0)
213         Ey2 = Ry*E(2,m,nd,0,0,0)
214         Ez2 = Rz*E(3,m,nd,0,0,0)
215
216         E(1,m,nd,0,0,1) = pf(2,m)*Ex2
217         E(2,m,nd,0,0,1) = pf(2,m)*Ey2
218         E(3,m,nd,0,0,1) = pf(2,m)*Ez2
219
220         Ex1 = E(1,m,nd,0,0,0)
221         Ey1 = E(2,m,nd,0,0,0)
222         Ez1 = E(3,m,nd,0,0,0)
223
224         E(1,m,nd,1,0,1) = pf(1,m)*Ex1
225         E(2,m,nd,1,0,1) = pf(1,m)*Ey1
226         E(3,m,nd,1,0,1) = pf(1,m)*Ez1
227
228  310   continue
229
230        do 370 Ib = 1,Lb
231
232         if( Ib.eq.1 )then
233          Ia1 = 1
234         else
235          Ia1 = 0
236         end if
237
238         do 360 Ia = Ia1,La
239
240c    ===>   Ip = 0; Ia = Ia1,La; Ib = 1,Lb   <===
241
242          do 320 m = 1,NPP
243
244           Ex2 = Rx*E(1,m,nd,0,Ia,Ib-1)
245           Ex3 =    E(1,m,nd,1,Ia,Ib-1)
246
247           Ey2 = Ry*E(2,m,nd,0,Ia,Ib-1)
248           Ey3 =    E(2,m,nd,1,Ia,Ib-1)
249
250           Ez2 = Rz*E(3,m,nd,0,Ia,Ib-1)
251           Ez3 =    E(3,m,nd,1,Ia,Ib-1)
252
253           E(1,m,nd,0,Ia,Ib) = pf(2,m)*Ex2 + Ex3
254           E(2,m,nd,0,Ia,Ib) = pf(2,m)*Ey2 + Ey3
255           E(3,m,nd,0,Ia,Ib) = pf(2,m)*Ez2 + Ez3
256
257  320     continue
258
259c    ===>   Ip = 1,Ia+Ib-2; Ia = Ia1,La; Ib = 1,Lb   <===
260
261          do 340 Ip = 1,Ia+Ib-2
262
263           do 330 m = 1,NPP
264
265            Ex1 =    E(1,m,nd,Ip-1,Ia,Ib-1)
266            Ex2 = Rx*E(1,m,nd,Ip  ,Ia,Ib-1)
267            Ex3 =    E(1,m,nd,Ip+1,Ia,Ib-1)
268
269            Ey1 =    E(2,m,nd,Ip-1,Ia,Ib-1)
270            Ey2 = Ry*E(2,m,nd,Ip  ,Ia,Ib-1)
271            Ey3 =    E(2,m,nd,Ip+1,Ia,Ib-1)
272
273            Ez1 =    E(3,m,nd,Ip-1,Ia,Ib-1)
274            Ez2 = Rz*E(3,m,nd,Ip  ,Ia,Ib-1)
275            Ez3 =    E(3,m,nd,Ip+1,Ia,Ib-1)
276
277            E(1,m,nd,Ip,Ia,Ib) = pf(1,m)*Ex1 + pf(2,m)*Ex2 + (Ip+1)*Ex3
278            E(2,m,nd,Ip,Ia,Ib) = pf(1,m)*Ey1 + pf(2,m)*Ey2 + (Ip+1)*Ey3
279            E(3,m,nd,Ip,Ia,Ib) = pf(1,m)*Ez1 + pf(2,m)*Ez2 + (Ip+1)*Ez3
280
281  330      continue
282
283  340     continue
284
285c    ===>   Ip = Ia+Ib-1,Ia+Ib; Ia = Ia1,La; Ib = 1,Lb   <===
286
287          Ip = Ia+Ib-1
288
289          do 350 m = 1,NPP
290
291           Ex1 =    E(1,m,nd,Ip-1,Ia,Ib-1)
292           Ex2 = Rx*E(1,m,nd,Ip  ,Ia,Ib-1)
293
294           Ey1 =    E(2,m,nd,Ip-1,Ia,Ib-1)
295           Ey2 = Ry*E(2,m,nd,Ip  ,Ia,Ib-1)
296
297           Ez1 =    E(3,m,nd,Ip-1,Ia,Ib-1)
298           Ez2 = Rz*E(3,m,nd,Ip  ,Ia,Ib-1)
299
300           E(1,m,nd,Ip,Ia,Ib) = pf(1,m)*Ex1 + pf(2,m)*Ex2
301           E(2,m,nd,Ip,Ia,Ib) = pf(1,m)*Ey1 + pf(2,m)*Ey2
302           E(3,m,nd,Ip,Ia,Ib) = pf(1,m)*Ez1 + pf(2,m)*Ez2
303
304           Ex1 = E(1,m,nd,Ip,Ia,Ib-1)
305           Ey1 = E(2,m,nd,Ip,Ia,Ib-1)
306           Ez1 = E(3,m,nd,Ip,Ia,Ib-1)
307
308           E(1,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ex1
309           E(2,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ey1
310           E(3,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ez1
311
312  350     continue
313
314  360    continue
315
316  370   continue
317
318       end if
319
320      else
321
322c Define the expansion coefficients for derivatives of the products of two CGF.
323
324c Define E(Ip,Ia,Ib) for Ip=0, Ia=0, Ib=0.
325
326       if( nd.eq.1 )then
327
328        do 1100 m = 1,NPP
329
330         c = -2.D0*( alpha(1,m)*alpha(2,m)/( alpha(1,m) + alpha(2,m) ) )
331
332         E(1,m,nd,0,0,0) = c*Rx*E(1,m,nd-1,0,0,0)
333         E(2,m,nd,0,0,0) = c*Ry*E(2,m,nd-1,0,0,0)
334         E(3,m,nd,0,0,0) = c*Rz*E(3,m,nd-1,0,0,0)
335
336 1100   continue
337
338       else
339
340        do 1105 m = 1,NPP
341
342         c = -2.D0*( alpha(1,m)*alpha(2,m)/( alpha(1,m) + alpha(2,m) ) )
343
344         E1 = c*(Rx*E(1,m,nd-1,0,0,0) + (nd-1)*E(1,m,nd-2,0,0,0))
345         E2 = c*(Ry*E(2,m,nd-1,0,0,0) + (nd-1)*E(2,m,nd-2,0,0,0))
346         E3 = c*(Rz*E(3,m,nd-1,0,0,0) + (nd-1)*E(3,m,nd-2,0,0,0))
347
348         E(1,m,nd,0,0,0) = E1
349         E(2,m,nd,0,0,0) = E2
350         E(3,m,nd,0,0,0) = E3
351
352 1105   continue
353
354       end if
355
356c Compute the prefactor for the 1st term of recursion formulas.
357
358       if( La.gt.0 .or. Lb.gt.0 )then
359        do 1110 m = 1,NPP
360         pf(1,m) = 0.5D0/( alpha(1,m) + alpha(2,m) )
361 1110   continue
362       end if
363
364c Define E(Ip,Ia,Ib) for Ip = 0,Ia+Ib; Ia = 1,La; Ib = 0.
365
366       if( La.gt.0 )then
367
368c Compute the prefactor for the 2nd term of the recursion formula.
369
370        do 1200 m = 1,NPP
371         pf(2,m) = -2.D0*alpha(2,m)*pf(1,m)
372 1200   continue
373
374c            ===>   Ip = 0,Ia; Ia = 1; Ib = 0   <===
375
376        do 1210 m = 1,NPP
377
378         Ex2 = Rx*E(1,m,nd,0,0,0) + nd*E(1,m,nd-1,0,0,0)
379         Ey2 = Ry*E(2,m,nd,0,0,0) + nd*E(2,m,nd-1,0,0,0)
380         Ez2 = Rz*E(3,m,nd,0,0,0) + nd*E(3,m,nd-1,0,0,0)
381
382         E(1,m,nd,0,1,0) = pf(2,m)*Ex2
383         E(2,m,nd,0,1,0) = pf(2,m)*Ey2
384         E(3,m,nd,0,1,0) = pf(2,m)*Ez2
385
386         Ex1 = E(1,m,nd,0,0,0)
387         Ey1 = E(2,m,nd,0,0,0)
388         Ez1 = E(3,m,nd,0,0,0)
389
390         E(1,m,nd,1,1,0) = pf(1,m)*Ex1
391         E(2,m,nd,1,1,0) = pf(1,m)*Ey1
392         E(3,m,nd,1,1,0) = pf(1,m)*Ez1
393
394 1210   continue
395
396        do 1260 Ia = 2,La
397
398c            ===>   Ip = 0; Ia = 2,La; Ib = 0   <===
399
400         do 1220 m = 1,NPP
401
402          Ex2 = Rx*E(1,m,nd,0,Ia-1,0) + nd*E(1,m,nd-1,0,Ia-1,0)
403          Ex3 =    E(1,m,nd,1,Ia-1,0)
404
405          Ey2 = Ry*E(2,m,nd,0,Ia-1,0) + nd*E(2,m,nd-1,0,Ia-1,0)
406          Ey3 =    E(2,m,nd,1,Ia-1,0)
407
408          Ez2 = Rz*E(3,m,nd,0,Ia-1,0) + nd*E(3,m,nd-1,0,Ia-1,0)
409          Ez3 =    E(3,m,nd,1,Ia-1,0)
410
411          E(1,m,nd,0,Ia,0) = pf(2,m)*Ex2 + Ex3
412          E(2,m,nd,0,Ia,0) = pf(2,m)*Ey2 + Ey3
413          E(3,m,nd,0,Ia,0) = pf(2,m)*Ez2 + Ez3
414
415 1220    continue
416
417c            ===>   Ip = 1,Ia-2; Ia = 2,La; Ib = 0   <===
418
419          do 1240 Ip = 1,Ia-2
420
421           do 1230 m = 1,NPP
422
423            Ex1 =    E(1,m,nd,Ip-1,Ia-1,0)
424            Ex2 = Rx*E(1,m,nd,Ip  ,Ia-1,0) + nd*E(1,m,nd-1,Ip,Ia-1,0)
425            Ex3 =    E(1,m,nd,Ip+1,Ia-1,0)
426
427            Ey1 =    E(2,m,nd,Ip-1,Ia-1,0)
428            Ey2 = Ry*E(2,m,nd,Ip  ,Ia-1,0) + nd*E(2,m,nd-1,Ip,Ia-1,0)
429            Ey3 =    E(2,m,nd,Ip+1,Ia-1,0)
430
431            Ez1 =    E(3,m,nd,Ip-1,Ia-1,0)
432            Ez2 = Rz*E(3,m,nd,Ip  ,Ia-1,0) + nd*E(3,m,nd-1,Ip,Ia-1,0)
433            Ez3 =    E(3,m,nd,Ip+1,Ia-1,0)
434
435            E(1,m,nd,Ip,Ia,0) = pf(1,m)*Ex1 + pf(2,m)*Ex2 + (Ip+1)*Ex3
436            E(2,m,nd,Ip,Ia,0) = pf(1,m)*Ey1 + pf(2,m)*Ey2 + (Ip+1)*Ey3
437            E(3,m,nd,Ip,Ia,0) = pf(1,m)*Ez1 + pf(2,m)*Ez2 + (Ip+1)*Ez3
438
439 1230     continue
440
441 1240    continue
442
443c            ===>   Ip = Ia-1,Ia; Ia = 2,La; Ib = 0   <===
444
445         Ip = Ia-1
446
447cbreak2
448         do m = 1,NPP
449
450          Ex1 = E(1,m,nd,Ip,Ia-1,0)
451          Ey1 = E(2,m,nd,Ip,Ia-1,0)
452          Ez1 = E(3,m,nd,Ip,Ia-1,0)
453
454          E(1,m,nd,Ip,Ia,0) =  pf(2,m) * (Rx*Ex1 +
455     +         nd*E(1,m,nd-1,Ip,Ia-1,0))
456          E(2,m,nd,Ip,Ia,0) =  pf(2,m) * (Ry*Ey1 +
457     +         nd*E(2,m,nd-1,Ip,Ia-1,0))
458          E(3,m,nd,Ip,Ia,0) =  pf(2,m) * (Rz*Ez1 +
459     +         nd*E(3,m,nd-1,Ip,Ia-1,0))
460
461
462          E(1,m,nd,Ip+1,Ia,0) = pf(1,m)*Ex1
463          E(2,m,nd,Ip+1,Ia,0) = pf(1,m)*Ey1
464          E(3,m,nd,Ip+1,Ia,0) = pf(1,m)*Ez1
465
466       enddo
467         do  m = 1,NPP
468
469            Ex1 =    E(1,m,nd,Ip-1,Ia-1,0)
470            Ey1 =    E(2,m,nd,Ip-1,Ia-1,0)
471            Ez1 =    E(3,m,nd,Ip-1,Ia-1,0)
472
473
474            E(1,m,nd,Ip,Ia,0) = E(1,m,nd,Ip,Ia,0)+pf(1,m)*Ex1
475            E(2,m,nd,Ip,Ia,0) = E(2,m,nd,Ip,Ia,0)+pf(1,m)*Ey1
476            E(3,m,nd,Ip,Ia,0) = E(3,m,nd,Ip,Ia,0)+pf(1,m)*Ez1
477
478
479         enddo
480
481 1260   continue
482
483       end if
484
485c Define E(Ip,Ia,Ib) for Ip=0,Ia+Ib, Ia=0,La, Ib=1,Lb.
486
487       if( Lb.gt.0 )then
488
489c Compute the prefactor for the 2nd term of the recursion formula.
490
491        do 1300 m = 1,NPP
492         pf(2,m) = 2.D0*alpha(1,m)*pf(1,m)
493 1300   continue
494
495c    ===>   Ip = 0,Ia+Ib; Ia = 0; Ib = 1   <===
496
497        do 1310 m = 1,NPP
498
499         Ex2 = Rx*E(1,m,nd,0,0,0) + nd*E(1,m,nd-1,0,0,0)
500         Ey2 = Ry*E(2,m,nd,0,0,0) + nd*E(2,m,nd-1,0,0,0)
501         Ez2 = Rz*E(3,m,nd,0,0,0) + nd*E(3,m,nd-1,0,0,0)
502
503         E(1,m,nd,0,0,1) = pf(2,m)*Ex2
504         E(2,m,nd,0,0,1) = pf(2,m)*Ey2
505         E(3,m,nd,0,0,1) = pf(2,m)*Ez2
506
507         Ex1 = E(1,m,nd,0,0,0)
508         Ey1 = E(2,m,nd,0,0,0)
509         Ez1 = E(3,m,nd,0,0,0)
510
511         E(1,m,nd,1,0,1) = pf(1,m)*Ex1
512         E(2,m,nd,1,0,1) = pf(1,m)*Ey1
513         E(3,m,nd,1,0,1) = pf(1,m)*Ez1
514
515 1310   continue
516
517        do 1370 Ib = 1,Lb
518
519         if( Ib.eq.1 )then
520          Ia1 = 1
521         else
522          Ia1 = 0
523         end if
524
525         do 1360 Ia = Ia1,La
526
527c    ===>   Ip = 0; Ia = Ia1,La; Ib = 1,Lb   <===
528
529          do 1320 m = 1,NPP
530
531           Ex2 = Rx*E(1,m,nd,0,Ia,Ib-1) + nd*E(1,m,nd-1,0,Ia,Ib-1)
532           Ex3 =    E(1,m,nd,1,Ia,Ib-1)
533
534           Ey2 = Ry*E(2,m,nd,0,Ia,Ib-1) + nd*E(2,m,nd-1,0,Ia,Ib-1)
535           Ey3 =    E(2,m,nd,1,Ia,Ib-1)
536
537           Ez2 = Rz*E(3,m,nd,0,Ia,Ib-1) + nd*E(3,m,nd-1,0,Ia,Ib-1)
538           Ez3 =    E(3,m,nd,1,Ia,Ib-1)
539
540           E(1,m,nd,0,Ia,Ib) = pf(2,m)*Ex2 + Ex3
541           E(2,m,nd,0,Ia,Ib) = pf(2,m)*Ey2 + Ey3
542           E(3,m,nd,0,Ia,Ib) = pf(2,m)*Ez2 + Ez3
543
544 1320     continue
545
546c    ===>   Ip = 1,Ia+Ib-2; Ia = Ia1,La; Ib = 1,Lb   <===
547
548          do 1340 Ip = 1,Ia+Ib-2
549
550           do 1330 m = 1,NPP
551
552            Ex1 =    E(1,m,nd,Ip-1,Ia,Ib-1)
553            Ex2 = Rx*E(1,m,nd,Ip  ,Ia,Ib-1) + nd*E(1,m,nd-1,Ip,Ia,Ib-1)
554            Ex3 =    E(1,m,nd,Ip+1,Ia,Ib-1)
555
556            Ey1 =    E(2,m,nd,Ip-1,Ia,Ib-1)
557            Ey2 = Ry*E(2,m,nd,Ip  ,Ia,Ib-1) + nd*E(2,m,nd-1,Ip,Ia,Ib-1)
558            Ey3 =    E(2,m,nd,Ip+1,Ia,Ib-1)
559
560            Ez1 =    E(3,m,nd,Ip-1,Ia,Ib-1)
561            Ez2 = Rz*E(3,m,nd,Ip  ,Ia,Ib-1) + nd*E(3,m,nd-1,Ip,Ia,Ib-1)
562            Ez3 =    E(3,m,nd,Ip+1,Ia,Ib-1)
563
564            E(1,m,nd,Ip,Ia,Ib) = pf(1,m)*Ex1 + pf(2,m)*Ex2 + (Ip+1)*Ex3
565            E(2,m,nd,Ip,Ia,Ib) = pf(1,m)*Ey1 + pf(2,m)*Ey2 + (Ip+1)*Ey3
566            E(3,m,nd,Ip,Ia,Ib) = pf(1,m)*Ez1 + pf(2,m)*Ez2 + (Ip+1)*Ez3
567
568 1330      continue
569
570 1340     continue
571
572c    ===>   Ip = Ia+Ib-1,Ia+Ib; Ia = Ia1,La; Ib = 1,Lb   <===
573
574          Ip = Ia+Ib-1
575
576          do 1350 m = 1,NPP
577
578           Ex1 =    E(1,m,nd,Ip-1,Ia,Ib-1)
579           Ex2 = Rx*E(1,m,nd,Ip  ,Ia,Ib-1) + nd*E(1,m,nd-1,Ip,Ia,Ib-1)
580
581           Ey1 =    E(2,m,nd,Ip-1,Ia,Ib-1)
582           Ey2 = Ry*E(2,m,nd,Ip  ,Ia,Ib-1) + nd*E(2,m,nd-1,Ip,Ia,Ib-1)
583
584           Ez1 =    E(3,m,nd,Ip-1,Ia,Ib-1)
585           Ez2 = Rz*E(3,m,nd,Ip  ,Ia,Ib-1) + nd*E(3,m,nd-1,Ip,Ia,Ib-1)
586
587           E(1,m,nd,Ip,Ia,Ib) = pf(1,m)*Ex1 + pf(2,m)*Ex2
588           E(2,m,nd,Ip,Ia,Ib) = pf(1,m)*Ey1 + pf(2,m)*Ey2
589           E(3,m,nd,Ip,Ia,Ib) = pf(1,m)*Ez1 + pf(2,m)*Ez2
590
591           Ex1 = E(1,m,nd,Ip,Ia,Ib-1)
592           Ey1 = E(2,m,nd,Ip,Ia,Ib-1)
593           Ez1 = E(3,m,nd,Ip,Ia,Ib-1)
594
595           E(1,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ex1
596           E(2,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ey1
597           E(3,m,nd,Ip+1,Ia,Ib) = pf(1,m)*Ez1
598
599 1350     continue
600
601 1360    continue
602
603 1370   continue
604
605       end if
606
607      end if
608
609      end
610