
* file aaxex  23-sep-1990

*###[ aaxex:
	subroutine aaxex(ccxi,cdxi,cexi,d0,xmm,xpi,level,ier)
***#[*comment:***********************************************************
*									*
*	Calculation of four point tensor integrals.  Just a wrapper	*
*	for ffxdx nowadays, see there for the real description.		*
*									*
*	Input:	xpi(20)	  real	   the same as in FF			*
*		level	  integer  rank of tensor integral		*
*	Output:	ccxi(30)  complex  cc0(1),cc1( ),[cc2( ),cc3( ) ]  i,j	*
*		cdxi(55)  complex  cd0(1),cd1(3),cd2(7),[cd3(13)]	*
*					i=1,2,3,4,5			*
*		cexi(35)  complex  ce0(1),ce1(4),ce2(10),ce3(20)	*
*		ier	  integer  FF error flag			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION xpi(20),d0,xmm
	DOUBLE COMPLEX ccxi(30),cdxi(55),cexi(35)
*
*	local variables
*
	DOUBLE COMPLEX caxi(5),cbxi(30)
	DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
*
*  #] declarations:
*  #[ call ffxex:
*
	call ffxex(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,mexi,
     +		d0,xmm,xpi,level,ier)
*
*  #] call ffxex:
*###] aaxex:
	end
*###[ ffxex:
	subroutine ffxex(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,
     +		mexi,d0,xmm,xpi,level,ier)
***#[*comment:***********************************************************
*									*
*	Calculation of five-point formfactors, notation is defined in	*
*	adapt.prc and covdec5.prc (WW), leaving out the d_ terms	*
*	adapted by GJvO 4-apr-1995					*
*									*
*	cexi(1)		E0						*
*	cexi(2-5)       E1i, coeff of pi				*
*	cexi(6-15)	E2i, coeff of p1p1,p2p2,p3p3,p4p4,(p1p2),	*
*			(p1p3),(p1p4),(p2p3),(p2p4),(p3p4)		*
*	cexi(16-35)	E2i, coeff of p1p1p1,p2p2p2,p3p3p3,p4p4p4,	*
*			(p1p1p2),(p1p1p3),(p1p1p4),(p2p2p1),(p2p2p3),	*
*			(p2p2p4),(p3p3p1),(p3p3p2),(p3p3p4),(p4p4p1),	*
*			(p4p4p2),(p4p4p3),(p1p2p3),(p1p2p4),(p1p3p4),	*
*			(p2p3p4)					*
*	cdxi(55)	5*(D0(1),D1i(3),D2i(7)), s_i missing		*
*	ccxi(30)	15*(C0(1),C1i(2))  ??????			*
*	cbxi(30)	not used yet					*
*	caxi(5)		not used yet					*
*	m[abcde]xi()	if c[abcde]xi() were written as a sum of stable	*
*			terms, the largest term in that sum, i.e., the	*
*			accuracy of c[abcde]xi() is precc*m[abcde]xi	*
*									*
*	Input:	xpi(20)	  real	   the same as in FF			*
*		level	  integer  rank of tensor integral		*
*	Output:	ccxi(30)  complex  cc0(1),cc1( ),[cc2( ),cc3( ) ]  i,j	*
*		cdxi(55)  complex  cd0(1),cd1(3),cd2(7),[cd3(13)]	*
*					i=1,2,3,4,5			*
*		cexi(35)  complex  ce0(1),ce1(5),ce2(10),ce3(20)	*
*		ier	  integer  FF error flag			*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION xpi(20),d0,xmm
	DOUBLE COMPLEX caxi(5),cbxi(30),ccxi(30),cdxi(55),cexi(35)
	DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
*
*	local variables
*
	integer i,j,dl,iplace(11,5),ier0,ier1
	DOUBLE PRECISION xpj(13),absc,big
        DOUBLE COMPLEX cd0i(5),cdxj(120),ccxj(140),cbxj(60),caxj(20),cc
	DOUBLE PRECISION mdxj(120),mcxj(140),mbxj(60),maxj(20)
	save iplace
*
*	common blocks
*
	include 'ff.h'
	include 'aa.h'
*
*	statement functions
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*
*	data
*
	data iplace /
     +          2,3,4,5, 07,08,09,15, 12,13, 17,
     +          1,3,4,5, 11,08,09,10, 14,13, 18,
     +          1,2,4,5, 06,12,09,10, 14,15, 19,
     +          1,2,3,5, 06,07,13,10, 11,15, 20,
     +          1,2,3,4, 06,07,08,14, 11,12, 16/
