c**********************************************************************
c***** This program computes the cross section for a charge-exchange
c***** reaction using the Glauber model. This is a general version 
c***** for isovector excitations with or without spin flip.
c**********************************************************************
      	implicit real*8(a-h,o-z)
      	complex*16 uim,clamb,cint(200),csumins(200),cse,cte,cso,cto,cc,
     * ctp,ctg
      	common / ra / faclog(500),dofac(500)
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
      	dimension tfp(2000),tft(2000),q(2000),rg(4),r1(200),r2(200),
     *  felast(2000),w(800),iw(102),sigma(200),ctp(4),ctg(4),
     *  par(2),cs0k(2),sumtot(200),reaction(30),rhop(200),rhot(200),
     *	tdp(200),tdt(200)
      	data par/2h -,2h +/
      	external fnor
      	external fdws
      	FACLOG(1)=0.
	FACLOG(2)=0.
      	dofac(1)=0.
      	dofac(2)=0.
      	fn=1.
      	DO 1 I=3,500
	FN=FN+1.
      	FACLOG(I)=FACLOG(I-1)+DLOG(FN)
      	dofac(i)=dofac(i-2)+dlog(fn)
   1	CONTINUE
c**********************************************************************
c****** Reaction data
c**********************************************************************
      	read(5,100)reaction
      	read(5,*)elab
      	read(5,*)a1,z1,a2,z2,a3,z3,a4,z4
      	read(5,*)aj1,aj2,aj3,aj4
      	read(5,*)ipi1,ipi2,ipi3,ipi4
      	read(5,*)t1,tz1,t2,tz2,t3,tz3,t4,tz4
      	read(5,*)signn,alfnn,af
      	read(5,*)qmin,qmax,dq
      	read(5,*)bmin,bmax,db
      	read(5,*)thmin,thmax,dth
        read(5,*)iext
      	jpi1=(ipi1+3)/2
      	jpi2=(ipi2+3)/2
      	jpi3=(ipi3+3)/2
      	jpi4=(ipi4+3)/2
      	ipit=ipi1*ipi4
      	ipip=ipi2*ipi3
c*************************************************************
c****** Calculation of the clebsh-gordan for the isospin
c*************************************************************
      	it=1
      	it1=nint(2.*t1)
      	it2=nint(2.*t2)
      	it3=nint(2.*t3)
      	it4=nint(2.*t4)
      	iz1=nint(2.*tz1)
      	iz2=nint(2.*tz2)
      	iz3=nint(2.*tz3)
      	iz4=nint(2.*tz4)
      	itt=2.*it
      	ittz=iz3-iz2
      	cleb1=cleb(it2,itt,it3,iz2,ittz,iz3)
      	cleb2=cleb(it1,itt,it4,iz1,-ittz,iz4)
      	somb=(itt+1.)**2/((it3+1.)*(it4+1.))
      	coeft=cleb1**2*cleb2**2*somb
        write(6,1100)
      	write(6,*)coeft,cleb1,cleb2
 1100   format(' coeft, cleb1, cleb2')
      	write(6,113)reaction,elab
      	write(6,101)
      	write(6,102)aj1,par(jpi1),t1,tz1,aj2,par(jpi2),t2,tz2
      	write(6,103)aj3,par(jpi3),t3,tz3,aj4,par(jpi4),t4,tz4
      	write(6,104)
      	write(6,105)signn,alfnn,af
      	write(6,111)
      	write(6,106)qmin,qmax,dq
      	write(6,107)bmin,bmax,db
      	write(6,108)thmin,thmax,dth
      	uim=cmplx(0.,1.)
      	pi=atan(1.)*4.
      	amu=931.5
      	hbarc=197.3
      	cs0k(1)=-sqrt(12.*pi)
      	cs0k(2)=sqrt(24.*pi/5.)
      	redm=a1*a2/(a1+a2)
      	ecm=elab*redm/a2
      	eta=z1*z2*1.44/hbarc*sqrt(amu*a2*0.5/elab)
      	appapt=sqrt(2.*redm*amu*ecm/hbarc**2)
      	appann=sqrt(amu*elab/(2.*a2*hbarc**2))
      	prom=(1./(2.*aj1+1.)/(2.*aj2+1.))*(appapt/appann)**2
      	write(6,*)ecm,eta,appapt,appann
      	nq=(qmax-qmin)/dq+1.001
      	nb=(bmax-bmin)/db+1.001
      	nth=(thmax-thmin)/dth+1.001
      	fnr=alfnn*signn*appann/(4.*pi)
      	fni=signn*appann/(4.*pi)
