1!$Id:$
2      subroutine int1dg(l,sw)
3
4!      * * F E A P * * A Finite Element Analysis Program
5
6!....  Copyright (c) 1984-2017: Regents of the University of California
7!                               All rights reserved
8
9!-----[--.----+----.----+----.-----------------------------------------]
10!      Purpose: Gauss quadrature for 1-d element
11
12!      Inputs:
13!         l     - Number of points
14
15!      Outputs:
16!         sw(1,*) - Gauss points
17!         sw(2,*) - Gauss weights
18!-----[--.----+----.----+----.-----------------------------------------]
19      implicit  none
20
21      integer   l
22      real*8    sw(2,*), t
23
24      save
25
26      if(l.eq.1) then
27
28        sw(1,1) = 0.0d0
29        sw(2,1) = 2.0d0
30
31      elseif(l.eq.2) then
32
33        sw(1,1) = -1.d0/sqrt(3.d0)
34        sw(1,2) = -sw(1,1)
35        sw(2,1) = 1.0d0
36        sw(2,2) = 1.0d0
37
38      elseif(l.eq.3) then
39
40        sw(1,1) = -sqrt(0.6d0)
41        sw(1,2) = 0.0d0
42        sw(1,3) = -sw(1,1)
43        sw(2,1) = 5.d0/9.d0
44        sw(2,2) = 8.d0/9.d0
45        sw(2,3) = sw(2,1)
46
47      elseif(l.eq.4) then
48
49        t       =  sqrt(4.8d0)
50        sw(1,1) = -sqrt((3.d0+t)/7.d0)
51        sw(1,2) = -sqrt((3.d0-t)/7.d0)
52        sw(1,3) = -sw(1,2)
53        sw(1,4) = -sw(1,1)
54        t       =  1.d0/3.d0/t
55        sw(2,1) =  0.5d0 - t
56        sw(2,2) =  0.5d0 + t
57        sw(2,3) =  sw(2,2)
58        sw(2,4) =  sw(2,1)
59
60      elseif(l.eq.5) then
61
62        t       =  sqrt(1120.0d0)
63
64        sw(1,1) = (70.d0+t)/126.d0
65        sw(1,2) = (70.d0-t)/126.d0
66
67        t       =  1.d0/(15.d0 * (sw(1,2) - sw(1,1)))
68
69        sw(2,1) = (5.0d0*sw(1,2) - 3.0d0)*t/sw(1,1)
70        sw(2,2) = (3.0d0 - 5.0d0*sw(1,1))*t/sw(1,2)
71        sw(2,3) =  2.0d0*(1.d0 - sw(2,1) - sw(2,2))
72        sw(2,4) =  sw(2,2)
73        sw(2,5) =  sw(2,1)
74
75        sw(1,1) = -sqrt(sw(1,1))
76        sw(1,2) = -sqrt(sw(1,2))
77        sw(1,3) =  0.0d0
78        sw(1,4) = -sw(1,2)
79        sw(1,5) = -sw(1,1)
80
81!     Compute points and weights
82
83      else
84
85        call gausspw(l,sw)
86
87      endif
88
89      end
90
91      subroutine gausspw (nn,sw)
92
93!-----[--.----+-!--.----+----.----+------------------------------------]
94!     Input:
95!     nn       = Number of Gauss Points to compute
96
97!     Outputs:
98!     sw(1,*)  = Gauss Point Coordinates
99!     sw(2,*)  = Gauss Point Weights
100!-----[--.----+-!--.----+----.----+------------------------------------]
101      implicit   none
102
103      integer    nn, n
104      real*8     sw(2,*)
105      real*8     fn, beta, cc, xt, dpn,pn1, flgama
106
107      fn   = dble(nn)
108      beta = exp(2.d0*flgama(1.d0) - flgama(2.d0))
109      cc   = 2.d0*beta
110      do n = 2,nn
111        cc = cc*4.d0*dble(n-1)**4
112     &       / (dble(2*n-1)*dble(2*n-3)*dble(2*n-2)**2)
113      end do ! n
114
115      do n = 1,nn
116
117!       Largest zero
118
119        if(n.eq.1) then
120          xt = 1.d0 - 2.78d0/(4.d0 + fn*fn)
121
122!       Second zero
123
124        elseif(n.eq.2) then
125          xt = xt - (4.1d0 + 0.246d0*(fn - 8.d0)/fn)*(1.d0 - xt)
126
127!       Third zero
128
129        elseif(n.eq.3) then
130          xt = xt - (1.67d0 + 0.3674d0*(fn - 8.d0)/fn)*(sw(1,1) - xt)
131
132!       Second last zero
133
134        elseif(n.eq.nn-1) then
135          xt = xt + (xt - sw(1,n-2))/0.766d0/(1.d0 + 0.639d0*(fn-4.d0)
136     &                                      /(1.d0 + 0.710d0*(fn-4.d0)))
137
138!       Last zero
139
140        elseif(n.eq.nn) then
141          xt = xt + (xt - sw(1,n-2))/1.67d0/(1.d0 + 0.22d0*(fn-8.d0)/fn)
142
143!       Intermediate roots
144
145        else
146          xt    = 3.d0*sw(1,n-1) - 3.d0*sw(1,n-2) + sw(1,n-3)
147        endif
148
149!       Find root using xt-value
150
151        call root (xt,nn,dpn,pn1)
152        sw(1,n) = xt
153        sw(2,n) = cc/(dpn*pn1)
154
155      end do ! n
156
157!     Reverse order of points
158
159      do n = 1, nn
160        sw(1,n) = - sw(1,n)
161      end do ! n
162
163      end
164
165      subroutine root (x,nn,dpn,pn1)
166
167!-----[--.----+-!--.----+----.----+------------------------------------]
168!      Improve approximate root x; in addition we also obtain
169!         dpn = derivative of p(n) at x
170!         pn1 = value of p(n-1) at x
171!-----[--.----+-!--.----+----.----+------------------------------------]
172      implicit   none
173
174      logical    notconv
175      integer    nn, iter
176      real*8     x,dpn,pn1, d,p,dp
177
178      real*8     eps
179      data       eps / 1.d-39 /
180
181      iter = 0
182
183      notconv = .true.
184      do while(notconv .and. iter.lt.50)
185        iter = iter + 1
186        call recur (p,dp,pn1,x,nn)
187        d = p/dp
188        x = x - d
189        if(abs(d).le.eps) then
190          notconv = .false.
191        endif
192      end do ! while
193      dpn = dp
194
195      end
196
197      subroutine recur (pn,dpn,pn1,x,nn)
198
199      implicit   none
200
201      integer    nn, n
202      real*8     pn,dpn,pn1,x, c
203      real*8     p1, p,dp, dp1, q,dq
204
205      p1  = 1.d0
206      p   = x
207      dp1 = 0.d0
208      dp  = 1.d0
209      do n = 2,nn
210        c   = 4.d0*dble(n-1)**4
211     &       / (dble(2*n-1)*dble(2*n-3)*dble(2*n-2)**2)
212        q   = x*p  - c*p1
213        dq  = x*dp - c*dp1 + p
214        p1  = p
215        p   = q
216        dp1 = dp
217        dp  = dq
218      end do ! n
219      pn  = p
220      dpn = dp
221      pn1 = p1
222
223      end
224
225      function flgama (w)
226
227      implicit none
228
229      integer  m, i
230      real*8   flgama, w, pi,x, p, fk, y, z,zz
231
232      pi =  acos(-1.d0)
233      x  =  w
234      fk = -1.d0
235
236!     w less eq 0.5
237
238      if (x .lt. 0.5d0) then
239        m = 1
240        p = pi/sin(x*pi)
241        x = 1.d0 - x
242      else
243        m = 0
244        p = 0.0d0
245      endif
246
247      do while(x + fk - 6.d0 .le.0.0d0)
248        fk = fk + 1.d0
249      end do ! while
250
251      z  = x + fk
252      zz = z*z
253
254      y  = (z - 0.5d0)*log(z) - z + 0.9189385332047d0 + (((((-4146.d0/zz
255     &        + 1820.d0)/zz - 1287.d0)/zz + 1716.d0)/zz - 6006d0)/zz
256     &        + 180180.d0)/z/2162160.d0
257
258      if(fk.gt.0.0d0) then
259        do i = 1,int(fk)
260          fk = fk - 1.d0
261          y  = y - log(x + fk)
262        end do ! i
263      endif
264
265      if(m.ne.0) then
266        if(p.le.0.0d0) then
267          write (*,2000) w
268          y = 0.d0
269        else
270          y = log(p) - y
271        endif
272      endif
273      flgama = y
274
2752000  format (2x,'gamma(',e11.4,') is negative')
276
277      end
278