*
*  #] declarations:
*  #[ init:
*
*	initialize to something ridiculous so that one immediately 
*	notices when it is accidentally used...
*
	big = 1/(1d20*xclogm)
	do i=1,5
	    caxi(i) = big
	enddo
	do i=1,30
	    cbxi(i) = big
	enddo
	do i=1,30
	    ccxi(i) = big
	enddo
	do i=1,55
	    cdxi(i) = big
	enddo
	do i=1,35
	    cexi(i) = big
	enddo
*
*  #] init:
*  #[ level 0: E0, and kinematical quantities for 5 point PV-red
*	E0-function (ff)
*
	ldot = .TRUE.
	ier1 = ier
	call ffxe0(cexi(1),cd0i,xpi,ier1)
	mexi(1) = absc(cexi(1))*DBLE(10)**mod(ier1,50)
	do i=1,5
	    cdxi(1+11*(i-1)) = cd0i(i)
	enddo
*
	if ( awrite ) then
	    print *,'    '
	    print *,'aaxex : level 0, imported from ff '
	    print *,'E0    = ',cexi(1),mexi(1)
	    print *,'D0(1) = ',cd0i(1)
	    print *,'D0(2) = ',cd0i(2)
	    print *,'D0(3) = ',cd0i(3)
	    print *,'D0(4) = ',cd0i(4)
	    print *,'D0(5) = ',cd0i(5)
	    print *,'xpi used:'
	    do i =1,15
		print *,i,xpi(i)
	    enddo
	    print *,'imported stuff via ff.h:'
	    print *,'    kin determinat = ',fdel4
	    print *,'dotpr(1,1)= ',fpij5(6,6)
	    print *,'dotpr(2,2)= ',fpij5(7,7)
	    print *,'dotpr(1,2)= ',fpij5(6,7)
	endif
*
	if (level .eq. 0) return
*
*  #] level 0: E0, and kinematical quantities for 5 point PV-red
*  #[ need D-functions till d-level=(level-1):
	dl=level-1
*
*	go trough the 5 different cancellation patterns
*
	if ( awrite ) then
	   print *,'     '
	   print *,'------>underlying D-functions up to level:',dl
	endif
	xpj(12) = 0
	xpj(13) = 0
	do j=1,5
*	    D(j) is the D0 with leg j missing.
	    do i=1,11
		xpj(i) = xpi(iplace(i,j))
	    enddo
*	    note that we recompute the D0 (or get it from cache)
	    ier0 = ier
	    call ffxdx(caxj(1+4*(j-1)),maxj(1+4*(j-1)),cbxj(1+12*(j-1))
     +		,mbxj(1+12*(j-1)),ccxj(1+28*(j-1)),mcxj(1+28*(j-1)),
     +		cdxj(1+24*(j-1)),mdxj(1+24*(j-1)),d0,xmm,xpj,dl,ier0)
	    ier1 = max(ier1,ier0)
	enddo
	ier = ier1
	if ( awrite ) then
	    print *,' '
	    print *,'ier = ',ier
	    print *,'---->end of D-function output--------------------'
	endif
	if ( atest ) then