c      	write(8,*)fnr,fni
      	jtmin=abs(aj1-aj4)+.001
      	jtmax=aj1+aj4+.001      	
      	jpmin=abs(aj3-aj2)+.001
      	jpmax=aj3+aj2+.001      	
c      	write(6,*)jpmin,jpmax,jtmin,jtmax
c**********************************************************************
c****** Elastic phase-shift lambda(b)
c**********************************************************************
      	read(5,*)ndenp,ndent
c     projectile data. 
      	nn=1
      	if(ndenp.eq.2)go to 2
      	read(5,*)r0p,a0p,bep
      	write(6,*)r0p,a0p,bep
      	call fdgau(r0p,a0p,bep,tfp)
      	go to 3
   2    continue
      	read(5,*)pc,pw,pz
      	write(6,*)pc,pw,pz
      	epsabs=0.0e+00
      	epsrel=1.0e-04
      	bp=10.*pc
      	a=0.
      	ifail=0
      	call d01ajf(fnor,a,bp,epsabs,epsrel,pnor,abserr,
     *  w,800,iw,102,ifail)
      	pr0=a2/pnor
      	write(6,809)pr0
  809	format(2x,'r0=',f8.4)
      	q1=qmin-dq
      	a=0.0d0
      	do 4 iq=1,nq
      	q1=q1+dq
      	ifail=0
c      	write(6,*)pr0,pc,pw,pz
      	call d01ajf(fdws,a,bp,epsabs,epsrel,ftri,abserr,
     *  w,800,iw,102,ifail)
      	tfp(iq)=ftri
c      	write(8,*)q1,ftri
   4	continue
   3	continue
c      target data
      	nn=2
      	if(ndent.eq.2)go to 5
      	read(5,*)r0t,a0t,bet
      	write(6,*)r0t,a0t,bet
      	call fdgau(r0t,a0t,bet,tft)
      	go to 6
   5	continue
      	read(5,*)tc,tw,tz
      	write(6,*)tc,tw,tz
      	epsabs=0.0e+00
      	epsrel=1.0e-04
      	a=0.d0
      	bt=5.*tc
      	ifail=0
      	call d01ajf(fnor,a,bt,epsabs,epsrel,tnor,abserr,
     *  w,800,iw,102,ifail)
      	tr0=a1/tnor
      	write(6,809)tr0
      	q1=qmin-dq 
      	a=0.0d0
      	do 7 iq=1,nq
      	q1=q1+dq
      	ifail=0
      	call d01ajf(fdws,a,bt,epsabs,epsrel,ftri,abserr,
     *  w,800,iw,102,ifail)
      	tft(iq)=ftri
c      	write(8,*)q1,ftri
    7	continue
    6	continue
      	bs=bmin-db
      	do 8 jb=1,nb
      	bs=bs+db
      	bm=(eta+sqrt(eta**2+(appapt*bs)**2))/appapt
      	b(jb)=bs
      	qs=qmin-dq
      	do 9 jq=1,nq
      	qs=qs+dq
      	q(jq)=qs
      	qb=q(jq)*bm
      	ifail=0
      	fj0=s17aef(qb,ifail)
      	ef=exp(-0.5*af*q(jq)**2)
      	felast(jq)=tfp(jq)*tft(jq)*q(jq)*fj0*ef
