Line data Source code
1 46 : subroutine vinta(x,a,y,nn,ia,id)
2 : c vinta
3 : c sets integral of the real variable a(1,n) ( =a(1,x(n)) into
4 : c real y(1,n), n=1,nn
5 : c
6 : c
7 : c the independent variable x, which need not be uniformly divided
8 : c or increasing with n (but must be monotonic), must be supplied
9 : c by the calling programme
10 : c
11 : c
12 : c
13 : c
14 : c Double precision version.
15 : c ++++++++++++++++++++++++
16 : c
17 : c Dated 10/3/90.
18 : c
19 : c Note: this double precision version of the routine has same name
20 : c as single precision version, but is distinguished by its file name
21 : c
22 : implicit double precision(a-h,o-z)
23 : dimension a(nn),y(nn),x(nn)
24 : c
25 : c common defining standard input and output
26 : c
27 : common/cstdio/ istdin, istdou, istdpr, istder
28 : c
29 46 : if(nn.eq.1) then
30 0 : if(istdpr.gt.0) write(istdpr,100)
31 0 : y(1)=0
32 0 : return
33 46 : else if(nn.eq.2) then
34 : c
35 : c use trapezoidal integration
36 : c
37 0 : y(1)=0
38 0 : y(1+id)=0.5d0*(x(2)-x(1))*(a(1)+a(1+ia))
39 0 : return
40 : end if
41 : c
42 46 : 1 wc=a(1)
43 46 : ir=1
44 46 : 3 is=ir*ia
45 46 : ird=ir*id
46 46 : n=nn
47 46 : m=n/2
48 46 : ie=n-m*2
49 46 : if(ie.eq.0) m=m-1
50 46 : w=1.0d0
51 46 : r=x(n)-x(1)
52 46 : if(r.eq.0.d0) go to 50
53 : c
54 46 : j3=1
55 46 : k3=1
56 46 : y(1)=0.d0
57 46 : ah2=0.d0
58 46 : vsimf=0.d0
59 : c
60 33764 : 4 do 7 i=1,m
61 33718 : i2=i*2
62 33718 : i1=i2-1
63 33718 : i3=i2+1
64 33718 : hn=(x(i2)-x(i1))/6.d0
65 33718 : hn1=(x(i3)-x(i2))/6.d0
66 33718 : if(r*hn.le.0.d0.or.r*hn1.le.0.d0) go to 51
67 33718 : wa=hn+hn1
68 33718 : wb=hn-hn1
69 33718 : wd=wa/hn
70 33718 : ah=wd*(hn+wb)
71 33718 : ah2=wa/hn1
72 33718 : ah1=wa*wd*ah2
73 33718 : ah2=ah2*(hn1-wb)
74 33718 : bh1=hn+3.d0*hn1
75 33718 : bh2=hn/hn1
76 33718 : bh=(bh1+hn)*hn/wa
77 33718 : bh1=bh1*bh2
78 33718 : bh2=-hn*hn*bh2/wa
79 33718 : wa=wc
80 33718 : k2=k3+is
81 33718 : k3=k2+is
82 33718 : wb=a(k2)
83 33718 : wc=a(k3)
84 33718 : 6 j2=j3+ird
85 33718 : j3=j2+ird
86 33718 : y(j2)=vsimf+bh*wa+bh1*wb+bh2*wc
87 33718 : vsimf=vsimf+ah*wa+ah1*wb+ah2*wc
88 33718 : y(j3)=vsimf
89 46 : 7 continue
90 : c
91 46 : if(ie.ne.0) return
92 : c
93 46 : wa=x(n-1)
94 46 : hn=(wa-x(n-2))/6.d0
95 46 : hn1=(x(n)-wa)/6.d0
96 46 : if(r*hn.le.0.d0.or.r*hn1.le.0.d0) go to 52
97 46 : wa=hn+hn1
98 46 : wd=hn1/hn
99 46 : ah=-hn1*hn1*wd/wa
100 46 : vsimf=vsimf+ah*wb
101 46 : wb=3.0d0*hn+hn1
102 46 : ah1=wd*wb
103 46 : ah2=hn1*(wb+hn1)/wa
104 46 : wa=wc
105 46 : wb=a(k3+is)
106 46 : 9 vsimf=vsimf+ah1*wa+ah2*wb
107 46 : y(j3+ird)=vsimf
108 : c
109 46 : return
110 : c
111 : c
112 : c diagnostics
113 0 : 50 if(istdpr.gt.0) write(istdpr,60) n,x(1)
114 0 : return
115 0 : 51 i=i2+1
116 0 : if(istdpr.gt.0)
117 0 : * write(istdpr,61) i1,x(i1),i2,x(i2),i,x(i),n,x(1),x(n)
118 0 : return
119 0 : 52 if(istdpr.gt.0)
120 0 : * write(istdpr,61) i2,x(i2),i1,x(i1),n,x(n),n,x(1),x(n)
121 0 : return
122 : c
123 : c
124 : c
125 : 60 format(//1x,10('*'),5x,'null range of independent variable in ',
126 : * 'intf/inta',5x,10('*')//21x,'x(1) = x(',i4,') =',1pe14.6/)
127 : 61 format(//1x,10('*'),5x,'independent variable not monotonic in',
128 : * ' intf/inta',5x,10('*')//16x,'x(',i4,') =',1pe13.5,
129 : * ', x(',i4,') =',1pe13.5,', x(',i4,') =',1pe13.5/
130 : * 16x,'number of mesh points =',i4, 3x,'x(1) =',1pe13.5,
131 : * 8x,'x(n) =',1pe13.5/)
132 : c
133 : 100 format(/' **** Error in vinta: nn = 1')
134 : end
|