*	    (although these should come from cache) (but don't yet!)
	    do i=1,5
		if ( xloss*10d0**(-mod(ier,50))*absc(cd0i(i)-
     +			cdxj(1+24*(i-1))) .gt. precx*absc(cd0i(i)) )
     +			then
		    print *,'aaxex: error: D0 does not agree with ',
     +			'recomputed: ',cd0i(i),cdxj(1+24*(i-1)),ier
		endif
	    enddo
	endif
*
*  #] need D-functions till d-level=(level-1):
*  #[ break to let ffzez tie in:
*
*	convert ??xj to ??xi
*
	call ffeji(cdxi,mdxi,ccxi,mcxi,cbxi,mbxi,caxi,maxi,
     +		cdxj,mdxj,ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
*
*	and call the real routine for the rest
*
	call ffxexp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,mexi,
     +		xpi,level,ier)
*
*  #] break to let ffzez tie in:
	end
	subroutine ffxexp(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,cexi,
     +		mexi,xpi,level,ier)
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier,level
	DOUBLE PRECISION xpi(20),d0,xmm
	DOUBLE COMPLEX caxi(5),cbxi(30),ccxi(30),cdxi(55),cexi(35)
	DOUBLE PRECISION maxi(5),mbxi(30),mcxi(30),mdxi(55),mexi(35)
*
*	local variables
*
	integer i,j,ier0,ier1,ij2k(4,4),m2ijk(3,20)
	DOUBLE PRECISION xi5(10),f1,f2,f3,f4,absc
        DOUBLE COMPLEX R(70),cd0i(5),cd1ij(3,5),ce2ij(4,4),ce3ijk(4,4,4)
     +		,cd2ijk(3,3,5),cd2i(5),cxy(70),cc,rg(4),cexj(39)
	DOUBLE PRECISION mdxj(120),mcxj(140),mbxj(60),maxj(20),
     +		del3ij(5,5),del4i(4)
	save ij2k,m2ijk
*
*	common blocks
*
	include 'ff.h'
	include 'aa.h'
*
*	statement functions
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*
*	data
*
	data ij2k / 6,10,11,12,
     +		   10, 7,13,14,
     +		   11,13, 8,15,
     +		   12,14,15, 9/
	data m2ijk /1,1,1,  2,2,2,
     +		    3,3,3,  4,4,4,
     +		    1,1,2,  1,1,3,
     +		    1,1,4,  2,2,1,
     +		    2,2,3,  2,2,4,
     +		    3,3,1,  3,3,2,
     +		    3,3,4,  4,4,1,
     +		    4,4,2,  4,4,3,
     +		    1,2,3,  1,2,4,
     +		    1,3,4,  2,3,4/
*
*  #] declarations:
*  #[ kinematical quatities for 5pv-red:
*
*	inverse kinematical matrix xi5  (4X4)
*
	call aaxi5(xi5,ier)
*
*	AAs f-functions:
*
	f1 = 2*fpij5(1,6)
	f2 = 2*fpij5(1,7)
	f3 = 2*fpij5(1,8)
	f4 = 2*fpij5(1,9)
*
*  #] kinematical quatities for 5pv-red:
*  #[ level 1: E11,E12,E13,E14,D0(I)
	do i=1,5
	    cd0i(i) = cdxi(1+11*(i-1))
	enddo
	call ffxe1(cexi(2),cexi(1),del3ij,del4i,cd0i,xpi,fpij5,fdel4,
     +		ier)
	if ( atest ) then
*
*	PV-reduction
*
	    R(1) = (f1*cexi(1) + cd0i(2) - cd0i(1))/2
	    R(2) = (f2*cexi(1) + cd0i(3) - cd0i(2))/2
	    R(3) = (f3*cexi(1) + cd0i(4) - cd0i(3))/2
	    R(4) = (f4*cexi(1) + cd0i(5) - cd0i(4))/2
	    cxy(2)=xi5(1)*R(1)+xi5(5)*R(2)+xi5(6) *R(3)+xi5(7) *R(4)
	    cxy(3)=xi5(5)*R(1)+xi5(2)*R(2)+xi5(8) *R(3)+xi5(9) *R(4)
	    cxy(4)=xi5(6)*R(1)+xi5(8)*R(2)+xi5(3) *R(3)+xi5(10)*R(4)
	    cxy(5)=xi5(7)*R(1)+xi5(9)*R(2)+xi5(10)*R(3)+xi5(4) *R(4)
	    do i=2,5
		if ( xloss*10d0**(-mod(ier,50))*absc(cexi(i)-cxy(i))
     +			.gt. precc*absc(cxy(i)) ) then
		    print *,'aaxex: error: E1 from ffxe1 is not correct'
     +			,i,cexi(i),cxy(i),ier
		endif
	    enddo
	endif
	if ( awrite ) then
	    print *,'     '
	    print *,'aaxex: level 1: id,nevent ',id,nevent
	    print *,'E11=',cexi(2),ier
	    print *,'E12=',cexi(3),ier
	    print *,'E13=',cexi(4),ier
	    print *,'E14=',cexi(5),ier
	endif

	if ( level.le.1 ) return
*
*  #] level 1:
*  #[ level 2: E21,E22,E23,E24,E25,E26,E27,E28,E29,E210
*
*	D11(I),D12(I),D13(I) need 5 diff D1-functions
*
	do i=1,5
	    j = 2 + (i-1)*11
	    cd1ij(1,i)=cdxi(j)
	    cd1ij(2,i)=cdxi(j+1)
	    cd1ij(3,i)=cdxi(j+2)
	enddo
*
*	GJ reduction:
*
	call ffxe2(ce2ij,cexi(2),cexi(1),cd1ij,cd0i,xpi,fpij5,
     +          del3ij,del4i,fdel4,ier)
	if ( atest ) then
*
*	PV-reduction
*
	R(11) = (f1*cexi(2) + cd1ij(1,2) + cd0i(1) )/2
	R(12) = (f2*cexi(2) + cd1ij(1,3) - cd1ij(1,2))/2
	R(13) = (f3*cexi(2) + cd1ij(1,4) - cd1ij(1,3))/2
	R(14) = (f4*cexi(2) + cd1ij(1,5) - cd1ij(1,4))/2
*
	R(15) = (f1*cexi(3) + cd1ij(1,2) - cd1ij(1,1))/2
	R(16) = (f2*cexi(3) + cd1ij(2,3) - cd1ij(1,2))/2
	R(17) = (f3*cexi(3) + cd1ij(2,4) - cd1ij(2,3))/2
	R(18) = (f4*cexi(3) + cd1ij(2,5) - cd1ij(2,4))/2
*
	R(19) = (f1*cexi(4) + cd1ij(2,2) - cd1ij(2,1))/2
	R(20) = (f2*cexi(4) + cd1ij(2,3) - cd1ij(2,2))/2
	R(21) = (f3*cexi(4) + cd1ij(3,4) - cd1ij(2,3))/2
	R(22) = (f4*cexi(4) + cd1ij(3,5) - cd1ij(3,4))/2
*
	R(23) = (f1*cexi(5) + cd1ij(3,2) - cd1ij(3,1))/2
	R(24) = (f2*cexi(5) + cd1ij(3,3) - cd1ij(3,2))/2
	R(25) = (f3*cexi(5) + cd1ij(3,4) - cd1ij(3,3))/2
	R(26) = (f4*cexi(5)            - cd1ij(3,4))/2
*
	cexi(6) = xi5(1)*R(11)+xi5(5)*R(12)+xi5(6) *R(13)+xi5(7) *R(14)
	cexi(7) = xi5(5)*R(15)+xi5(2)*R(16)+xi5(8) *R(17)+xi5(9) *R(18)
	cexi(8) = xi5(6)*R(19)+xi5(8)*R(20)+xi5(3) *R(21)+xi5(10)*R(22)
	cexi(9) = xi5(7)*R(23)+xi5(9)*R(24)+xi5(10)*R(25)+xi5(4) *R(26)
	cexi(10)= xi5(5)*R(11)+xi5(2)*R(12)+xi5(8) *R(13)+xi5(9) *R(14)
	cexi(11)= xi5(6)*R(11)+xi5(8)*R(12)+xi5(3) *R(13)+xi5(10)*R(14)
	cexi(12)= xi5(7)*R(11)+xi5(9)*R(12)+xi5(10)*R(13)+xi5(4) *R(14)
	cexi(13)= xi5(6)*R(15)+xi5(8)*R(16)+xi5(3) *R(17)+xi5(10)*R(18)
	cexi(14)= xi5(7)*R(15)+xi5(9)*R(16)+xi5(10)*R(17)+xi5(4) *R(18)
	cexi(15)= xi5(7)*R(19)+xi5(9)*R(20)+xi5(10)*R(21)+xi5(4) *R(22)
*
	if ( atest ) then
*
*	    redundancy check
*
	cxy(10) = xi5(1)*R(15)+xi5(5)*R(16)+xi5(6)*R(17)+xi5(7) *R(18)
	cxy(11) = xi5(1)*R(19)+xi5(5)*R(20)+xi5(6)*R(21)+xi5(7) *R(22)
	cxy(12) = xi5(1)*R(23)+xi5(5)*R(24)+xi5(6)*R(25)+xi5(7) *R(26)
	cxy(13) = xi5(5)*R(19)+xi5(2)*R(20)+xi5(8)*R(21)+xi5(9) *R(22)
	cxy(14) = xi5(5)*R(23)+xi5(2)*R(24)+xi5(8)*R(25)+xi5(9) *R(26)
	cxy(15) = xi5(6)*R(23)+xi5(8)*R(24)+xi5(3)*R(25)+xi5(10)*R(26)
	    do i=10,15
		if ( absc(cxy(i)-cexi(i))  .gt. .1d-3  ) then
		    print *,'redundancy check at level 2 failed'
		endif
	    enddo
*###] :
	endif
*
*	check against GJ
*
	do i=1,4
	    do j=1,4
		if ( xloss*10d0**(-mod(ier,50))*absc(ce2ij(i,j) -
     +			cexi(ij2k(i,j))) .gt. precc*absc(ce2ij(i,j)) )
     +			then
		    print *,'ffxdxp: error: GJ does not agree with PV:',
     +			i,j,ce2ij(i,j),cexi(ij2k(i,j)),ce2ij(i,j) -
     +			cexi(ij2k(i,j)),ier
		endif
	    enddo
	enddo
	endif
*
*	copy to AAs arrays
*
	do i=1,4
	    do j=1,4
		cexi(ij2k(j,i)) = ce2ij(j,i)
	    enddo
	enddo
*
	if ( awrite ) then
	    print *,'     '
	    print *,'aaxex: level 2: id,nevent ',id,nevent
	    print *,'E21 =',cexi(6),ier
	    print *,'E22 =',cexi(7),ier
	    print *,'E23 =',cexi(8),ier
	    print *,'E24 =',cexi(9),ier
	    print *,'E25 =',cexi(10),ier
	    print *,'E26 =',cexi(11),ier
	    print *,'E27 =',cexi(12),ier
	    print *,'E28 =',cexi(13),ier
	    print *,'E28 =',cexi(14),ier
	    print *,'E210=',cexi(15),ier
	endif
*
	if (level .le. 2) return
*
*  #] level 2:
*###[ : level 3 :
*
*	first recast the D2's to a more useful form
*
	do 15 i=1,5
	    j = 5 +(i-1)*11
	    cd2ijk(1,1,i) = cdxi(j)
	    cd2ijk(2,2,i) = cdxi(j+1)
	    cd2ijk(3,3,i) = cdxi(j+2)
	    cd2ijk(1,2,i) = cdxi(j+3)
	    cd2ijk(2,1,i) = cdxi(j+3)
	    cd2ijk(1,3,i) = cdxi(j+4)
	    cd2ijk(3,1,i) = cdxi(j+4)
	    cd2ijk(2,3,i) = cdxi(j+5)
	    cd2ijk(3,2,i) = cdxi(j+5)
	    cd2i(i) = cdxi(j+6)
 15	continue
*
*	FF function
*
	call ffxe3(ce3ijk,ce2ij,cexi(2),cexi(1), cd2ijk,cd2i,cd1ij,cd0i,
     +		xpi,fpij5, del3ij,del4i,fdel4, ier)
*
*	copy FF structs to AA
*
	do i=16,35
	    cexi(i) = ce3ijk(m2ijk(1,i),m2ijk(2,i),m2ijk(3,i))
	enddo
	if ( atest ) then
*
*
*		PV-reduction
*	g-terms
	rg(1)=1/2d0*( f1*cexi(6)+f2*cexi(10)+f3*cexi(11)+f4*cexi(12) )
	rg(2)=1/2d0*( f1*cexi(10)+f2*cexi(7)+f3*cexi(13)+f4*cexi(14) )
	rg(3)=1/2d0*( f1*cexi(11)+f2*cexi(13)+f3*cexi(8)+f4*cexi(15) )
	rg(4)=1/2d0*( f1*cexi(12)+f2*cexi(14)+f3*cexi(15)+f4*cexi(9) )
*
	cexj(36)=xpi(1)*cexj(2)-1/2d0*cd0i(1) -rg(1)
	cexj(37)=xpi(1)*cexj(3)+1/2d0*cd1ij(1,1)-rg(2)
	cexj(38)=xpi(1)*cexj(4)+1/2d0*cd1ij(2,1)-rg(3)
	cexj(39)=xpi(1)*cexj(5)+1/2d0*cd1ij(3,1)-rg(4)
*
*	terms ~pipi
*	1)
	R(31)=1/2d0*(f1*cexj(6)+cd2ijk(1,1,2)-cd0i(1)) -2*cexj(36)
	R(32)=1/2d0*(f2*cexj(6)+cd2ijk(1,1,3)-cd2ijk(1,1,2))
	R(33)=1/2d0*(f3*cexj(6)+cd2ijk(1,1,4)-cd2ijk(1,1,3))
	R(34)=1/2d0*(f4*cexj(6)+cd2ijk(1,1,5)-cd2ijk(1,1,4))