c      	write(8,*)q(jq),felast(jq)
   9	continue
      	ifail=0
      	call d01gaf(q,felast,nq,ans,error,ifail)
c      	if(ans.lt.0.)write(6,*)ans,bs
c      	write(7,*)bs,ans,error
      	sigc=eta*log(appapt*bs)
      	tata=(exp(-ans*fni/appann))**2
c      	write(7,*)b(jb),tata
c      	clamb(jb)=ans*(-uim*fnr+fni)/appann
      	clamb(jb)=ans*(fnr+uim*fni)/(uim*appann)-2.*uim*sigc
   8	continue
      	do ilab=1,10
      	ilabel(ilab)=10939
      	enddo


c*************************************************************
c****** Calculation of the nucleon-nucleon transition amplitude 
c****** in terms of the Love and Franey's parametrization 
c****** in the momentum space
c*************************************************************
	qqmax=2.*appann
c      	etta=-2.*pi*hbarc**2/amu/sqrt(1.+elab/a2/2./amu)
      	etta=-4.*pi*hbarc**2/sqrt(amu**2+(hbarc*appann)**2)
      	nqq=(qqmax-qmin)/dq+1.001
      	do is=0,1
      	ikfin=2*is
      	do ik=0,ikfin,2
      	ifile=is+ik/2+1
      	ifin=3+ik/2
	do 22 i=1,ifin
     	read(ifile,*)rg(i),ser,ter,sor,tor
	read(ifile,*)sei,tei,soi,toi
	cse=ser+uim*sei
	cte=ter+uim*tei
	cso=sor+uim*soi
	cto=tor+uim*toi
     	ctp(i)=(1-is)*(cse-3.*cte-cso+3.*cto)/16.+is*(
     * ((2-ik)/2)*(-cse-cte+cso+cto)/16.-
     * (ik/2)*(-cte+cto)/4.)
       	ctg(i)=(1-is)*(cse-3.*cte+cso-3.*cto)/16.+is*(
     * ((2-ik)/2)*(-cse-cte-cso-cto)/16.)
       if(iext.eq.1)ctg(i)=ctg(i)-is*
     * (ik/2)*(-cte-cto)/4.
  22	continue
      	qq=qmin-dq
	do 23 iq=1,nqq
	qq=qq+dq
	qg=sqrt(4.*appann**2-qq**2)
      	cc=cmplx(0.,0.)
	do 24 m=1,ifin
 	factorp=((2-ik)/2)*4.*pi*rg(m)**3/(1.+(qq*rg(m))**2)+
     * (ik/2)*32.*pi*rg(m)**7*qq**2/(1.+(qq*rg(m))**2)**3
 	factorg=((2-ik)/2)*4.*pi*rg(m)**3/(1.+(qg*rg(m))**2)+
     * (ik/2)*32.*pi*rg(m)**7*qg**2/(1.+(qg*rg(m))**2)**3
	cc=cc+factorp*ctp(m)+factorg*ctg(m)
c      	write(6,*)rg(m),ctp(m),factorp,ctg(m),factorg
   24	continue
      	fnnr(ifile,iq)=real(cc)/etta
      	fnni(ifile,iq)=dimag(cc)/etta
   23	continue
	enddo
      	qq=qmin-dq
      	enddo
      	ifden=10
      	do 300 is=0,1
      	do 302 ijp=jpmin,jpmax
     	il1min=iabs(is-ijp)
      	il1max=is+ijp
      	do 303 il1=il1min,il1max
      	if(1-2*mod(il1,2).ne.ipip)go to 303
      	ifden=ifden+1
      	ilabel(ifden-10)=ijp+is*10+il1*100+1000
c****** **************************************************************
c****** If the transition densities are normalized to 1, they should be
c****** previously multiplied by A_p or A_t respectively
c****** **************************************************************
c****** Reads transition density of the projectile ***********
      	write(6,*)ifden,il1,is,ijp
      	read(ifden,'(i)')icartel
      	if(icartel.eq.ilabel(ifden-10))go to 310
      	write(6,3100)
 3100	format(' projectile density does not correspond to q.numbers')
      	go to 1000
 310    read(ifden,*)rmin2,rmax2,dr2
      	nr2=(rmax2-rmin2)/dr2+1.001
      	read(ifden,114)(rhop(ir),ir=1,nr2)
