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