*	2)
	R(35)=1/2d0*(f1*cexj(7)+cd2ijk(1,1,2)-cd2ijk(1,1,1))
	R(36)=1/2d0*(f2*cexj(7)+cd2ijk(2,2,3)-cd2ijk(1,1,2)) -2*cexj(37)
	R(37)=1/2d0*(f3*cexj(7)+cd2ijk(2,2,4)-cd2ijk(2,2,3))
	R(38)=1/2d0*(f4*cexj(7)+cd2ijk(2,2,5)-cd2ijk(2,2,4))
*	3)
	R(39)=1/2d0*(f1*cexj(8)+cd2ijk(2,2,2)-cd2ijk(2,2,1))
	R(40)=1/2d0*(f2*cexj(8)+cd2ijk(2,2,3)-cd2ijk(2,2,2))
	R(41)=1/2d0*(f3*cexj(8)+cd2ijk(3,3,4)-cd2ijk(2,2,3)) -2*cexj(38)
	R(42)=1/2d0*(f4*cexj(8)+cd2ijk(3,3,5)-cd2ijk(3,3,4))
*	4)
	R(43)=1/2d0*(f1*cexj(9)+cd2ijk(3,3,2)-cd2ijk(3,3,1))
	R(44)=1/2d0*(f2*cexj(9)+cd2ijk(3,3,3)-cd2ijk(3,3,2))
	R(45)=1/2d0*(f3*cexj(9)+cd2ijk(3,3,4)-cd2ijk(3,3,3))
	R(46)=1/2d0*(f4*cexj(9)         -cd2ijk(3,3,4) ) -2*cexj(39)
