program BISprogram
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     This program uses the BIS routine to calculate the fiber stress  c
c     concentration factors that are changed away from 1 when one or   c
c     more fibers are broken. A broken fiber has stress factor 0, and  c
c     in a system without broken bonds all fibers have stress factor 1.c
c     When fibers are broken, their (unbroken) neighbors become        c
c     subject to increased stress.                                     c
c                                                                      c
c     In order to compile the program you need to link to the two      c
c     subroutines from Numerical recipes ludcmp and lubksb in the      c
c     double precision version (files dludcmp.for and dlubks.for).     c
c                                                                      c
c     This program is a wrap around the BIS subroutine in order to     c
c     test the subroutine and illustrate how it functions. Try running c
c     the program with various fibers broken (indicated by the '1's in c
c     the input file 'breaks.dat'). The subroutine can be lifted out   c
c     and used in the program for statistics of fracture. If you do so,c
c     make sure that the 'statistics of fracture' program defines      c
c     fm, breaklist, and nbreaks appropriately.                        c
c                                                                      c
c                                    Elsebeth Schroder, April 11 2000  c
c23456789012345678901234567890123456789012345678901234567890123456789012
c
      implicit none
      integer maxfibers
      parameter(maxfibers=401)
      real*8 fm(-(maxfibers-1)/2:(maxfibers-1)/2),
     *       breaklist(-(maxfibers-1)/2:(maxfibers-1)/2)
      integer fiber, nbreaks
c
c     read for a list of fibers whether they are broken (1)
c     or unbroken (0), count the number of breaks.
c
      nbreaks=0
      open(1,file='breaks.dat',type='old')
      do fiber=-(maxfibers-1)/2,(maxfibers-1)/2
         read(1,*) breaklist(fiber)
         if (breaklist(fiber).eq.1) then
            nbreaks=nbreaks+1
         endif
      enddo
      close(1)     
c
c     calculate the fiber stress concentration factors fm
c
      call BISrout(fm,breaklist,nbreaks)
c
c     write the fiber stress concentration factors to file
c
      open(2,file='factors.dat',type='unknown')
      do fiber=-(maxfibers-1)/2,(maxfibers-1)/2
         write(2,*) fiber, fm(fiber)
      enddo
      close(2)
c      
      stop
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine BISrout(fm,breaklist,nbreaks)
      implicit none
      integer maxfibers, nbreaks
      parameter(maxfibers=401)
      real*8 fm(-(maxfibers-1)/2:(maxfibers-1)/2),
     *       breaklist(-(maxfibers-1)/2:(maxfibers-1)/2)
      integer fiber, fiber2, break, break2
      integer indx(1:nbreaks)
      real*8 b(1:nbreaks), gamma(1:nbreaks,1:nbreaks), d
c     
c     set up the system of linear equations to be solved
c
      break=0
      do fiber=-(maxfibers-1)/2,(maxfibers-1)/2
         if(breaklist(fiber).eq.1) then
            break=break+1
            b(break)=1d0
            break2=0
            do fiber2=-(maxfibers-1)/2,(maxfibers-1)/2
               if (breaklist(fiber2).eq.1) then
                  break2=break2+1
                  gamma(break,break2)=-1.d0/(4d0*(fiber-fiber2)**2-1)
               endif
            enddo
         endif
      enddo
c
c     use Numerical Recipies subroutines to solve the equations
c
      call ludcmp(gamma,nbreaks,nbreaks,indx,d)
      call lubksb(gamma,nbreaks,nbreaks,indx,b)
c
c     calculate the fiber stress concentration factors fm
c
      do fiber=-(maxfibers-1)/2,(maxfibers-1)/2 
         break2=0
         fm(fiber)=0d0
         do fiber2=-(maxfibers-1)/2,(maxfibers-1)/2 
            if(breaklist(fiber2).eq.1) then
               break2=break2+1
               fm(fiber)=fm(fiber)+b(break2)/(4d0*(fiber-fiber2)**2-1)
            endif
         enddo
         fm(fiber)=fm(fiber)+1d0
      enddo
c
      return
      end