-
Notifications
You must be signed in to change notification settings - Fork 9
/
vcalc.f
65 lines (53 loc) · 1.54 KB
/
vcalc.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
subroutine vcalc(r,theta,chi,r0,rmax,rstep,nrgrd,nthgrd,nchgrd,
+ vtable,vpes)
implicit double precision(a-h,o-z)
dimension vtable(0:nrgrd*nthgrd*nchgrd-1)
parameter(wno2k=0.6950356d0,wno2kj=1.196266d-2)
c parameter(maxrpt=nrgrd-1,mxthpt=180,mxchpt=180)
c character*30 argum
c ... the subroutine takes r / bohr, theta / deg, and chi / deg
maxrpt=nrgrd-1
mxthpt=nthgrd-1
mxchpt=nchgrd-1
if(r.lt.r0) r=r0
if(r.gt.rmax)r=rmax
ir=int((r-r0)/rstep)
ith=int(theta)
ich=int(chi)
if(ich.gt.mxchpt)ich=mxchpt
if(ich.lt.0)ich=0
if(ith.gt.mxthpt)ith=mxthpt
if(ith.lt.0)ith=0
c write(6,*)ir,ith,ich
indrtc=(ir*nthgrd+ith)*nchgrd+ich
v0=vtable(indrtc)
c write(6,*)v0,v0*wno2k
c ... linear interpolation
if(ir.eq.maxrpt) then
gradr=0.d0
delr=0.d0
else
indrtc=((ir+1)*nthgrd+ith)*nchgrd+ich
gradr=(vtable(indrtc)-v0)/rstep
delr=r-(r0+ir*rstep)
endif
if(ith.eq.mxthpt) then
gradth=0.d0
delth=0.d0
else
indrtc=(ir*nthgrd+(ith+1))*nchgrd+ich
gradth=(vtable(indrtc)-v0)
delth=theta-dfloat(ith)
endif
if(ich.eq.mxchpt) then
gradch=0.d0
delch=0.d0
else
indrtc=(ir*nthgrd+ith)*nchgrd+ich+1
gradch=(vtable(indrtc)-v0)
delch=chi-dfloat(ich)
endif
vpes=v0+gradr*delr+gradth*delth+gradch*delch
c write(6,*)vpes,vpes*wno2k,vpes*wno2k*wno2kj
return
end