*
*	terms ~p1pi
*	1)
	R(47)=1/2d0*(f1*cexj(10)+cd2ijk(1,1,2)+cd1ij(1,1)) -cexj(37)
	R(48)=1/2d0*(f2*cexj(10)+cd2ijk(1,2,3)-cd2ijk(1,1,2)) -cexj(36)
	R(49)=1/2d0*(f3*cexj(10)+cd2ijk(1,2,4)-cd2ijk(1,2,3))
	R(50)=1/2d0*(f4*cexj(10)+cd2ijk(1,2,5)-cd2ijk(1,2,4))
*	2)
	R(51)=1/2d0*(f1*cexj(11)+cd2ijk(1,2,2)+cd1ij(2,1) ) -cexj(38)
	R(52)=1/2d0*(f2*cexj(11)+cd2ijk(1,2,3)-cd2ijk(1,2,2))
	R(53)=1/2d0*(f3*cexj(11)+cd2ijk(1,3,4)-cd2ijk(1,2,3)) -cexj(36)
	R(54)=1/2d0*(f4*cexj(11)+cd2ijk(1,3,5)-cd2ijk(1,3,4))
*	3)
	R(55)=1/2d0*(f1*cexj(12)+cd2ijk(1,3,2)+cd1ij(3,1) ) -cexj(39)
	R(56)=1/2d0*(f2*cexj(12)+cd2ijk(1,3,3)-cd2ijk(1,3,2))
	R(57)=1/2d0*(f3*cexj(12)+cd2ijk(1,3,4)-cd2ijk(1,3,3))
	R(58)=1/2d0*(f4*cexj(12)         -cd2ijk(1,3,4) ) -cexj(36)