c        write(6,*)(rhop(ir),ir=1,nr2)
c****** Fourier transforms ***********************************
      	qs=qmin-dq
      	do 305 iq=1,nqq
      	qs=qs+dq
      	q(iq)=qs
        rr2=rmin2-dr2
      	do ir=1,nr2
      	rr2=rr2+dr2
      	argu=qs*rr2
      	r2(ir)=rr2
      	fj2=besself(il1,argu)
      	tdp(ir)=rhop(ir)*r2(ir)**2*fj2
      	end do
      	ifail=0
      	call d01gaf(r2,tdp,nr2,ansp,error,ifail)
      	trp(ifden-10,iq)=4.*pi*ansp
c      	write(18,*)q(iq),trp(ifden-10,iq)
  305 	continue
  303	continue
  302	continue
      	do 301 ijt=jtmin,jtmax
      	il2min=iabs(is-ijt)
      	il2max=is+ijt
      	do 304 il2=il2min,il2max
      	if(1-2*mod(il2,2).ne.ipit)go to 304
      	ifden=ifden+1
      	ilabel(ifden-10)=ijt+is*10+il2*100+5000
c****** Reads transition density of the target ***************
      	write(6,*)ifden,il2,is,ijt
      	read(ifden,*)icartel
      	if(icartel.eq.ilabel(ifden-10))go to 320
      	write(6,3200)
 3200	format(' target density does not correspond to q.numbers')
      	go to 1000
  320  	read(ifden,*)rmin1,rmax1,dr1
      	nr1=(rmax1-rmin1)/dr1+1.001
      	read(ifden,114)(rhot(ir),ir=1,nr1)
c        write(6,*)(rhot(ir),ir=1,nr1)
  998   continue
	qs=qmin-dq
      	do 306 iq=1,nqq
      	qs=qs+dq
      	q(iq)=qs
      	rr1=rmin1-dr1
      	do ir=1,nr1
      	rr1=rr1+dr1
      	argu=qs*rr1
      	r1(ir)=rr1
      	fj1=besself(il2,argu)
      	tdt(ir)=rhot(ir)*rr1**2*fj1
      	end do
      	ifail=0
      	call d01gaf(r1,tdt,nr1,anst,error,ifail)
      	trt(ifden-10,iq)=4.*pi*anst
c      write(19,*)q(iq),trt(ifden-10,iq)
 306	continue
 304	continue
 301	continue
 300	continue
c*********************************************
c********* begins sum on jt, jp, ltr, mtr
c*********************************************
      	do jth=1,nth
      	sumtot(jth)=0.d0
      		end do
      	do 30 ijt=jtmin,jtmax
      	jttw=2*ijt
c      	write(6,*)jpmin,jpmax
      	do 32 ijp=jpmin,jpmax
c      	write(6,*)ijp
      	jptw=2*ijp
      	ltrmin=iabs(ijt-ijp)
      	ltrmax=ijt+ijp
      	sombjtp=(jttw+1)*(jptw+1)
      	do 34 ltr=ltrmin,ltrmax
      	xmodltr=0.d0
c      	write(6,*)ijp
      	ltrtw=2*ltr
      	do 36 mtr=-ltr,ltr
      	if(1-2*mod(abs(mtr),2).ne.ipip*ipit)go to 36
c      	write(6,*)ijp
      	mtrtw=2*mtr
c************************************************
c********* begins summation to be squared
c************************************************
      	do 38 is=0,1
      	do jth=1,nth
      	csumins(jth)=(0.,0.)
      		end do
      	if(is.eq.1)go to 1993
