c
c------------------------------------------------------gauss_qn---
c
subroutine gauss_qn (m,ep,om)
c
c Evaluations of Gauss points using normal quadrature.
c
implicit none
c
include '_param.for'
c
integer m,k,i,ii,ir1,ir0,expon
real*8 coef(50),root(0:50,2),toler,epsilon(50),omeg(50)
real*8 ep(50),om(50)
real*8 n
c
n=m
c
toler=1.0d0/10.0d0**14
if(n.gt.50.or.n.lt.2) then
WRITE(F_W,3000)
stop 1111
end if
c
c Initiations.
c
ir1=1
ir0=2
root(0,1)=-1.0d0
root(1,1)= 0.0d0
root(2,1)= 1.0d0
c
c Computation of polynomial coefficients and roots.
c
do 1 k=2,n
do 2 i=0,n/2
call coefic(k,i,coef(i+1))
2 continue
do 3 i=1,k
call roots(k,coef,root(i-1,ir1),root(i,ir1),
@ toler,root(i,ir0))
3 continue
root(k+1,ir0)=1.0d0
root(0,ir0)=-1.0d0
ii=ir0
ir0=ir1
ir1=ii
1 continue
c
c Computation of weights.
c
do 4 i=1,n
epsilon(i)=root(i,ir1)
omeg(i)=0.0d0
do 5 k=0,n/2
expon=n-2*k
if(expon.gt.1)then
omeg(i)=omeg(i)+coef(k+1)*expon*epsilon(i)**(expon-1)
else
omeg(i)=omeg(i)+coef(k+1)*expon
end if
5 continue
omeg(i)=omeg(i)/2.0d0**n
omeg(i)=2.0d0/(omeg(i)**2.0d0*(1.0d0-epsilon(i)**2.0d0))
ep(i)=epsilon(i)
om(i)=omeg(i)
4 continue
return
c
c Printing formats.
c
3000 format(5x,'ERROR: ',1x,
@ 'The number of Gauss points used is not allowed')
end
c
c-------------------------------------------------------coefic--
c
subroutine coefic(n,k,value)
integer inter,k,n,n1,n2,n3,i
real*8 sig,value
value=1.0d0
inter=k/2
sig=1.0d0
if(inter*2-k.ne.0) sig=-1.0d0
n2=n-k
n1=2*n2
n3=n-2*k
if(n1.eq.0)then
value=1.0d0
else
do 1 i=n2+1,n1
value=value*REAL(i,8)
1 continue
end if
if(n3.ne.0) then
do 2 i=1,n3
value=value/REAL(i,8)
2 continue
end if
if(k.ne.0) then
do 3 i=1,k
value=value/REAL(i,8)
3 continue
end if
value=value*sig
return
end
c
c-------------------------------------------------------------evalua--
c
subroutine evalua(ncoef,coef,x,value)
implicit none
integer ncoef,k,expon
real*8 coef(ncoef),x,value
value=0
do 1 k=0,ncoef/2
expon=ncoef-2*k
if(expon.eq.0) then
value=value+coef(k+1)
else
value=value+coef(k+1)*x**REAL(expon,8)
end if
1 continue
return
end
c
c-------------------------------------------------------------roots---
c
subroutine roots(ncoef,coef,x00,xff,error,xi)
implicit none
integer ncoef
real*8 x00,xff,x0,xf,xi,error,coef(ncoef),value0,valuef,valuei,
@ errac
x0=x00
xf=xff
call evalua(ncoef,coef,x0,value0)
call evalua(ncoef,coef,xf,valuef)
errac=dabs(xf-x0)
do while (errac.gt.error)
xi=(xf+x0)/2.0d0
call evalua(ncoef,coef,xi,valuei)
if(dabs(valuei).lt.error) return
if(valuei*value0.lt.0.0d0) then
xf=xi
valuef=valuei
else
x0=xi
value0=valuei
end if
errac=DABS(xf-x0)
end do
return
end
c
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
