c-------------------------------------------------------
      subroutine bwr(jr1,in,ik,beta)

      implicit double precision (a-h,o-z)
       double precision beta
c      real*8 x,y,z,dx,dy,dz,cphi,sphi,beta
c      real*8 eps,eps1,one,pi,pi2,rxy,perpx,perpy,perpz
c      real*8 rp,rxyp,rxyp2,seta,sm,sm1,a1,a2,t1,t2
c example  call bwr(n,2,n,beta)
      parameter(maxa=500,pi=3.1415926535d0,pi2=2.d0*pi,pi12=pi/2.d0,
     *  eps=1.d-8, eps1=1.d-9, one=1.d0-1.d-12)

      common /matrix/ x(maxa),y(maxa),z(maxa),bangle(maxa)
      common /ddx/ dx(maxa),dy(maxa),dz(maxa)
      common /wr/jwr,www(maxa),wri(maxa,maxa)
      dimension cphi(maxa),sphi(maxa),ibad(maxa)

      beta=0.0
      imax=ik-in+3
      ibas=in-2

      do 11 ii=1,imax
        i=ibas+ii
        if(i.eq.0) i=jr1
        if(i.gt.jr1) i=i-jr1
        rxy=dsqrt(dx(i)*dx(i)+dy(i)*dy(i))
        if(rxy.ge.eps) then
          cphi(i)=dx(i)/rxy
          sphi(i)=dy(i)/rxy
          ibad(i)=0
        else
          cphi(i)=1.
          sphi(i)=0.
          ibad(i)=1
        endif
   11 continue

      mmax=imax-1

      do 1 mm=1,mmax
        m=ibas+mm
        m1=m+1
        if(m.eq.0) m=jr1
        if(m.gt.jr1) m=m-jr1
        if(m1.gt.jr1) m1=m1-jr1
        if(ibad(m).eq.1) then
          perpx=            -dz(m)*dy(m1)
          perpy=dz(m)*dx(m1)-eps  *dz(m1)
          perpz=eps  *dy(m1)
        elseif(ibad(m1).eq.1) then
          perpx=dy(m)*dz(m1)
          perpy=dz(m)*eps   -dx(m)*dz(m1)
          perpz=            -dy(m)*eps
        else
          perpx=dy(m)*dz(m1)-dz(m)*dy(m1)
          perpy=dz(m)*dx(m1)-dx(m)*dz(m1)
          perpz=dx(m)*dy(m1)-dy(m)*dx(m1)
        endif

        if(dabs(perpz).lt.eps1) then
         www(m1)=0.d0
          goto 1
        endif
        perpzr=perpz
        si=sign(1.d0,perpzr)

        rxyp2=perpx*perpx+perpy*perpy
        rxyp=sqrt(rxyp2)
        if(rxyp.lt.eps1) then
         www(m1)=0.d0
          goto 1
        endif
        rp=sqrt(rxyp2+perpz*perpz)
        seta=rxyp/rp
        cchi=perpx/rxyp*si
        schi=perpy/rxyp*si

        sm1=sphi(m1)*cchi - cphi(m1)*schi
        sm =sphi(m) *cchi - cphi(m) *schi

        a1=seta*sm1
        if(a1.gt.one) then
          t1=pi12
        elseif(a1.lt.-one) then
          t1=-pi12
        else
          t1=asin(a1)
        endif
        a2=seta*sm
        if(a2.gt.one) then
          t2=pi12
        elseif(a2.lt.-one) then
          t2=-pi12
        else
          t2=asin(a2)
        endif
       www(m1)=t1-t2
        beta=beta+www(m1)
    1 continue

      beta=beta/pi2

      return
      end