c************************************************
c*******   without spin-flip
c************************************************
      	ik=0
      	if(1-2*mod(ijp,2).ne.ipip.or.1-2*mod(ijt,2).ne.ipit
     * )go to 38
      	sum=0.d0
      	do 40 m1=-ijp,ijp
      	m2=mtr-m1
      	if(iabs(m2).gt.ijt)go to 40      	      
      	term=cleb(jptw,jttw,ltrtw,2*m1,2*m2,mtrtw)*
     1	      blm(ijp,m1)*blm(ijt,m2)/sqrt(sombjtp)
      	sum=sum+term
   40	continue
      	if(abs(sum).lt.1.e-08)go to 38
      	call bint(mtr,ijp,ijt,is,ijp,ijt,ik,cint)
      	do jth=1,nth
      	csumins(jth)=sum*cint(jth)
      		end do
      	go to 1995
c***********************************************************
c*********** with spin-flip
c***********************************************************
 1993	istw=2*is
c      	write(6,*)ijt,ijp,ltr,mtr,is
      	il1min=iabs(is-ijp)
      	il1max=is+ijp
      	il2min=iabs(is-ijt)
      	il2max=is+ijt
c      	write(6,*)il1min,il1max,il2min,il2max,ipip,ipit,is,ijp,ijt
      	do 42 il1=il1min,il1max
      	if(1-2*mod(il1,2).ne.ipip)go to 42
      	l1tw=2*il1
      	do 44 il2=il2min,il2max
      	if(1-2*mod(il2,2).ne.ipit)go to 44
      	l2tw=2*il2
      	do 46 ik=0,2,2
      	ktw=2*ik
      	hatk=sqrt(ktw+1.)*cs0k(ik/2+1)*(-1.)**(il2+ik)
      	il12min=max(iabs(il1-il2),iabs(ltr-ik))
      	il12max=min(il1+il2,ltr+ik)
c      	write(6,*)ijt,ijp,ltr,mtr,is,il1,il2,ik,il12min,il12max
      	if(il12min.gt.il12max)go to 46
      	sum=0.d0
      	do 48 il12=il12min,il12max
      	l12tw=2*il12
      	anine=hatk*sqrt(l12tw+1.)
     1        *rac9(l1tw,l2tw,l12tw,istw,istw,ktw,jptw,jttw,ltrtw)
      	do 50 im=-ik,ik,2
      	if(mod(im+ik,2).ne.0)go to 50
      	m12=mtr-im
      	if(abs(m12).gt.il12)go to 50
      	mtw=2*im
      	m12tw=2*m12
      	do 52 im2=-il2,il2
      	m1=m12-im2
      	if(abs(m1).gt.il1.or.mod(m1+il1,2).ne.0.or.mod(im2+il2,2).
     * ne.0)go to 52
      	m1tw=2*m1
      	m2tw=2*im2
      	c1=cleb(l1tw,l2tw,l12tw,m1tw,m2tw,m12tw)
      	c2=cleb(l12tw,ktw,ltrtw,m12tw,mtw,mtrtw)
      	bl2=blm(il2,im2)
      	bl1=blm(il1,m1)
      	bk=blm(ik,im)
      	write(6,*)ijt,ijp,ltr,mtr,is,il1,il2,ik,il12,im,m12,im2,m1
      	write(6,*)c1,c2,bl1,bl2,bk,anine
      	termems=cleb(l1tw,l2tw,l12tw,m1tw,m2tw,m12tw)*blm(il2,im2)*anine
     1	  *cleb(l12tw,ktw,ltrtw,m12tw,mtw,mtrtw)*blm(il1,m1)*blm(ik,im)
      	sum=sum+termems
c      	if(abs(termems).lt.1.e-08)go to 52
      	write(6,*)termems,sum
   52	continue
   50	continue
   48	continue
      	if(abs(sum).lt.1.d-08)go to 46
      	call bint(mtr,il1,il2,is,ijp,ijt,ik,cint)
      	write(6,*)sum,cint(20)
      	do jth=1,nth
      	csumins(jth)=csumins(jth)+sum*cint(jth)
      		end do
   46	continue
   44	continue
   42	continue
 1995  	do jth=1,nth
      	xmod=sombjtp*cdabs(csumins(jth))**2
      	sumtot(jth)=sumtot(jth)+xmod
      		end do
      	xmodltr=xmodltr+sombjtp*cdabs(csumins(2))**2
     	write(6,*)ijp,ijt,ltr,mtr,is,xmodltr
   38	continue
   36	continue
   34	continue
   32	continue
   30	continue
        write(6,*)' cross section (mb/sr)'
      	th=thmin-dth
      	do jth=1,nth
      	th=th+dth
      	sigma(jth)=sumtot(jth)*prom*coeft*10.
      	write(6,*)th,sigma(jth)
      		end do
  21	continue
  100	format(30a1)
  113	format(2x,'reaction data',30a1,3x,'Elab=',f9.4) 
  112	format(2x,'1=target',2x,'2=projectile',2x,'3=ejectile',2x,
     *  '4=residual')
  110	format(2x,'j transf.in the target=',i3,3x,'j transf. in ', 
     *  'the projectile=',i3,3x,'spin transf=',i3,3x,'isospin transf.=',
     *  3i,3x,'tz transf.=',i3)
  101	format(2x,'initial and final quantum numbers in the form', 
     *  ' (j,parity,t,tz)')
  102	format(2x,'target:(',f3.1,a2,',',f3.1,',',f4.1,')',3x,
     *  'projectile:(',f3.1,a2,',',f3.1,',',f4.1,')')
  103	format(2x,'ejectile:(',f3.1,a2,',',f3.1,',',f4.1,')',3x,
     *  'residual:(',f3.1,a2,',',f3.1,',',f4.1,')')
  104	format(2x,'nucleon-nucleon scattering data')
  105	format(2x,'sigmann=',f8.4,2x,'alphann=',f8.4,2x,'af=',f8.4)
  111	format(2x,'step and integration limits')
  106	format(2x,'qmin=',f8.4,2x,'qmax=',f8.4,2x,'dq=',f8.4)
  107	format(2x,'bmin=',f8.4,2x,'bmax=',f8.4,2x,'db=',f8.4)
  108	format(2x,'thmin=',f8.4,2x,'thmax=',f8.4,2x,'dth=',f8.4)
  114   format(10e12.5)
 1000  	stop
      	end


c********************************************************************
      	subroutine bint(mtr,il1,il2,is,ijp,ijt,ik,cint)
      	implicit real*8(a-h,o-z)
      	complex*16 rmt,clamb,cint,uim,cfintb
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
      	dimension rmt(2000),fintbr(2000),fintbi(2000),cint(200)
c*************************************************************
c****** Integration in b
c*************************************************************
      	call cnu(mtr,il1,il2,is,ijp,ijt,ik,rmt)
      	write(6,*)mtr,il1,il2,is,ijp,ijt,ik,rmt(100)
      	th=thmin-dth
      	do 21 jth=1,nth
      	th=th+dth
      	thr=th*pi/180.
      	do 19 jb=1,nb
      	arg=2.*appapt*b(jb)*sin(thr/2.)
      	fjmt=besseln(abs(mtr),arg)
      	if(mtr.lt.0)fjmt=fjmt*(-1)**mtr
      	cfintb=b(jb)*rmt(jb)*fjmt*cdexp(-clamb(jb))
      	fintbr(jb)=dreal(cfintb)
      	fintbi(jb)=dimag(cfintb)
c      	write(8,*)b(jb),rmt(jb)
c    	write(9,*)b(jb),fintbi(jb)
  19 	continue
      	ifail=0
      	call d01gaf(b,fintbr,nb,ansbr,errorr,ifail)
      	ifail=0
      	call d01gaf(b,fintbi,nb,ansbi,errori,ifail)
      	cint(jth)=cmplx(ansbr,ansbi)
c      	write(8,*)ansbr,errorr,ifail
c      	write(9,*)ansbi,errori,ifail
  21	continue
      	return
      	end


