#include "externals.h"
#include "types.h"


*###[ ffxb1:
	subroutine ffxb1(cb1,cb0,ca0i,xp,xm1,xm2,piDpj,ier)
***#[*comment:***********************************************************
*									*
*	Calculate   1	/	   d^n Q Q(mu)				*
*		 ------	| ------------------------ = B1*p(mu)		*
*		 i pi^2	/ (Q^2-m1^2)((Q+p)^2-m2^2)			*
*									*
*	Input:	cb0	   complex	scalar twopoint function	*
*		ca0i(2)	   complex	scalar onepoint function with	*
*						m1,m2			*
*		xp	   real		p.p in B&D metric		*
*		xm1,2	   real		m_1^2,m_2^2			*
*		piDpj(3,3) real		dotproducts between s1,s2,p	*
*		ier	   integer	digits lost so far		*
*	Output:	cb1	   complex	B1				*
*		ier	   integer	digits lost			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	RealType xp,xm1,xm2,piDpj(3,3)
	ComplexType cb1,cb0,ca0i(2)
*
*	local variables
*
	RealType dm1m2
*
*	common blocks
*
#include "ff.h"
*
*  #] declarations:
*  #[ get differences:
	dm1m2 = xm1 - xm2
*  #] get differences:
*  #[ call ffxb1a:
	call ffxb1a(cb1,cb0,ca0i,xp,xm1,xm2,dm1m2,piDpj,ier)
*  #] call ffxb1a:
*###] ffxb1:
	end


*###[ ffxb1a:
	subroutine ffxb1a(cb1,cb0,ca0i,xp,xm1,xm2,dm1m2,piDpj,ier)
***#[*comment:***********************************************************
*									*
*	Calculate   1	/	   d^n Q Q(mu)				*
*		 ------	| ------------------------ = B1*p(mu)		*
*		 i pi^2	/ (Q^2-m1^2)((Q+p)^2-m2^2)			*
*									*
*	Input:	cb0	   complex	scalar twopoint function	*
*		ca0i(2)	   complex	scalar onepoint function with	*
*						m1,m2			*
*		xp	   real		p.p in B&D metric		*
*		xm1,2	   real		m_1^2,m_2^2			*
*		piDpj(3,3) real		dotproducts between s1,s2,p	*
*		ier	   integer	digits lost so far		*
*	Output:	cb1	   complex	B1				*
*		ier	   integer	digits lost			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	RealType xp,xm1,xm2,dm1m2,piDpj(3,3)
	ComplexType cb1,cb0,ca0i(2)
*
*	local variables
*
	logical lneg
	RealType xmax,absc,s,s1,h,slam,bnd101,bnd105,bnd110,
     +		xma,xmb,x,ax,xlogm,small,dmbma,xprec,xlam,ts2Dp,
     +		xmxp,xlo3,dfflo3
	ComplexType cs(5),cc,csom
	RealType ffbnd,dfflo1
	external ffbnd,dfflo1,dfflo3
	save xprec,bnd101,bnd105,bnd110
*
*	common blocks
*
#include "ff.h"
*
*	statement function
*
	absc(cc) = abs(Re(cc)) + abs(Im(cc))
*
*	data
*
	data xprec /0D0/
*
*  #] declarations:
*  #[ p^2 != 0:
	if ( xp .ne. 0 ) then
* 	#[ normal case:
	    if ( dm1m2 .ne. 0 ) then
		cs(1) = -ca0i(2)
		cs(2) = +ca0i(1)
	    else
		cs(1) = 0
		cs(2) = 0
	    endif
	    cs(3) = +Re(2*piDpj(1,3))*cb0
	    cb1 = cs(1) + cs(2) + cs(3)
	    xmax = max(absc(cs(2)),absc(cs(3)))
	    if ( absc(cb1) .ge. xloss*xmax ) goto 110
* 	#] normal case:
* 	#[ almost equal masses:
	    if ( abs(dm1m2) .le. xloss*xm1 ) then
		cs(2) = Re(dm1m2/xm1)*cs(2)
		cs(1) = -xm2*dfflo1(-dm1m2/xm2,ier)
		cb1 = cs(1) + cs(2) + cs(3)
		xmax = max(absc(cs(2)),absc(cs(3)))
		if ( absc(cb1) .ge. xloss*xmax ) goto 110
*		for the perfectionist (not me (today)):
*		if d0=0 and mu~m1(~m2), then the terms of order
*		(m1^2-m2^2) also cancel.  To patch this I need d0 and mu
	    endif
* 	#] almost equal masses:
* 	#[ p2 -> 0:
	    if ( xloss**2*max(xm1,xm2) .gt. abs(xp) ) then
		if ( xm2.gt.xm1 ) then
		    xma = xm1
		    xmb = xm2
		    ts2Dp = +2*piDpj(2,3)
		    lneg = .FALSE.
		else
		    xma = xm2
		    xmb = xm1
		    ts2Dp = -2*piDpj(1,3)
		    lneg = .TRUE.
		endif
	    else
		goto 100
	    endif
