cmBremE2.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cmbreme2 (up, v, nc)

Function/Subroutine Documentation

cmbreme2()

subroutine cmbreme2	(	real*8	up,
		real*8	v,
		integer	nc
	)

Definition at line 38 of file cmBremE2.f.

References d, d0, false, ksamppw(), ksamppwmrg(), ksamppwx(), parameter(), rndc(), and true.

Referenced by cmbreme().

      real*8 up       ! input Upsilon value
      real*8 v        ! output. sampled fractional value
          integer nc      ! output. number of rejections. average trial
!                                   number is 1.16 for up=0.01 and 1.26 for
!                                   up=10.
!
      integer n1, n2, n, i, na
      parameter(n1 = 3, n2 = 5, n= 12)
      real*8 c1(n1), x10(n1), p1(n1), x1p(n1+1), coef1(n1)  ! for 1/(2+3Uz)^3
!
      real*8 c2(n2), x20(n2), p2(n2), x2p(n2+1), coef2(n2) ! for z^2K2/3(z)
!
      real*8 coef(n), p(n), xp(n+1)
      real*8 z, fz, ft, ck23, u
      logical first

      data p1/0.04433d0, 1.5d0, 2.8125d0/
      data c1/0.1195d0,  1.5625d-2, 3.05175d-5/

!
      data x20/0.01d0, 1.d0, 2.511886d0, 7.4131d0,
     *        19.95d0/
      data c2/2.3096d-3, 4.9463d-1, 0.41808d0,
     *        1.5448d-2, 2.427d-7/
      data p2/-1.3298d0, -0.446d0, 1.04103, 5.92368d0,
     *        18.007d0/
     

      data first/.true./
      save first, x1p, coef1, x2p, coef2, x10, x20
      
!

      x10(1) = 0.01d0/up
      x10(2) = 2.d0/3.d0/up
      x10(3) = 10.d0/up
!

      call ksamppwx(c1, p1, x10, n1, x1p(2), coef1 )
      x1p(1) = 0.
      x1p(n1+1) =  30./up
!          write(*,*)' x1p, coef1=', x1p, coef1
!   
      if(first) then
         call ksamppwx(c2, p2, x20, n2, x2p(2), coef2 )
         first = .false.
         x2p(1) = 0.
         x2p(n2+1) = 30.
!             write(*,*) ' x2p, coef2=', x2p, coef2
      endif 
    call ksamppwmrg(coef1, p1, x1p, n1,
     *  coef2, p2, x2p, n2,
     *  coef,  p,  xp,  na)
!       write(*,*)'n, coef, p, xp=',  na, coef, p, xp
    do i =1, 1000
       call ksamppw(i, coef, p, xp, na, z, fz)
!           write(*,*) ' z, fz=', z, fz
       ft = z*z * ck23(z) / (2.+ 3.* up*z)**3
       call rndc(u)
       if(u .lt. ft/fz) goto 10
    enddo
 10 continue    
    nc = i
    v =  3.*up*z/(2.+ 3.*up*z)

Here is the call graph for this function:

Here is the caller graph for this function: