! wifin from ai 10.2Dr subroutine poten(xp,v,r1,r2,th) integer,parameter :: npropin = 246 double precision,intent(in) :: r1,r2,th double precision,intent(out):: v double precision,intent(in) :: xp(1:npropin) double precision :: reoh,thetae,b1,roh,alphaoh double precision :: phh2,t0,ut,x0,ut2,x02,xs1,xs2,xst,rs,rm,rr1,rr2,xep1,xep2,xep3 double precision :: rhh,vhh,vpb1,vpb2,v0,vp1,vp2,vp3,vp,vps1,vps2,y1,y2,y12,y22,voh1,voh2 integer(4) :: i ! reoh=0.95860d0 ! thetae=104.3475d0 thetae=104.4800d0 b1=2.15d0 ! roh=0.9578d0 alphaoh=0.0d0 phh2=6.20164303995d0 thetae=thetae*.314159265358979312d01*.00555555555555555555d0 xs1=(r1+r2)*0.5d0-reoh xs2=(r1-r2)*0.5d0 xst=dcos(th)-dcos(thetae) rr1=r1-reoh rr2=r2-reoh alpha=2.2668d0 ! alpha=1.5d0 xep1=(dexp(-2.d0*alpha*rr1)-2.d0*dexp(-alpha*rr1))* & (0.4389830771267D+05 & -0.8606164763375D+03*rr1 & +0.1746549380317D+05*rr1**2 & -0.3326053761996D+05*rr1**3 & +0.4988743744461D+04*rr1**4 & +0.8751979626215D+04*rr1**5 & -0.3770596827415D+04*rr1**6 & +0.4270854122634D+03*rr1**7) + 0.4389830771267D+05 xep2=(dexp(-2.d0*alpha*rr2)-2.d0*dexp(-alpha*rr2))* & (0.4389830771267D+05 & -0.8606164763375D+03*rr2 & +0.1746549380317D+05*rr2**2 & -0.3326053761996D+05*rr2**3 & +0.4988743744461D+04*rr2**4 & +0.8751979626215D+04*rr2**5 & -0.3770596827415D+04*rr2**6 & +0.4270854122634D+03*rr2**7) + 0.4389830771267D+05 xep3=dexp(-b1*((r1-reoh)**2+(r2-reoh)**2)) rhh=dsqrt(r1**2+r2**2-2.d0*r1*r2*dcos(th)) vhh=0.820894739261131734D+06*dexp(-phh2*rhh) v0 = xp( 1) *xs1**0 *xs2**0 *xst**0 vp1=+xp( 2) *xs1**0 *xs2**0 *xst**1 & +xp( 3) *xs1**1 *xs2**0 *xst**0 & +xp( 4) *xs1**0 *xs2**0 *xst**2 & +xp( 5) *xs1**0 *xs2**2 *xst**0 & +xp( 6) *xs1**1 *xs2**0 *xst**1 & +xp( 7) *xs1**2 *xs2**0 *xst**0 & +xp( 8) *xs1**0 *xs2**0 *xst**3 & +xp( 9) *xs1**0 *xs2**2 *xst**1 & +xp( 10) *xs1**1 *xs2**0 *xst**2 & +xp( 11) *xs1**1 *xs2**2 *xst**0 & +xp( 12) *xs1**2 *xs2**0 *xst**1 & +xp( 13) *xs1**3 *xs2**0 *xst**0 & +xp( 14) *xs1**0 *xs2**0 *xst**4 & +xp( 15) *xs1**0 *xs2**2 *xst**2 & +xp( 16) *xs1**0 *xs2**4 *xst**0 & +xp( 17) *xs1**1 *xs2**0 *xst**3 & +xp( 18) *xs1**1 *xs2**2 *xst**1 & +xp( 19) *xs1**2 *xs2**0 *xst**2 & +xp( 20) *xs1**2 *xs2**2 *xst**0 & +xp( 21) *xs1**3 *xs2**0 *xst**1 & +xp( 22) *xs1**4 *xs2**0 *xst**0 & +xp( 23) *xs1**0 *xs2**0 *xst**5 & +xp( 24) *xs1**0 *xs2**2 *xst**3 & +xp( 25) *xs1**0 *xs2**4 *xst**1 & +xp( 26) *xs1**1 *xs2**0 *xst**4 & +xp( 27) *xs1**1 *xs2**2 *xst**2 & +xp( 28) *xs1**1 *xs2**4 *xst**0 & +xp( 29) *xs1**2 *xs2**0 *xst**3 & +xp( 30) *xs1**2 *xs2**2 *xst**1 & +xp( 31) *xs1**3 *xs2**0 *xst**2 & +xp( 32) *xs1**3 *xs2**2 *xst**0 & +xp( 33) *xs1**4 *xs2**0 *xst**1 & +xp( 34) *xs1**5 *xs2**0 *xst**0 & +xp( 35) *xs1**0 *xs2**0 *xst**6 & +xp( 36) *xs1**0 *xs2**2 *xst**4 & +xp( 37) *xs1**0 *xs2**4 *xst**2 & +xp( 38) *xs1**0 *xs2**6 *xst**0 & +xp( 39) *xs1**1 *xs2**0 *xst**5 & +xp( 40) *xs1**1 *xs2**2 *xst**3 & +xp( 41) *xs1**1 *xs2**4 *xst**1 & +xp( 42) *xs1**2 *xs2**0 *xst**4 & +xp( 43) *xs1**2 *xs2**2 *xst**2 & +xp( 44) *xs1**2 *xs2**4 *xst**0 & +xp( 45) *xs1**3 *xs2**0 *xst**3 & +xp( 46) *xs1**3 *xs2**2 *xst**1 & +xp( 47) *xs1**4 *xs2**0 *xst**2 & +xp( 48) *xs1**4 *xs2**2 *xst**0 & +xp( 49) *xs1**5 *xs2**0 *xst**1 & +xp( 50) *xs1**6 *xs2**0 *xst**0 & +xp( 51) *xs1**0 *xs2**0 *xst**7 & +xp( 52) *xs1**0 *xs2**2 *xst**5 & +xp( 53) *xs1**0 *xs2**4 *xst**3 & +xp( 54) *xs1**0 *xs2**6 *xst**1 & +xp( 55) *xs1**1 *xs2**0 *xst**6 & +xp( 56) *xs1**1 *xs2**2 *xst**4 & +xp( 57) *xs1**1 *xs2**4 *xst**2 & +xp( 58) *xs1**1 *xs2**6 *xst**0 & +xp( 59) *xs1**2 *xs2**0 *xst**5 & +xp( 60) *xs1**2 *xs2**2 *xst**3 & +xp( 61) *xs1**2 *xs2**4 *xst**1 & +xp( 62) *xs1**3 *xs2**0 *xst**4 & +xp( 63) *xs1**3 *xs2**2 *xst**2 & +xp( 64) *xs1**3 *xs2**4 *xst**0 & +xp( 65) *xs1**4 *xs2**0 *xst**3 & +xp( 66) *xs1**4 *xs2**2 *xst**1 & +xp( 67) *xs1**5 *xs2**0 *xst**2 & +xp( 68) *xs1**5 *xs2**2 *xst**0 & +xp( 69) *xs1**6 *xs2**0 *xst**1 & +xp( 70) *xs1**7 *xs2**0 *xst**0 & +xp( 71) *xs1**0 *xs2**0 *xst**8 & +xp( 72) *xs1**0 *xs2**2 *xst**6 & +xp( 73) *xs1**0 *xs2**4 *xst**4 & +xp( 74) *xs1**0 *xs2**6 *xst**2 & +xp( 75) *xs1**0 *xs2**8 *xst**0 & +xp( 76) *xs1**1 *xs2**0 *xst**7 & +xp( 77) *xs1**1 *xs2**2 *xst**5 & +xp( 78) *xs1**1 *xs2**4 *xst**3 & +xp( 79) *xs1**1 *xs2**6 *xst**1 & +xp( 80) *xs1**2 *xs2**0 *xst**6 & +xp( 81) *xs1**2 *xs2**2 *xst**4 & +xp( 82) *xs1**2 *xs2**4 *xst**2 & +xp( 83) *xs1**2 *xs2**6 *xst**0 & +xp( 84) *xs1**3 *xs2**0 *xst**5 & +xp( 85) *xs1**3 *xs2**2 *xst**3 & +xp( 86) *xs1**3 *xs2**4 *xst**1 & +xp( 87) *xs1**4 *xs2**0 *xst**4 & +xp( 88) *xs1**4 *xs2**2 *xst**2 & +xp( 89) *xs1**4 *xs2**4 *xst**0 & +xp( 90) *xs1**5 *xs2**0 *xst**3 & +xp( 91) *xs1**5 *xs2**2 *xst**1 & +xp( 92) *xs1**6 *xs2**0 *xst**2 vp2=+xp( 93) *xs1**6 *xs2**2 *xst**0 & +xp( 94) *xs1**7 *xs2**0 *xst**1 & +xp( 95) *xs1**8 *xs2**0 *xst**0 & +xp( 96) *xs1**0 *xs2**0 *xst**9 & +xp( 97) *xs1**0 *xs2**2 *xst**7 & +xp( 98) *xs1**0 *xs2**4 *xst**5 & +xp( 99) *xs1**0 *xs2**6 *xst**3 & +xp(100) *xs1**0 *xs2**8 *xst**1 & +xp(101) *xs1**1 *xs2**0 *xst**8 & +xp(102) *xs1**1 *xs2**2 *xst**6 & +xp(103) *xs1**1 *xs2**4 *xst**4 & +xp(104) *xs1**1 *xs2**6 *xst**2 & +xp(105) *xs1**1 *xs2**8 *xst**0 & +xp(106) *xs1**2 *xs2**0 *xst**7 & +xp(107) *xs1**2 *xs2**2 *xst**5 & +xp(108) *xs1**2 *xs2**4 *xst**3 & +xp(109) *xs1**2 *xs2**6 *xst**1 & +xp(110) *xs1**3 *xs2**0 *xst**6 & +xp(111) *xs1**3 *xs2**2 *xst**4 & +xp(112) *xs1**3 *xs2**4 *xst**2 & +xp(113) *xs1**3 *xs2**6 *xst**0 & +xp(114) *xs1**4 *xs2**0 *xst**5 & +xp(115) *xs1**4 *xs2**2 *xst**3 & +xp(116) *xs1**4 *xs2**4 *xst**1 & +xp(117) *xs1**5 *xs2**0 *xst**4 & +xp(118) *xs1**5 *xs2**2 *xst**2 & +xp(119) *xs1**5 *xs2**4 *xst**0 & +xp(120) *xs1**6 *xs2**0 *xst**3 & +xp(121) *xs1**6 *xs2**2 *xst**1 & +xp(122) *xs1**7 *xs2**0 *xst**2 & +xp(123) *xs1**7 *xs2**2 *xst**0 & +xp(124) *xs1**8 *xs2**0 *xst**1 & +xp(125) *xs1**9 *xs2**0 *xst**0 & +xp(126) *xs1**0 *xs2**2 *xst**8 & +xp(127) *xs1**0 *xs2**4 *xst**6 & +xp(128) *xs1**0 *xs2**6 *xst**4 & +xp(129) *xs1**0 *xs2**8 *xst**2 & +xp(130) *xs1**1 *xs2**0 *xst**9 & +xp(131) *xs1**1 *xs2**2 *xst**7 & +xp(132) *xs1**1 *xs2**4 *xst**5 & +xp(133) *xs1**1 *xs2**6 *xst**3 & +xp(134) *xs1**1 *xs2**8 *xst**1 & +xp(135) *xs1**2 *xs2**0 *xst**8 & +xp(136) *xs1**2 *xs2**2 *xst**6 & +xp(137) *xs1**2 *xs2**4 *xst**4 & +xp(138) *xs1**2 *xs2**6 *xst**2 & +xp(139) *xs1**3 *xs2**0 *xst**7 & +xp(140) *xs1**3 *xs2**2 *xst**5 & +xp(141) *xs1**3 *xs2**4 *xst**3 & +xp(142) *xs1**3 *xs2**6 *xst**1 & +xp(143) *xs1**4 *xs2**0 *xst**6 & +xp(144) *xs1**4 *xs2**2 *xst**4 & +xp(145) *xs1**4 *xs2**4 *xst**2 & +xp(146) *xs1**5 *xs2**0 *xst**5 & +xp(147) *xs1**5 *xs2**2 *xst**3 & +xp(148) *xs1**5 *xs2**4 *xst**1 & +xp(149) *xs1**6 *xs2**0 *xst**4 & +xp(150) *xs1**6 *xs2**2 *xst**2 & +xp(151) *xs1**6 *xs2**4 *xst**0 & +xp(152) *xs1**7 *xs2**0 *xst**3 & +xp(153) *xs1**7 *xs2**2 *xst**1 & +xp(154) *xs1**8 *xs2**0 *xst**2 & +xp(155) *xs1**8 *xs2**2 *xst**0 & +xp(156) *xs1**9 *xs2**0 *xst**1 & +xp(157) *xs1**10*xs2**0 *xst**0 & +xp(158) *xs1**0 *xs2**2 *xst**9 & +xp(159) *xs1**0 *xs2**4 *xst**7 & +xp(160) *xs1**0 *xs2**6 *xst**5 & +xp(161) *xs1**0 *xs2**8 *xst**3 & +xp(162) *xs1**0 *xs2**10*xst**1 & +xp(163) *xs1**1 *xs2**0 *xst**10 & +xp(164) *xs1**1 *xs2**2 *xst**8 & +xp(165) *xs1**1 *xs2**4 *xst**6 & +xp(166) *xs1**1 *xs2**6 *xst**4 & +xp(167) *xs1**1 *xs2**8 *xst**2 & +xp(168) *xs1**1 *xs2**10*xst**0 & +xp(169) *xs1**2 *xs2**0 *xst**9 & +xp(170) *xs1**2 *xs2**2 *xst**7 & +xp(171) *xs1**2 *xs2**4 *xst**5 & +xp(172) *xs1**2 *xs2**6 *xst**3 & +xp(173) *xs1**2 *xs2**8 *xst**1 & +xp(174) *xs1**3 *xs2**0 *xst**8 vp3=+xp(175) *xs1**3 *xs2**2 *xst**6 & +xp(176) *xs1**3 *xs2**4 *xst**4 & +xp(177) *xs1**3 *xs2**6 *xst**2 & +xp(178) *xs1**3 *xs2**8 *xst**0 & +xp(179) *xs1**4 *xs2**0 *xst**7 & +xp(180) *xs1**4 *xs2**2 *xst**5 & +xp(181) *xs1**4 *xs2**4 *xst**3 & +xp(182) *xs1**4 *xs2**6 *xst**1 & +xp(183) *xs1**5 *xs2**0 *xst**6 & +xp(184) *xs1**5 *xs2**2 *xst**4 & +xp(185) *xs1**5 *xs2**4 *xst**2 & +xp(186) *xs1**5 *xs2**6 *xst**0 & +xp(187) *xs1**6 *xs2**0 *xst**5 & +xp(188) *xs1**6 *xs2**2 *xst**3 & +xp(189) *xs1**6 *xs2**4 *xst**1 & +xp(190) *xs1**7 *xs2**0 *xst**4 & +xp(191) *xs1**7 *xs2**2 *xst**2 & +xp(192) *xs1**8 *xs2**0 *xst**3 & +xp(193) *xs1**8 *xs2**2 *xst**1 & +xp(194) *xs1**9 *xs2**0 *xst**2 & +xp(195) *xs1**9 *xs2**2 *xst**0 & +xp(196) *xs1**10*xs2**0 *xst**1 & +xp(197) *xs1**11*xs2**0 *xst**0 & +xp(198) *xs1**0 *xs2**0 *xst**2 & +xp(199) *xs1**0 *xs2**2 *xst**0 & +xp(200) *xs1**0 *xs2**4 *xst**8 & +xp(201) *xs1**0 *xs2**6 *xst**6 & +xp(202) *xs1**0 *xs2**8 *xst**4 & +xp(203) *xs1**0 *xs2**10*xst**2 & +xp(204) *xs1**0 *xs2**12*xst**0 & +xp(205) *xs1**1 *xs2**0 *xst**11 & +xp(206) *xs1**1 *xs2**2 *xst**9 & +xp(207) *xs1**1 *xs2**4 *xst**7 & +xp(208) *xs1**1 *xs2**6 *xst**5 & +xp(209) *xs1**1 *xs2**8 *xst**3 & +xp(210) *xs1**1 *xs2**10*xst**1 & +xp(211) *xs1**2 *xs2**0 *xst**10 & +xp(212) *xs1**2 *xs2**2 *xst**8 & +xp(213) *xs1**2 *xs2**4 *xst**6 & +xp(214) *xs1**2 *xs2**6 *xst**4 & +xp(215) *xs1**2 *xs2**8 *xst**2 & +xp(216) *xs1**2 *xs2**10*xst**0 & +xp(217) *xs1**3 *xs2**0 *xst**9 & +xp(218) *xs1**3 *xs2**2 *xst**7 & +xp(219) *xs1**3 *xs2**4 *xst**5 & +xp(220) *xs1**3 *xs2**6 *xst**3 & +xp(221) *xs1**3 *xs2**8 *xst**1 & +xp(222) *xs1**4 *xs2**0 *xst**8 & +xp(223) *xs1**4 *xs2**2 *xst**6 & +xp(224) *xs1**4 *xs2**4 *xst**4 & +xp(225) *xs1**4 *xs2**6 *xst**2 & +xp(226) *xs1**4 *xs2**8 *xst**0 & +xp(227) *xs1**5 *xs2**0 *xst**7 & +xp(228) *xs1**5 *xs2**2 *xst**5 & +xp(229) *xs1**5 *xs2**4 *xst**3 & +xp(230) *xs1**5 *xs2**6 *xst**1 & +xp(231) *xs1**6 *xs2**0 *xst**6 & +xp(232) *xs1**6 *xs2**2 *xst**4 & +xp(233) *xs1**6 *xs2**4 *xst**2 & +xp(234) *xs1**6 *xs2**6 *xst**0 & +xp(235) *xs1**7 *xs2**0 *xst**5 & +xp(236) *xs1**7 *xs2**2 *xst**3 & +xp(237) *xs1**7 *xs2**4 *xst**1 & +xp(238) *xs1**8 *xs2**0 *xst**4 & +xp(239) *xs1**8 *xs2**2 *xst**2 & +xp(240) *xs1**8 *xs2**4 *xst**0 & +xp(241) *xs1**9 *xs2**0 *xst**3 & +xp(242) *xs1**9 *xs2**2 *xst**1 & +xp(243) *xs1**10*xs2**0 *xst**2 & +xp(244) *xs1**10*xs2**2 *xst**0 & +xp(245) *xs1**11*xs2**0 *xst**1 & +xp(246) *xs1**12*xs2**0 *xst**0 vp=vp1+vp2+vp3 v=v0+vp*xep3+vhh+xep1+xep2 ! v=v0+vp*xep3+xep1+xep2 end subroutine poten SUBROUTINE potv(V,R1,R2,xcos) IMPLICIT DOUBLE PRECISION (A-H,O-Y),logical(z) COMMON /MASS/ XMASS(3),G1,G2,xmassr(3) dimension pv(20),u(20) dimension ht(20) ! DATA RZ/0.9576257 /,RHO/75.48992362 / DATA RZ/0.958600d0 /,the/104.48d0 / DATA TOANG/0.5291772d0/ DATA X1/1.0D0/,X0/0.0D0/,TINY/9.0D-15/,X2/2.0D0/ DATA pv/1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0/ dimension rij(3) real*8 :: xp(1:300) integer :: i_t pi=dacos(-1.d0) IF (G1 .EQ. X0) THEN Q1 = R1 Q2 = R2 THETA = ACOS(XCOS) ELSE IF (G2 .EQ. X0) THEN XX = R1 * G1 YY = R1 * (X1 - G1) IF (R2 .EQ. X0 .OR. XCOS .GE. (X1 - TINY)) THEN Q1 = ABS(XX - R2) Q2 = (YY + R2) COST = -X1 ELSE IF (XCOS .LE. (TINY - X1)) THEN Q1 = (XX + R2) Q2 = ABS(YY + R2) COST = X1 ELSE Q1 = SQRT(XX*XX + R2*R2 - X2*XX*R2*XCOS) Q2 = SQRT(YY*YY + R2*R2 + X2*YY*R2*XCOS) COST = (Q1**2 + Q2**2 - R1**2) / (X2 * Q1 * Q2) ENDIF THETA = ACOS(COST) ELSE F1= X1/G1 F2= X1/G2 F12= X1 - F1*F2 P1= R1*(X1-F1)/(G2*F12) P2= R2*(X1-F2)/(G1*F12) S1= R1-P1 S2= R2-P2 Q1= SQRT(P1*P1 + S2*S2 + X2*P1*S2*XCOS)/(X1-G1) Q2= SQRT(P2*P2 + S1*S1 + X2*P2*S1*XCOS)/(X1-G2) Q3= SQRT(P1*P1 + P2*P2 - X2*P1*P2*XCOS) COST = (Q1*Q1 + Q2*Q2 - Q3*Q3)/(X2*Q1*Q2) THETA = ACOS(COST) ENDIF ! BODC if(q1.gt.4.0) q1=4.0 if(q2.gt.4.0) q2=4.0 CALL bodc16(V16,Q1,Q2,THETA) CALL bodc18(V18,Q1,Q2,THETA) xmaso16=15.990526 xmaso17=16.9947425 xmaso18=17.9947714 xmasd=2.013553214 xmash=1.00727647 vad18=(xmaso18*xmaso16*(v16-v18)/(xmaso18-xmaso16))/xmaso18+ *(2*xmash*(xmaso18*v18-xmaso16*v16)/(2*(xmaso18-xmaso16)))/xmash THETAeq = the*3.141592654/180.d0 Req= RZ/TOANG Y1 = R1 - Req Y2 = Cos(THETA) - Cos(THETAeq) Y3 = R2 - Req S1 = (Y1 + Y3)/2.0d0 S2 = Y2 S3 = (Y1 - Y3)/2.0d0 open(unit=32,status='old',form='formatted',file='pot.fit25.18') read(32,*) npropin do i = 1,npropin read(32,*) i_t, xp(i) end do close(unit=32) if (npropin>size(xp)) then write(6,"('wifin: Too many parameters in pot.fit: ', * i8,' max = ',i8)") npropin,size(xp) stop 'wifin: Too many parameters in pot.fit' endif ! att1=Q1*0.5291772d0 att2=Q2*0.5291772d0 ! call poten(xp(1:npropin),vp,att1,att2,THETA) v=vp/219474.624d0+Vad18 end SUBROUTINE BODC16(V,R1,R2,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FT(69), cv(69) DATA C0/0.0D0/,SCALE/1.0D-6/ data nv/69/ DATA ZERO/0.0D0/ DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/ DATA cv/2785.4109260806,-90.0463360612,48.5226385100, *201.6994168988,-19.4625405248,217.8639531832, *-48.9272565544,-76.0803927294,25.9675092231,-65.2206704082, *60.7284555026,-386.5173198637,-2.0352745279,3.5405842322, *-32.1830475213,126.8885218953,60.2606780157,134.6300248947, *-14.0991224631,560.8792812623,27.8464660146,-67.1464124186, *-181.5695164733,43.7287769150,-150.7468776877,-103.2652779984, *-159.6432245839,-24.8956181482,185.4825335893,-231.4775497546, *-51.9208759775,3.7066140140,-212.4422129941,207.1884490703, *383.1659239100,50.8641728660,35.1754939408,127.2280510484, *-154.4544312699,55.5787967758,282.6216945851,116.5606405651, *5.4433254144,-107.1461094167,-173.8895717556,-26.5859674990, *-560.8756697840,-237.7109212157,143.9462552048,-592.3478209334, *0.0,-198.6468835005,-19.9674473372,-14.1731270087,193.4510720304, *4.6347021028,32.9502486772,-221.1685318724,26.4090449111, *-268.432837385,-147.1422366151,133.5465868568,363.9554096142, *673.9006484856,214.9454642643,40.7735822438,65.2246188257, *173.0708970426,1.9795259929/ THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/toang S1 = R1 - Req S2 = DCOS(THETA) - DCOS(THETAeq) S3 = R2 - Req Y1 = (S1 + S3)/2.0d0 Y2 = S2 Y3 = (S1 - S3)/2.0d0 DO 88 I5=1,NV 88 FT(I5)=0.0D0 FT(1)=1.0d0 FT(2)=Y1 FT(3)=Y2 FT(4)=Y1**2 FT(5)=Y2**2 FT(6)=Y3**2 FT(7)=Y1*Y2 FT(8)=Y1**3 FT(9)=Y2**3 FT(10)=Y1**2*Y2 FT(11)=Y2**2*Y1 FT(12)=Y3**2*Y1 FT(13)=Y3**2*Y2 FT(14)=Y1**4 FT(15)=Y2**4 FT(16)=Y3**4 FT(17)=Y1**3*Y2 FT(18)=Y2**3*Y1 FT(19)=Y1**2*Y2**2 FT(20)=Y1**2*Y3**2 FT(21)=Y2**2*Y3**2 FT(22)=Y3**2*Y1*Y2 FT(23)=Y1**5 FT(24)=Y2**5 FT(25)=Y1**4*Y2 FT(26)=Y2**4*Y1 FT(27)=Y3**4*Y1 FT(28)=Y3**4*Y2 FT(29)=Y1**3*Y2**2 FT(30)=Y1**3*Y3**2 FT(31)=Y2**3*Y1**2 FT(32)=Y2**3*Y3**2 FT(33)=Y1**2*Y2*Y3**2 FT(34)=Y1*Y2**2*Y3**2 FT(35)=Y1**6 FT(36)=Y2**6 FT(37)=Y3**6 FT(38)=Y1**5*Y2 FT(39)=Y2**5*Y1 FT(40)=Y1**4*Y2**2 FT(41)=Y2**4*Y1**2 FT(42)=Y2**4*Y3**2 FT(43)=Y3**4*Y2**2 FT(44)=Y1**4*Y3**2 FT(45)=Y3**4*Y1**2 FT(46)=Y3**4*Y1*Y2 FT(47)=Y1**3*Y2**3 FT(48)=Y2**3*Y1**2*Y3**2 FT(49)=Y1**3*Y3**2*Y2 FT(50)=Y2**3*Y3**2*Y1 FT(51)=Y1**2*Y2**2*Y1**2 FT(52)=Y1**7 FT(53)=Y2**7 FT(54)=Y1**6*Y2 FT(55)=Y2**6*Y1 FT(56)=Y3**6*Y1 FT(57)=Y3**6*Y2 FT(58)=Y1**5*Y2**2 FT(59)=Y1**5*Y3**2 FT(60)=Y2**5*Y1**2 FT(61)=Y2**5*Y3**2 FT(62)=Y1**4*Y2*Y3**2 FT(63)=Y1**4*Y2**3 FT(64)=Y2**4*Y1*Y3**2 FT(65)=Y2**4*Y1**3 FT(66)=Y3**4*Y1*Y2**2 FT(67)=Y3**4*Y2*Y1**2 FT(68)=Y3**4*Y1**3 FT(69)=Y3**4*Y2**3 V=ZERO DO 40 I=1,NV 40 V=V+CV(I)*FT(I) V=C0+SCALE*V RETURN END SUBROUTINE BODC18(V,R1,R2,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FT(117), cv(117) DATA C0/0.0D0/,SCALE/1.0D-6/ data nv/117/ DATA ZERO/0.0D0/ DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/ DATA cv/2520.8349350538,-91.1209324309,49.9198880161, *201.5986877584,-27.0761946057,213.0632778654,-42.3461775982, *-89.5081758955,3.1590606555,-82.9610483797,91.2462747271, *-345.6845385214,2.9125633339,-9.2743637228,19.1249589270, *126.8439173732,102.4748093797,29.1206937761,-0.4903235370, *548.5647508851,57.4425043137,-65.3115364766,-32.7831814238, *111.9709105628,-278.4048997815,-231.4930911427,-323.4608284981, *-98.0602324965,181.1676886771,-282.2900662875,167.5274873935, *47.8266275957,-213.7402930366,20.8555467237,169.6474418937, *-62.1263131255,71.9040085206,72.1511749317,123.9432784777, *173.6777915640,58.5107154936,-27.8257872910,99.6002434065, *-135.2045098416,145.2475743775,-23.3633794048,-599.5375919902, *-540.4989859645,275.5466454040,-502.8009780288,167.2307803194, *-109.3097362364,-77.1423552840,492.0803282409,313.6574722992, *-22.4885664677,141.5175472800,-1134.9499411146,88.6368503668, *-587.9451428738,-247.5329792460,-4.9333140627,322.2004639448, *916.1768044740,479.0204491692,93.9437933859,107.5055433254, *0.0,63.1640462956,191.6712842623,1.6032078252,78.8696064123, *1.0174252672,-377.6585243760,-192.1212462644,653.6712620983, *0.0,261.5131272695,146.4973633837,0.0,-137.1867814770, *-98.7015480745,0.0,0.0,0.0,0.0,-97.1438825386,0.0,0.0,0.0,0.0, *0.0,22.7027546667,0.0,0.0,0.0,-341.9933632553,0.0,0.0, *-111.8520186501,0.0,0.0,0.0,0.0,-48.2716070600,0.0,0.0,0.0,0.0, *-5.0447226741,0.0,0.0,0.0,0.0,0.0,0.0,144.1528855579/ THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/toang S1 = R1 - Req S2 = DCOS(THETA) - DCOS(THETAeq) S3 = R2 - Req Y1 = (S1 + S3)/2.0d0 Y2 = S2 Y3 = (S1 - S3)/2.0d0 c Potential in Y1,Y2,Y3. Full set of c derivatives for each order (1-4) and upper terms for others. DO 88 I5=1,NV 88 FT(I5)=0.0D0 FT(1)=1.0d0 FT(2)=Y1 FT(3)=Y2 FT(4)=Y1**2 FT(5)=Y2**2 FT(6)=Y3**2 FT(7)=Y1*Y2 FT(8)=Y1**3 FT(9)=Y2**3 FT(10)=Y1**2*Y2 FT(11)=Y2**2*Y1 FT(12)=Y3**2*Y1 FT(13)=Y3**2*Y2 FT(14)=Y1**4 FT(15)=Y2**4 FT(16)=Y3**4 FT(17)=Y1**3*Y2 FT(18)=Y2**3*Y1 FT(19)=Y1**2*Y2**2 FT(20)=Y1**2*Y3**2 FT(21)=Y2**2*Y3**2 FT(22)=Y3**2*Y1*Y2 FT(23)=Y1**5 FT(24)=Y2**5 FT(25)=Y1**4*Y2 FT(26)=Y2**4*Y1 FT(27)=Y3**4*Y1 FT(28)=Y3**4*Y2 FT(29)=Y1**3*Y2**2 FT(30)=Y1**3*Y3**2 FT(31)=Y2**3*Y1**2 FT(32)=Y2**3*Y3**2 FT(33)=Y1**2*Y2*Y3**2 FT(34)=Y1*Y2**2*Y3**2 FT(35)=Y1**6 FT(36)=Y2**6 FT(37)=Y3**6 FT(38)=Y1**5*Y2 FT(39)=Y2**5*Y1 FT(40)=Y1**4*Y2**2 FT(41)=Y2**4*Y1**2 FT(42)=Y2**4*Y3**2 FT(43)=Y3**4*Y2**2 FT(44)=Y1**4*Y3**2 FT(45)=Y3**4*Y1**2 FT(46)=Y3**4*Y1*Y2 FT(47)=Y1**3*Y2**3 FT(48)=Y2**3*Y1**2*Y3**2 FT(49)=Y1**3*Y3**2*Y2 FT(50)=Y2**3*Y3**2*Y1 FT(51)=Y1**2*Y2**2*Y1**2 FT(52)=Y1**7 FT(53)=Y2**7 FT(54)=Y1**6*Y2 FT(55)=Y2**6*Y1 FT(56)=Y3**6*Y1 FT(57)=Y3**6*Y2 FT(58)=Y1**5*Y2**2 FT(59)=Y1**5*Y3**2 FT(60)=Y2**5*Y1**2 FT(61)=Y2**5*Y3**2 FT(62)=Y1**4*Y2*Y3**2 FT(63)=Y1**4*Y2**3 FT(64)=Y2**4*Y1*Y3**2 FT(65)=Y2**4*Y1**3 FT(66)=Y3**4*Y1*Y2**2 FT(67)=Y3**4*Y2*Y1**2 FT(68)=Y3**4*Y1**3 FT(69)=Y3**4*Y2**3 FT(70)=Y1**3*Y3**2*Y2**2 FT(71)=Y1**8 FT(72)=Y2**8 FT(73)=Y3**8 FT(74)=Y1**7*Y2 FT(75)=Y2**7*Y1 FT(76)=Y1**6*Y2**2 FT(77)=Y1**6*Y3**2 FT(78)=Y2**6*Y1**2 FT(79)=Y2**6*Y3**2 FT(80)=Y3**6*Y1**2 FT(81)=Y3**6*Y2**2 FT(82)=Y3**6*Y2*Y1 FT(83)=Y1**5*Y3**2*Y2 FT(84)=Y1**5*Y2**3 FT(85)=Y2**5*Y1**3 FT(86)=Y2**5*Y3**2*Y1 FT(87)=Y1**4*Y2**4 FT(88)=Y1**4*Y3**4 FT(89)=Y1**4*Y2**2*Y3**2 FT(90)=Y2**4*Y3**4 FT(91)=Y2**4*Y3**2*Y1**2 FT(92)=Y3**4*Y1**2*Y2**2 FT(93)=Y3**4*Y1*Y2**3 FT(94)=Y3**4*Y1**3*Y2 FT(95)=Y1**3*Y2**3*Y3**2 FT(96)=Y2**9 FT(97)=Y2**9*Y1*Y3**4 FT(98)=Y2**9*Y1**5 FT(99)=Y2**10 FT(100)=Y2**10*Y1**2*Y3**2 FT(101)=Y2**10*Y3**4 FT(102)=Y2**10*Y1**4 FT(103)=Y2**11 FT(104)=Y2**11*Y1*Y3**2 FT(105)=Y2**11*Y1**3 FT(106)=Y2**12 FT(107)=Y2**12*Y1**2 FT(108)=Y2**12*Y3**2 FT(109)=Y2**13 FT(110)=Y2**13*Y1 FT(111)=Y2**14 FT(112)=Y2**8*Y1**6 FT(113)=Y2**8*Y3**6 FT(114)=Y2**8*Y1**4*Y3**2 FT(115)=Y2**8*Y3**4*Y1**2 FT(116)=Y2**7*Y1**7 FT(117)=Y2**7*Y1**5*Y3**2 V=ZERO DO 40 I=1,NV 40 V=V+CV(I)*FT(I) C SCALE AND SHIFT THE ZERO V=C0+SCALE*V RETURN END