*
*	    We found a situation in which p2 is much smaller than
*	    the masses.
*
	    dmbma = abs(dm1m2)
	    if ( xma.eq.0 ) then
		xlogm = 1
	    elseif ( dmbma .gt. xloss*xmb ) then
		xlogm = log(xmb/xma)
	    else
		xlogm = dfflo1(-dmbma/xma,ier)
	    endif
	    xlam =  (dmbma-xp)**2 - 4*xma*xp
	    if ( xlam.gt.0 ) then
*		#[ real roots:
		slam = sqrt(xlam)
		small = xp*(-2*(xma+xmb) + xp)/(slam+dmbma)
		h = slam+2*piDpj(1,2)
		cs(1) = xlogm*xma*(4*xmb*(small-xp) + (small-xp)**2)/(2*
     +			(slam+dmbma)*h)
		if ( xprec.ne.precx ) then
		    xprec = precx
		    bnd101 = ffbnd(2,1,xinfac)
		    bnd105 = ffbnd(2,5,xinfac)
		    bnd110 = ffbnd(2,10,xinfac)
		endif
		x = xp/slam
		ax = abs(x)
		if ( ax.gt.bnd110 ) then
		    s = x*(xinfac(12) + x*(xinfac(13) + x*(xinfac(14) +
     +			x*(xinfac(15) + x*xinfac(16) ))))
		else
		    s = 0
		endif
		if ( ax.gt.bnd105 ) then
		    s = x*(xinfac(7) + x*(xinfac(8) + x*(xinfac(9) +
     +			x*(xinfac(10) + x*(xinfac(11) + s )))))
		endif
		if ( ax.gt.bnd101) then
		    s = x*(xinfac(3) + x*(xinfac(4) + x*(xinfac(5) +
     +			x*(xinfac(6) + s ))))
		endif
		s = x**2*(.5D0 + s)
		h = ts2Dp + slam
		s1 = 2*xp/h*(s + x)
		h = -4*xp**2*xmb/(slam*h**2) - s + s1
		if ( abs(h) .lt. .1 ) then
		    cs(2) = dmbma*slam/xp*dfflo1(h,ier)
		else
		    goto 100
		endif
		if ( lneg ) then
		    cs(1) = -cs(1)
		    cs(2) = -cs(2)
		endif
		cs(3) = -Re(xp)*cb0
		cb1 = cs(1) + cs(2) + cs(3)
		xmax = max(absc(cs(2)),absc(cs(3)))
		if ( absc(cb1) .gt. xloss*xmax) goto 110
*
*		this still occurs in the case xp << dmamb << xma,
*		with a cancellation of order dmamb/xma between cs1 and
*		cs2; as the standard model does not contain these kind
*		of doublets I leave this as an exercise for the
*		reader...
*
*		#] real roots:
	    else
*		#[ imaginary roots:
*		#] imaginary roots:
	    endif
* 	#] p2 -> 0:
* 	#[ give up:
*
*	    give up...
*
  100	    continue
  110	    continue
* 	#] give up:
	    cb1 = cb1*(1/Re(2*xp))
*  #] p^2 != 0:
*  #[ p^2=0, m1 != m2:
	elseif ( dm1m2 .ne. 0 ) then
	    cs(1) = +Re(xm2/(2*dm1m2**2))*(ca0i(2)+Re(xm2)/2)
	    cs(2) = -Re(xm1/(2*dm1m2**2))*(ca0i(1)+Re(xm1)/2)
	    cs(3) = +ca0i(2)*(1/Re(dm1m2))
	    cb1 = cs(1) + cs(2) + cs(3)
	    xmax = max(absc(cs(1)),absc(cs(2)),absc(cs(3)))
	    if ( absc(cb1).ge.xloss**2*xmax ) goto 120
*
*	    m1 ~ m2, see b21.frm
*
	    if ( abs(dm1m2).lt.xloss*xm1 ) then
		xlogm = dfflo1(dm1m2/xm1,ier)
	    else
		xlogm = log(xm2/xm1)
	    endif
	    cs(1) = -(xm1/dm1m2)/2
	    cs(2) = -xlogm/2*(xm1/dm1m2)**2
	    cs(3) = +1/Re(4) - ca0i(1)*Re(1/(2*xm1))
	    cs(4) = xlogm/2
	    csom = cs(1) + cs(2) + cs(3) + cs(4)
	    xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)))
	    if ( xmxp.lt.xmax ) then
		xmax = xmxp
		cb1 = csom
		if ( absc(cb1).gt.xloss**2*xmax ) goto 120
	    endif
*
*	    better
*
	    xlo3 = dfflo3(dm1m2/xm1,ier)
	    cs(1) = -(dm1m2/xm1)**2/4
	    cs(2) = -(dm1m2/xm1)/2
	    cs(3) = -xlo3/(dm1m2/xm1)**2/2
	    cs(4) = xlo3/2
	    cs(5) = 1/Re(2) - ca0i(1)*Re(1/(2*xm1))
	    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
		cb1 = csom
		if ( absc(cb1).gt.xloss**2*xmax ) goto 120
	    endif
*
*	    give up
*
  120	    continue
*  #] p^2=0, m1 != m2:
*  #[ p^2=0, m1 == m2:
	else
	    cb1 = -cb0/2
	endif
*  #] p^2=0, m1 == m2:
*###] ffxb1a:
	end
