*###[ ffcb2:
	subroutine ffcb2(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,piDpj,ier)
***#[*comment:***********************************************************
*									*
*	Calculate    1   /    d^n Q Q(mu) Q(nu)				*
*		  ------ | ------------------------			*
*		  i pi^2 / (Q^2-ma^2)((Q+p)^2-mb^2)			*
*						 p mu			*
*		    = B2p*p(mu)*p(nu) + B2d*delta    /p^2		*
*						 p nu			*
*									*
*	Input:	cb1	   complex	vector twopoint function	*
*		cb0	   complex	scalar twopoint function	*
*		ca0i(2)	   complex	scalar onepoint function with	*
*						m1,m2			*
*		xp	   complex	p.p in B&D metric		*
*		xma,2	   complex	m_1^2,m_2^2			*
*		piDpj(3,3) complex	dotproducts between s1,s2,p	*
*		ier	   integer	digits lost so far		*
*	Output:	cb2p	   complex	coefficient of p(mu)*p(nu)	*
*		cb2d	   complex	coefficient of delta()/p^2	*
*		ier	   integer	digits lost			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	DOUBLE COMPLEX xp,xma,xmb,piDpj(3,3)
	DOUBLE COMPLEX cb2p,cb2d,cb1,cb0,ca0i(2)
*
*	local variables
*
	integer ier0,i,j
	DOUBLE COMPLEX dmap,dmbp,dmamb,cc
	DOUBLE PRECISION rm1,rm2,rp,rpiDpj(3,3),sprec
*
*	common blocks
*
	include 'ff.h'
*
*  #] declarations:
*  #[ real case:
	if ( DIMAG(xma).eq.0 .and. DIMAG(xmb).eq.0 ) then
	    rm1 = DBLE(xma)
	    rm2 = DBLE(xmb)
	    rp  = DBLE(xp)
	    do 20 j=1,3
		do 10 i=1,3
		    rpiDpj(i,j) = DBLE(piDpj(i,j))
   10		continue
   20	    continue
	    sprec = precx
	    precx = precc
	    call ffxb2(cb2p,cb2d,cb1,cb0,ca0i,rp,rm1,rm2,rpiDpj,
     +	    	ier)
	    precx = sprec
	    return
	endif
*  #] real case:
*  #[ get differences:
	ier0 = 0
	dmamb = xma - xmb
	dmap = xma - xp
	dmbp = xmb - xp
*  #] get differences:
*  #[ call ffcb2a:
	call ffcb2a(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,dmap,dmbp,dmamb,
     +		piDpj,ier)
*  #] call ffcb2a:
*###] ffcb2:
	end
*###[ ffcb2a:
	subroutine ffcb2a(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,
     +		dmap,dmbp,dmamb,piDpj,ier)
***#[*comment:***********************************************************
*									*
*	see ffcb2, plus differences.					*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	DOUBLE COMPLEX xp,xma,xmb,dmap,dmbp,dmamb,piDpj(3,3)
	DOUBLE COMPLEX cb2p,cb2d,cb1,cb0,ca0i(2)
*
*	local variables
*
	integer i,j,ier0,init,ithres
	logical llogmm,lreal
	DOUBLE PRECISION absc,xmax,xmxp,rloss
	DOUBLE PRECISION rm1,rm2,rp,rm1p,rm2p,rm1m2,rpiDpj(3,3),sprec
	DOUBLE COMPLEX delsp,xlam,xlo3,xlogmm,zfflo1,zfflo3
	DOUBLE COMPLEX cc,cs(6),cb21,cb22,csom
	DOUBLE COMPLEX cqi(3),cqiqj(3,3),qiDqj(3,3)
	save init
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*
*  #] declarations:
*  #[ real case:
	if ( DIMAG(xma).eq.0 .and. DIMAG(xmb).eq.0 ) then
	    lreal = .TRUE.
	elseif ( nschem.le.4 ) then
	    lreal = .TRUE.
	    if ( init.eq.0 ) then
	    	init = 1
	    	print *,'ffcb2a: nschem <= 4, ignoring complex masses:',
     +	    		nschem
	    endif
	elseif ( nschem.le.6 ) then
	    if ( init.eq.0 ) then
	    	init = 1
	    	print *,'ffcb2a: nschem = 5,6 complex masses near ',
     +	    		'threshold: ',nschem
	    endif
	    cqi(1) = xma
	    cqi(2) = xmb
	    cqi(3) = xp
	    cqiqj(1,2) = dmamb
	    cqiqj(2,1) = -cqiqj(1,2)
	    cqiqj(1,3) = dmap
	    cqiqj(3,1) = -cqiqj(1,3)
	    cqiqj(2,3) = dmbp
	    cqiqj(3,2) = -cqiqj(2,3)
	    cqiqj(1,1) = 0
	    cqiqj(2,2) = 0
	    cqiqj(3,3) = 0
	    call ffthre(ithres,cqi,cqiqj,3,1,2,3)
	    if ( ithres.eq.0 .or. ithres.eq.1 .and. nschem.eq.5 ) then
		lreal = .TRUE.
	    else
		lreal = .FALSE.
	    endif
	else
	    lreal = .FALSE.
	endif
	if ( lreal ) then
	    rm1 = DBLE(xma)
	    rm2 = DBLE(xmb)
	    rp  = DBLE(xp)
	    rm1p  = DBLE(dmap)
	    rm2p  = DBLE(dmbp)
	    rm1m2 = DBLE(dmamb)
	    do 20 j=1,3
		do 10 i=1,3
		    rpiDpj(i,j) = DBLE(piDpj(i,j))
   10		continue
   20	    continue
	    sprec = precx
	    precx = precc
	    call ffxb2a(cb2p,cb2d,cb1,cb0,ca0i,rp,rm1,rm2,rm1p,rm2p,
     +	    	rm1m2,rpiDpj,ier)
	    precx = sprec
	    return
	endif
