!=======================================================================
! Program dancoff - calculate Dancoff factor for pebble bed HTR
! J.L. Kloosterman
! Delft University of Technology
! Reactor Institute Delft
! Mekelweg 15, Delft
! Email: J.L.Kloosterman@tudelft.nl
! Website: www.JanLeenKloosterman.nl
! Phone: +31 15 278 1191
!=======================================================================
      program dancoff
      parameter (pi=3.14159265)
      real intra,inter,infdan
!-----read the radius of the fuel grain (cm) (usually 0.025 cm)
      read(5,*) rg1
!-----read the number of fuel grains per pebble
      read(5,*) grn
!-----read the radius of the fuel zone of the pebble (cm) (usually 2.5 cm)
      read(5,*) rp1
!-----read the radius of the graphite moderator shell (cm) (usually 3 cm)
      read(5,*) rp2
!-----sigma of graphite is 0.4097 1/cm
      sigma = 0.4097
      rg2   = rp1 / (grn**(1./3.))
!-----calculate dancoff factor
      dan1 = intra(rg1,rg2,rp1,rp2,sigma)
      dan2 = inter(rg1,rg2,rp1,rp2,sigma)
      cinf = grain(rg1,rg2,sigma)
      write(6,1100) dan1,dan2,cinf,dan1+dan2
!=======================================================================
 1100 format(' Intra Dancoff factor ',f8.4,/, &
             ' Inter Dancoff factor ',f8.4,/, &
             ' Infinite med Dancoff ',f8.4,/, &
             ' Total Dancoff factor ', f8.4)
!=======================================================================
      end

      real function grain(r1,r2,sigma)
!-----calculates the inf medium fuel grain Dancoff factor
!-----Eq. A.12 in thesis Evert Bende
      call trans(r1,r2,sigma,tio,toi,too)
      grain = tio*toi/(1.0-too)
      return
      end

      real function intra(rg1,rg2,rp1,rp2,sigma)
!-----calculates the intra-pebble Dancoff factor
!-----Eq. A.24 in thesis Evert Bende
      call trans(rg1,rg2,sigma,tio,toi,too)
      star  = -alog(too)/chord(rg2)
      intra = grain(rg1,rg2,sigma)*(1.0-pf(rp1*star))
      return
      end

      real function inter(rg1,rg2,rp1,rp2,sigma)
!-----calculates the inter-pebble Dancoff factor
!-----Eq. A.26 and A.28 in thesis Evert Bende
      call trans(rg1,rg2,sigma,tio,toi,too)
      star  = -alog(too)/chord(rg2)
      prob  = pf(rp1*star)
!-----Eq. A.26 in thesis Evert Bende
      tii   = 1.0-(4./3.)*rp1*star*prob
      fac   = (1.0-tii)/(1.0-tii*tio*toi)
      inter = grain(rg1,rg2,sigma)*grain(rp1,rp2,sigma)*prob*fac
      return
      end

      subroutine trans(r1,r2,sigma,tio,toi,too)
!-----calculates transmission probabilities in a sphere
!-----Eqs. A.3, A.5, A.6, and A.7 in thesis Bende
      nint = 100
      r12  = r1**2
      r22  = r2**2
      dy   = r1/real(nint)
      tio  = 0.0
      do i=1,nint
         y   = (i-0.5)*dy
         y2  = y**2
         u   = sqrt(r22-y2)-sqrt(r12-y2)
         tio = tio + y*exp(-sigma*u)*dy
      enddo
      tio = 2.0*tio
      toi = tio/r22
      tio = tio/r12
      a   = sigma*sqrt(r22-r12)
      too = (1.-(1.+2.*a)*exp(-2.*a)) / (2.*r22*(sigma**2))
      return
      end

      real function chord(r)
!-----calculates mean chord length of a sphere
!-----under Eq. A11 thesis Bende
      chord = 4.0*r/3.0
      return
      end

      real function pf(rlam)
!-----calculates the first flight escape probability for a sphere
!-----Eq. A.9 in thesis Evert Bende
      twor = 2.*rlam
      pf   = .75*(1.-(1.-(1.+twor)*exp(-twor))/(twor*rlam))/rlam
      return
      end