c*************************************************************
      	subroutine cnu(mt,il1,il2,is,ijp,ijt,ik,rmt)
      	implicit real*8(a-h,o-z)
      	complex*16 uim,rmt,cse,cte,cso,cto,cc,ctp,ctg,ct,clamb
      	common / ra / faclog(500),dofac(500)
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
      	dimension q(2000),fintqr(2000),fintqi(2000),rmt(2000)
c*************************************************************
c****** integration in dq
c*************************************************************
c      	write(18,*)mt,il1,il2,is,ijp,ijt,ik
      	isk=is+(ik/2)+1
      	icartel1=100*il1+10*is+ijp+1000
      	icartel2=100*il2+10*is+ijt+5000
      	do 5 ilab=1,10
      	if(icartel1.ne.ilabel(ilab))go to 3
      	ilabel1=ilab
c      	write(18,*)icartel1,ilabel1,ilabel(ilabel1)
  3     if(icartel2.ne.ilabel(ilab))go to 5
      	ilabel2=ilab
c      	write(18,*)icartel2,ilabel2,ilabel(ilabel2)
  5	continue
      	nqq=(qqmax-qmin)/dq+1.001
      	bs=bmin-db
      	do 17 ib=1,nb
      	bs=bs+db
      	bm=(eta+sqrt(eta**2+(appapt*bs)**2))/appapt
      	qs=qmin-dq
      	do 16 iq=1,nqq
      	qs=qs+dq
      	q(iq)=qs
      	argu=bm*q(iq)
      	fj=besseln(abs(mt),argu)
      	if(mt.lt.0)fj=fj*(-1)**mt
      	fintq=trp(ilabel1,iq)*trt(ilabel2,iq)*q(iq)
c      	if(trp(ilabel1,iq).eq.0.or.trt(ilabel2,iq).eq.0.)write(18,*)
c     * ilabel1,ilabel2
      	ftqr=fintq*fnnr(isk,iq)
      	ftqi=fintq*fnni(isk,iq)
      	fintqr(iq)=ftqr*fj
      	fintqi(iq)=ftqi*fj
   16	continue
      	ifail=0
      	call d01gaf(q,fintqr,nqq,ansqr,errorqr,ifail)
      	ifail=0
      	call d01gaf(q,fintqi,nqq,ansqi,errorqi,ifail)
      	rmt(ib)=ansqr+uim*ansqi
c      	write(25,*)bs,ansqr,errorqr
c      	write(26,*)bs,ansqi,errorqi
   17    continue
 	return
      	end


C*******************************************************************
        function blm(l,m)
      	implicit real*8(a-h,o-z)
      	common / ra / faclog(500),dofac(500)
      	pi=atan(1.)*4.
      	blm=0.d0
c*************************************************************
c****** Blm geometrical factors (without any phase factor)
c*************************************************************
      	sq1=dsqrt((2.*l+1.)/(4.*pi))
      	if(mod((l+m),2).ne.0)go to 10
      	fact=exp(faclog(l+m+1)+faclog(l-m+1))
      	dfact=exp(dofac(l+m+1)+dofac(l-m+1))
      	sq2=sqrt(fact)
      	blm=sq1*sq2/dfact
c      	write(6,*)l,m,blm
  10	continue
        return
      	end


      	subroutine fdgau(r0,a0,be,tf)
c********************************************************************
c****** subroutine for the calculation of the density in the momentum
c****** space using a gaussian shape for the density
c********************************************************************
      	implicit real*8(a-h,o-z)
      	complex*16 uim,clamb
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
      	dimension tf(2000)
      	q=qmin-dq
      	do 1 i=1,nq
      	q=q+dq
      	tf(i)=(pi*a0**2)**1.5*r0*exp((-1)*(q*a0)**2/4.)