*  #] real case:
*  #[ p^2 != 0:
	if ( xp .ne. 0 ) then
* 	#[ normal case:
	    call ffclmb(xlam,-xp,-xmb,-xma,dmbp,dmap,dmamb)
	    delsp = -xlam/4
*
*	    the first one is simple...
*
	    cs(1) = 2*cb1*piDpj(1,3)
	    cs(2) = ca0i(2)
	    cb2p = cs(1) + cs(2)
*
*	    the next one ain't.
*
	    cs(1) = ca0i(2)
	    cs(2) = 2*xma*cb0
	    cs(3) = -2*piDpj(1,3)*cb1
	    cs(4) = xma+xmb
	    cs(5) = -xp/3
	    cb2d = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
	    xmax = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),absc(cs(5)))
	    if ( absc(cb2d) .ge. xloss*xmax ) goto 110
  110	    continue
* 	#] give up:
	    cb2p = cb2p*DBLE(1/(2*xp))
	    cb2d = cb2d*(1/DBLE(6))
*  #] p^2 != 0:
*  #[ p^2=0:
	elseif ( dmamb .ne. 0 ) then
	    if ( init.eq.0 ) then
		init = 1
		print *,' '
		print *,'ffcb2a: note: in this case p^2=0 B21 is ',
     +			'returned rather than B2p which is undefined'
		print *,' '
	    endif
	    if ( dmamb .ne. 0 ) then
*		#[ B21:
		llogmm = .FALSE.
*
*		B21 (see thesis, b21.frm)
*
		cs(1) = xma**2/3/dmamb**3*ca0i(1)
		cs(2) = (-xma**2 + xma*xmb - xmb**2/3)/dmamb**3*ca0i(2)
		cs(3) = (5*xma**3/18 - xma*xmb**2/2 + 2*xmb**3/9)
     +			/dmamb**3
		cb21 = cs(1)+cs(2)+cs(3)
		xmax = max(absc(cs(2)),absc(cs(3)))
		if ( absc(cb21).gt.xloss**2*xmax ) goto 160
*
*		ma ~ mb
*
		if ( absc(dmamb).lt.xloss*absc(xma) ) then
		    xlogmm = zfflo1(dmamb/xma,ier)
		else
		    xlogmm = log(xmb/xma)
		endif
		llogmm = .TRUE.
		cs(1) = (xma/dmamb)/6
		cs(2) = (xma/dmamb)**2/3
		cs(3) = (xma/dmamb)**3*xlogmm/3
		cs(4) = -2/DBLE(9) + ca0i(1)/(3*xma)
		cs(5) = -xlogmm/3
		csom = cs(1)+cs(2)+cs(3)+cs(4)+cs(5)
		xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
     +			absc(cs(5)))
		if ( xmxp.lt.xmax ) then
		    xmax = xmxp
		    cb21 = csom
		    if ( absc(cb21).gt.xloss**2*xmax ) goto 160
		endif
*
*		and last try
*
		xlo3 = zfflo3(dmamb/xma,ier)
		cs(1) = (dmamb/xma)**2/6
		cs(2) = (dmamb/xma)/3
		cs(3) = xlo3/(3*(dmamb/xma)**3)
*same		cs(4) = -2/DBLE(9) + ca0i(1)/(3*xma)
		cs(5) = -xlo3/3
		csom = cs(1)+cs(2)+cs(3)+cs(4)+cs(5)
		xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
     +			absc(cs(5)))
		if ( xmxp.lt.xmax ) then
		    xmax = xmxp
		    cb21 = csom
		    if ( absc(cb21).gt.xloss**2*xmax ) goto 160
		endif
*
*		give up
*
  160		continue
*		#] B21:
*		#[ B22:
*
*		B22
*
		cs(1) = +xma/(4*dmamb)*ca0i(1)
		cs(2) = -xmb/(4*dmamb)*ca0i(2)
		cs(3) = (xma+xmb)/8
		cb22 = cs(1) + cs(2) + cs(3)
		xmax = max(absc(cs(2)),absc(cs(3)))
		if ( absc(cb22).gt.xloss*xmax ) goto 210
*
*		second try, close together
*
		if ( .not.llogmm ) then
		    if ( abs(dmamb).lt.xloss*absc(xma) ) then
		    	xlogmm = zfflo1(dmamb/xma,ier)
		    else
		    	xlogmm = log(xmb/xma)
		    endif
		endif
		cs(1) = dmamb*( -1/DBLE(8) - ca0i(1)/(4*xma) )
		cs(2) = dmamb*xlogmm/4
		cs(3) = xma*(xma/dmamb)/4*xlogmm
		cs(4) = xma*( 1/DBLE(4) + ca0i(1)/(2*xma) )
		cs(5) = -xma*xlogmm/2
		csom = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
		xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
     +			absc(cs(5))) 
		if ( xmxp.lt.xmax ) then
		    xmax = xmxp
		    cb22 = csom
		endif
		if ( absc(cb22).gt.xloss*xmax ) goto 210
*
*		give up
*
  210		continue
*		#] B22:
	    else
*
*		ma=mb: simple
*
		cb21 = cb0/3
		cb22 = xma/2*(cb0 + 1)
	    endif
	    cb2d = cb22
	    cb2p = cb21
	endif
*  #] p^2=0:
*###] ffcb2a:
	end