*
*	terms ~p2pi
*	1)
	R(59)=1/2d0*(f1*cexj(13)+cd2ijk(1,2,2)-cd2ijk(1,2,1))
	R(60)=1/2d0*(f2*cexj(13)+cd2ijk(2,2,3)-cd2ijk(1,2,2)) -cexj(38)
	R(61)=1/2d0*(f3*cexj(13)+cd2ijk(2,3,4)-cd2ijk(2,2,3)) -cexj(37)
	R(62)=1/2d0*(f4*cexj(13)+cd2ijk(2,3,5)-cd2ijk(2,3,4))
*	2)
	R(63)=1/2d0*(f1*cexj(14)+cd2ijk(1,3,2)-cd2ijk(1,3,1))
	R(64)=1/2d0*(f2*cexj(14)+cd2ijk(2,3,3)-cd2ijk(1,3,2)) -cexj(39)
	R(65)=1/2d0*(f3*cexj(14)+cd2ijk(2,3,4)-cd2ijk(2,3,3))
	R(66)=1/2d0*(f4*cexj(14)         -cd2ijk(2,3,4) ) -cexj(37)
*
*	terms ~p3pi
*	1)
	R(67)=1/2d0*(f1*cexj(15)+cd2ijk(2,3,2)-cd2ijk(2,3,1))
	R(68)=1/2d0*(f2*cexj(15)+cd2ijk(2,3,3)-cd2ijk(2,3,2))
	R(69)=1/2d0*(f3*cexj(15)+cd2ijk(3,3,4)-cd2ijk(2,3,3)) -cexj(39)
	R(70)=1/2d0*(f4*cexj(15)         -cd2ijk(3,3,4) ) -cexj(38)

	cexj(16)=xi5(1)*R(31)+xi5(5)*R(32)+xi5(6)*R(33)+xi5(7) *R(34)
	cexj(17)=xi5(5)*R(35)+xi5(2)*R(36)+xi5(8)*R(37)+xi5(9) *R(38)
	cexj(18)=xi5(6)*R(39)+xi5(8)*R(40)+xi5(3)*R(41)+xi5(10)*R(42)
	cexj(19)=xi5(7)*R(43)+xi5(9)*R(44)+xi5(10)*R(45)+xi5(4)*R(46)
	cexj(20)=xi5(5)*R(31)+xi5(2)*R(32)+xi5(8)*R(33)+xi5(9) *R(34)
	cexj(21)=xi5(6)*R(31)+xi5(8)*R(32)+xi5(3)*R(33)+xi5(10)*R(34)
	cexj(22)=xi5(7)*R(31)+xi5(9)*R(32)+xi5(10)*R(33)+xi5(4)*R(34)
	cexj(23)=xi5(1)*R(35)+xi5(5)*R(36)+xi5(6)*R(37)+xi5(7) *R(38)
	cexj(24)=xi5(6)*R(35)+xi5(8)*R(36)+xi5(3)*R(37)+xi5(10)*R(38)
	cexj(25)=xi5(7)*R(35)+xi5(9)*R(36)+xi5(10)*R(37)+xi5(4)*R(38)
	cexj(26)=xi5(1)*R(39)+xi5(5)*R(40)+xi5(6)*R(41)+xi5(7) *R(42)
	cexj(27)=xi5(5)*R(39)+xi5(2)*R(40)+xi5(8)*R(41)+xi5(9) *R(42)
	cexj(28)=xi5(7)*R(39)+xi5(9)*R(40)+xi5(10)*R(41)+xi5(4)*R(42)
	cexj(29)=xi5(1)*R(43)+xi5(5)*R(44)+xi5(6)*R(45)+xi5(7) *R(46)
	cexj(30)=xi5(5)*R(43)+xi5(2)*R(44)+xi5(8)*R(45)+xi5(9) *R(46)
	cexj(31)=xi5(6)*R(43)+xi5(8)*R(44)+xi5(3)*R(45)+xi5(10)*R(46)
	cexj(32)=xi5(6)*R(47)+xi5(8)*R(48)+xi5(3)*R(49)+xi5(10)*R(50)
	cexj(33)=xi5(7)*R(47)+xi5(9)*R(48)+xi5(10)*R(49)+xi5(4)*R(50)
	cexj(34)=xi5(7)*R(51)+xi5(9)*R(52)+xi5(10)*R(53)+xi5(4)*R(54)
	cexj(35)=xi5(7)*R(59)+xi5(9)*R(60)+xi5(10)*R(61)+xi5(4)*R(62)