c      	write(8,*)q,tf(i)
      	if(be.eq.0.)go to 1
      	tf(i)=tf(i)+pi**1.5*be*r0*a0**5*exp((-0.25)*(q*a0)**2)*
     1  (1.5-(q*a0)**2/4.)
    1	continue
      	return
      	end


      	function fnor(x)
c********************************************************************
c****** subroutine for the calculation of the normalization of the
c****** density in the case of Fermi distrib. Normalizes to the mass.
c********************************************************************
      	implicit real*8(a-h,o-z)
      	complex*16 uim,clamb
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
      	if(nn.eq.1)go to 1
      	w=tw
      	c=tc
      	z=tz
      	go to 2
   1	continue
      	w=pw
      	c=pc
      	z=pz
   2	continue
      	den=(1+w*(x/c)**2)/(1.+exp((x-c)/z))
      	fnor=den*pi*4.*x**2
      	return
      	end

      	function fdws(x)
c********************************************************************
c****** subroutine for the calculation of the density in the momentum
c****** space using a Fermi shape for the density
c********************************************************************
      	implicit real*8(a-h,o-z)
      	complex*16 uim,clamb
      	common pw,pc,pz,tw,tc,tz,r0p,a0p,bep,r0t,a0t,bet,clamb(2000),
     * b(2000),pdel,tdel,nq,nb,nth,qmin,qmax,dq,bmin,bmax,db,thmin,
     * thmax,dth,appapt,alfnn,eta,pi,q1,pr0,tr0,nn,uim,appann,hbarc,
     * elab,amu,fnnr(3,2000),fnni(3,2000),trp(10,2000),trt(10,2000),
     * ilabel(10),qqmax
c      	write(6,*)pr0,pc,pw,pz
      	if(nn.eq.1)go to 1
      	r0=tr0
      	w=tw
      	c=tc
      	z=tz
      	go to 2
   1	continue
      	r0=pr0
      	w=pw
      	c=pc
      	z=pz
   2	continue
      	if(q1.eq.0.)go to 25
      	fac=4.*pi*r0*x*sin(x*q1)/q1
      	go to 26
  25	continue
      	fac=4.*pi*r0*x**2
  26	continue
      	fac1=exp((x-c)/z)
      	fdws=fac/(1.+fac1)
c      	write(6,*)r0,c,w,z,pi
      	if(w.eq.0)go to 3
      	fdws=fdws+fac*w*(x/c)**2/(1.+fac1)
   3	continue
c      	write(8,*)q1,x,fdws
      	return
      	end

      	function besseln(n,argu)
c********************************************************************
c****** subroutine for the calculation of Bessel functions of 
c****** integer order
c********************************************************************
      	implicit real*8(a-h,o-z)
      	dimension fj(0:10)
      	ifail=0
      	fj(0)=s17aef(argu,ifail)
      	ifail=0
      	fj(1)=s17aff(argu,ifail)
      	if(n.lt.2)go to 2
      	if(argu.eq.0.)go to 25
      	do 1 i=2,n
      	fj(i)=2.*(i-1)/argu*fj(i-1)-fj(i-2)
    1	continue
      	go to 2
   25	continue
      	fj(n)=0.0d0
    2	continue
      	besseln=fj(n)
c      	write(7,*)n,argu,besseln
      	return
      	end

      	function besself(n,argu)
c********************************************************************
c****** subroutine for the calculation of Bessel functions of 
c****** fractional order
c********************************************************************
      	implicit real*8(a-h,o-z)
      	dimension fj(0:10)
      	ifail=0
      	if(argu.eq.0.)go to 25
      	fj(0)=sin(argu)/argu
      	fj(1)=sin(argu)/argu**2-cos(argu)/argu
      	if(n.lt.2)go to 2
      	do 1 i=2,n
      	fj(i)=(2.*i-1.)/argu*fj(i-1)-fj(i-2)
    1	continue
      	go to 2
   25	continue
      	fj(0)=1.
      	if(n.lt.1)go to 2
      	fj(n)=0.0d0
    2	continue
      	besself=fj(n)
c      	write(7,*)n,argu,besseln
      	return
      	end