source: PSPA/parmelaPSPA/trunk/gaus.f @ 315

      subroutine gaus(r1,r2,z1,z2,opt,er,ez)
c---calculate er and ez at the r and z location given in fldcom
c---by gauss quadrature integration over the double
c---interval from r1 to r2 and from z1 to z2.
c---if opt.gt.0, determine number of integration points as follows.
c---let rat = max(cr/cz,cz/cr), where cr=r2-r1, and cz=z2-z1
c---     if rat.le.2, use 2 x 2 point array
c---     if rat.gt.2, use 2 x 4 point array
c---     if rat.gt.4, use 2 x 6 point array
c---if opt.eq.0, use 2 x 2 point array regardless of rat
c---if opt.lt.0, use one point.
c--- rat for ratio
c--------------------------------------------------------------------------
      save
c
      include 'ucom.h'
c
      dimension r(6),z(6),wr(6),wz(6),xx(3,3),wx(3,3)
c--------------------------------------------------------------------------
c*
      data ((xx(i,j),i=1,3),j=1,3)/.2113248654d0,0.,0.,.0694318442d0,
     1  .3300094782d0,0.,.0337652429d0,.1693953068d0,.380690407d0/
      data ((wx(i,j),i=1,3),j=1,3)/.5,0.,0.,.1739274226d0,.3260725774d0,
     1  0., .08566224619d0,.180380786500d0,.233956967300d0/
      cr=r2-r1
      cz=z2-z1
      rfac=2./((r2**2-r1**2)*(z2-z1))
      crczrfac=cr*cz*rfac
      ir=1
      jz=1
      m=1
      if (opt.eq.0.) go to 20
      if (opt.lt.0.) go to 60
c---determine number of integration points
      rat=abs(cz/cr)
      l=0
      if (rat.ge.1.)go to 10
      rat=1./rat
      l=1
   10 continue
      if (rat.gt.2.)m=2
      if (rat.gt.4.)m=3
      if (l.eq.0) jz=m
      if (l.eq.1) ir=m
cdir$ IVDEP
   20 continue
      do 30 i=1,ir
      k=2*i-1
      r(k)=r1+cr*xx(i,ir)
      r(k+1)=r2-cr*xx(i,ir)
      wr(k)=wx(i,ir)
      wr(k+1)=wx(i,ir)
   30 continue
cdir$ IVDEP
      do 40 j=1,jz
      k=2*j-1
      z(k)=z1+cz*xx(j,jz)
      z(k+1)=z2-cz*xx(j,jz)
      wz(k)=wx(j,jz)
      wz(k+1)=wx(j,jz)
   40 continue
      ser=0.
      sez=0.
      kr=2*ir
      kz=2*jz
      do 50 i=1,kr
      do 50 j=1,kz
      call flds(r(i),z(j),er1,ez1)
      ser=ser+wr(i)*wz(j)*er1*r(i)
      sez=sez+wr(i)*wz(j)*ez1*r(i)
   50 continue
      er=crczrfac*ser
      ez=crczrfac*sez
      return
   60 continue
      rs=sqrt(.5*(r1**2+r2**2))
      zs=.5*(z1+z2)
      call flds (rs,zs,er,ez)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*