*

*###[ : redundancy check
	cxy(20)=xi5(1)*R(47)+xi5(5)*R(48)+xi5(6)*R(49)+xi5(7) *R(50)
	cxy(21)=xi5(1)*R(51)+xi5(5)*R(52)+xi5(6)*R(53)+xi5(7) *R(54)
	cxy(22)=xi5(1)*R(55)+xi5(5)*R(56)+xi5(6)*R(57)+xi5(7) *R(58)
	cxy(23)=xi5(5)*R(47)+xi5(2)*R(48)+xi5(8)*R(49)+xi5(9) *R(50)
	cxy(24)=xi5(5)*R(59)+xi5(2)*R(60)+xi5(8)*R(61)+xi5(9) *R(62)
	cxy(25)=xi5(5)*R(63)+xi5(2)*R(64)+xi5(8)*R(65)+xi5(9) *R(66)
	cxy(26)=xi5(6)*R(51)+xi5(8)*R(52)+xi5(3)*R(53)+xi5(10)*R(54)
	cxy(27)=xi5(6)*R(59)+xi5(8)*R(60)+xi5(3)*R(61)+xi5(10)*R(62)
	cxy(28)=xi5(6)*R(67)+xi5(8)*R(68)+xi5(3)*R(69)+xi5(10)*R(70)
	cxy(29)=xi5(7)*R(55)+xi5(9)*R(56)+xi5(10)*R(57)+xi5(4)*R(58)
	cxy(30)=xi5(7)*R(63)+xi5(9)*R(64)+xi5(10)*R(65)+xi5(4)*R(66)
	cxy(31)=xi5(7)*R(67)+xi5(9)*R(68)+xi5(10)*R(69)+xi5(4)*R(70)
	cxy(32)=xi5(5)*R(51)+xi5(2)*R(52)+xi5(8)*R(53)+xi5(9) *R(54)
	cxy(32)=xi5(1)*R(59)+xi5(5)*R(60)+xi5(6)*R(61)+xi5(7) *R(62)
	cxy(33)=xi5(5)*R(55)+xi5(2)*R(56)+xi5(8)*R(57)+xi5(9) *R(58)
	cxy(33)=xi5(1)*R(63)+xi5(5)*R(64)+xi5(6)*R(65)+xi5(7) *R(66)
	cxy(34)=xi5(6)*R(55)+xi5(8)*R(56)+xi5(3)*R(57)+xi5(10)*R(58)
	cxy(34)=xi5(1)*R(67)+xi5(5)*R(68)+xi5(6)*R(69)+xi5(7) *R(70)
	cxy(35)=xi5(6)*R(63)+xi5(8)*R(64)+xi5(3)*R(65)+xi5(10)*R(66)
	cxy(35)=xi5(5)*R(67)+xi5(2)*R(68)+xi5(8)*R(69)+xi5(9) *R(70)
	do 16 i=20,35
	   if ( absc(cxy(i)-cexj(i)) .gt. .1d-3  ) then
	      print *,'ffxdxp: redundancy check at level 3 failed ',
     +	      	i,cxy(i),cexj(i),cxy(i)-cexj(i),ier
	      return
	   endif
 16	continue
*	
*	check against FF
*	
	do 17 i=16,35
	    if ( xloss*10d0**(-mod(ier,50))*absc(cexi(i)-cexj(i)) .gt. 
     +	    		precc*absc(cexi(i)) ) then
	    	print *,'ffxexp: error: FF disagrees with PV: ',i,
     +	    		cexi(i),cexj(i),cexi(i)-cexj(i),ier
	    endif
   17	continue
*###] :
	endif
	if ( awrite ) then
*###[ : print level 3
	    print *,'     '
	    print *,'aaxex : level 3 '
	    print *,'E31 =',cexi(16)
	    print *,'E32 =',cexi(17)
	    print *,'E33 =',cexi(18)
	    print *,'E34 =',cexi(19)
	    print *,'E35 =',cexi(20)
	    print *,'E36 =',cexi(21)
	    print *,'E37 =',cexi(22)
	    print *,'E38 =',cexi(23)
	    print *,'E39 =',cexi(24)
	    print *,'E310=',cexi(25)
	    print *,'E311=',cexi(26)
	    print *,'E312=',cexi(27)
	    print *,'E313=',cexi(28)
	    print *,'E314=',cexi(29)
	    print *,'E315=',cexi(30)
	    print *,'E316=',cexi(31)
	    print *,'E317=',cexi(32)
	    print *,'E318=',cexi(33)
	    print *,'E319=',cexi(34)
	    print *,'E320=',cexi(35)
*###] :
	endif
*
	if (level .eq. 3) return
*
*###] : level 3 :
	if (level .gt. 3) then
	    print *,'higher than third rank not yet implemented'
	    stop
	endif

*###] ffxex:
	end
*###[ ffeji:
	subroutine ffeji(cdxi,mdxi,ccxi,mcxi,cbxi,mbxi,caxi,maxi,
     +		cdxj,mdxj,ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
***#[*comment:***********************************************************
*									*
*	Copy the fourpoint arrays [cm][dcba]xj to the five point arrays	*
*	[cm][dcba]xi.							*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
*
	integer level
	DOUBLE PRECISION mdxi(55),mcxi(30),mbxi(30),maxi(5),
     +		mdxj(120),mcxj(140),mbxj(60),maxj(20)
	DOUBLE COMPLEX cdxi(55),ccxi(30),cbxi(30),caxi(5),
     +		cdxj(120),ccxj(140),cbxj(60),caxj(20)
*
	integer i,j,k
*
	include 'aa.h'
*
*  #] declarations:
*  #[ copy necessary parts to E0 arrays:
*
*	1)D-output: reduce the array cdxj(5*24) to cdxi(5*11)
*		    d's are calculated only to (level-1)
	do j=1,5
	    do i=1,11
		cdxi(i+(j-1)*11) = cdxj(i+(j-1)*24)
	    enddo
	enddo
*
*	2)C-output: reduce the array ccxj(20*7) to ccxi(10*3)
*		    c's are calculated only to (level-2)
*
	do j=1,4
	    do i=1,3
		ccxi(0+i+(j-1)*3) = ccxj(0+i+(j-1)*7)
	    enddo
	enddo
	do j=1,3
	    do i=1,3
		ccxi(12+i+(j-1)*3) = ccxj(35+i+(j-1)*7)
	    enddo
	enddo
	do j=1,2
	    do i=1,3
		ccxi(21+i+(j-1)*3) = ccxj(70+i+(j-1)*7)
	    enddo
	enddo
	ccxi(28) = ccxj(106)
	ccxi(29) = ccxj(107)
	ccxi(30) = ccxj(108)
*
*	check the symmetry in B0(i,j)
*	if ( atest ) then
*	    do 13 i=1,4
*		j=4+i
*		k=8+i
*		if (    ( cbxj(i)     - cbxj(i+1*12) ) .ne. 0. .or.
*     +			( cbxj(j)     - cbxj(i+2*12) ) .ne. 0. .or.
*     +			( cbxj(k)     - cbxj(i+3*12) ) .ne. 0. .or.
*     +			( cbxj(j+1*12)- cbxj(j+2*12) ) .ne. 0. .or.
*     +			( cbxj(k+1*12)- cbxj(j+3*12) ) .ne. 0. .or.
*     +			( cbxj(k+2*12)- cbxj(k+3*12) ) .ne. 0. ) then
*		    print *,'error in B0-calculations in aaxcx.for'
*		endif
* 13	    continue
*	endif
*
*  #] copy necessary parts to E0 arrays:
*###] ffeji:
	end
