	!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
	!x  B+ -> Rho^0 K^+ Q2-body  | W.F. Wang  2016-07-06 x!
	!x     -> K^+ pi^+pi^-       |    Updated 2018-07-10 x!
	!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx!
	!                                                     !
	subroutine setcon    
	implicit double precision(a-h,o-z)    
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
	common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h   !for h=kaon			
	common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta    !for Resonance
	common/ctmp/xyz,xmc,xtmp,xvv,xqq,xoo,xss,xhh,x00    
	xlamb=0.250   
	xlam =0.193
	xmw  =80.4
	pi   =3.14159265359
	xmB  =5.27932                 !Mass of Bu meson PDG
	vmbq =4.8                     !Mass of b-quark
	fB   =0.190                   ! DC of B/Bs meson 
	wB   =0.4 !-0.04              !Omega_B0
!	xNB  =0.0                     !Auto-generated
	ckA = 0.836                   !CKM A
	ckl = 0.22453                 !CKM lambda  <=PDG
	ckr = (0.122)/(1.-ckl*ckl/2.) !CKM rho
!	cke =-(0.355)/(1.-ckl*ckl/2.) !CKM eta 
	xmh = 0.493677  !Mass for h=kaon
	xm0h= 1.90      !m0^h for h
	fh  = 0.160     !D.C. for h 
      a1h =-0.15       !Gen moms
      a2h = 0.25
      a4h =-0.015
      rhoh= (xmh/xm0h)
      et3h= 0.015
      w3h =-3.0 
	xmR = 0.77526      !Mass 
	xGaR= 149.1/1000.  !Full 
	fR  = 0.216        !for R
c	fTR = 0.           !Ratio needed
	a1R = 0.00
	a2R = 0.25
	a4R = 0.75
	xeta=-0.60
c	eta = 0.0       
	xyz = 0.0
!	xmc = 
!	xtmp=
	xvv =-1.0
!	xqq =
!	xoo =
!	xss =
!	xhh =
!	x00 =  
	END
	
	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX    
	program main
	implicit double precision(a-h,o-z)    
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(20),XU(20),ACC !for VEGAS
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
	common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h   !			
	common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta    !
	common/ctmp/xyz,xmc,xtmp,xvv,xqq,xoo,xss,xhh,x00    !
	common/rsut/xzbr,xzcp,xztem
	character*8::now  !Time now
	integer::i
	external funNB,DAtot
	!--------------
	call setcon	  !Call constants	
	!--------------NB
	ncall=20000
	call vegas(1,funNB,vnb,sd,chi2a)
	xNB=fb/2./sqrt(2.*3.)/vnb ! N_B for B meson
!	write(*,*) xNB
	!--------------NB	
	!--------------xu
	do i=8,13 
	xu(i) =3.9999999   !xu(8..13)
      end do  
	!--------------xu

	write(*,*)
	write(*,*)  ' Q2B  B+ = rho0 K+       '
	write(*,*)  

	xu(7)=0.322914  !XU#  0.322914->3GeV  
	xl(7)=0.0028    !#XL  0.0028  ->m_{pi+pi}

	call time(now); write(*,*) '--Now-- ',now 
	NCALL=50000
	call vegas(13,DAtot,vbr,sd,chi2a)  !Time test
	!-------------

	NCALL=250000;  call funmain
	NCALL=500000;  call funmain

	end                                      !---- END of main ----!	 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX	
	subroutine funmain
	implicit double precision(a-h,o-z)
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(20),XU(20),ACC !for VEGAS
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
	common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h   			
	common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta   
	common/ctmp/xyz,xmc,xtmp,xvv,xqq,xoo,xss,xhh,x00   
	common/rsut/xzbr,xzcp,xztem
	character*8::now
	external DAtot

	call time(now); write(*,*) '--Now-- ',now 
	
	cke =-(0.355)/(1.-ckl*ckl/2.) 
	call vegas(13,DAtot,vbr1,sd,chi2a)
	write(*,*) ' vbr1:', vbr1
	!----
	cke = (0.355)/(1.-ckl*ckl/2.)
	call vegas(13,DAtot,vbr2,sd,chi2a)
	write(*,*) ' vbr2:', vbr2

	xzbr=(vbr1+vbr2)/2.  
	xzcp=(vbr1-vbr2)/(vbr1+vbr2)

	write(*,*)  '    '
	write(*,*)  '    [Q2B]',xzbr
	write(*,*)  '     Acp', xzcp

	END
	!&&&&&&&&&&&&&&&&&&&&&&&&&&&&!
	double precision function DAtot(x,wgt)  ! for all
   	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
	common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h  			
	common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta   
      dimension x(13)
	double precision v1ab(2),v1abp(2),v1cd(2),v1cdp(2)
	double precision v1ef(2),v1efp(2),v1gh(2),v1ghp(2)
	double precision v1mn(2),v1mnp(2),v1op(2),v1opp(2)
	
	xB=x(1);z =x(2);x3=x(3); xBp=x(4); zp =x(5); x3p=x(6)
	bB=x(8);b =x(9);b3=x(10);bBp=x(11);bp =x(12);b3p=x(13)
		
	et= x(7)    ! <--------------<< eta  !

	call Figab(xB, z, x3, bB, b, b3, et,v1ab )
	call Figab(xBp,zp,x3p,bBp,bp,b3p,et,v1abp)
	call Figcd(xB, z, x3, bB, b, b3, et,v1cd )
	call Figcd(xBp,zp,x3p,bBp,bp,b3p,et,v1cdp)
	call Figef(xB, z, x3, bB, b, b3, et,v1ef )
	call Figef(xBp,zp,x3p,bBp,bp,b3p,et,v1efp)
	call Figgh(xB, z, x3, bB, b, b3, et,v1gh )
	call Figgh(xBp,zp,x3p,bBp,bp,b3p,et,v1ghp)
	call Figmn(xB, z, x3, bB, b, b3, et,v1mn )
	call Figmn(xBp,zp,x3p,bBp,bp,b3p,et,v1mnp)
	call Figop(xB, z, x3, bB, b, b3, et,v1op )
	call Figop(xBp,zp,x3p,bBp,bp,b3p,et,v1opp)

	!------------ Re terms
	AAA= v1ab(1)+ v1cd(1)+ v1ef(1)+ v1gh(1)+ v1mn(1)+ v1op(1)
	BBB=v1abp(1)+v1cdp(1)+v1efp(1)+v1ghp(1)+v1mnp(1)+v1opp(1)
	!------------ Im terms
	CCC= v1ab(2)+ v1cd(2)+ v1ef(2)+ v1gh(2)+ v1mn(2)+ v1op(2)
	DDD=v1abp(2)+v1cdp(2)+v1efp(2)+v1ghp(2)+v1mnp(2)+v1opp(2)

	!------------
	xmka=0.1396 !here for pion
	xw  =xmB*dsqrt(et)  ! 
	xw2 =xw*xw          ! xw2=s
	xmB2=xmB*xmB        ! mB^2

	fac1=dsqrt(xw2-4.* xmka**2.)/2.                           !=|p_1|=|p_2|
	fac2=((xmB2-(xw+xmh)**2.)*(xmB2-(xw-xmh)**2.))**0.5/2./xw !=|p_3|
	!GF*GF/2.*1.5192675E12=103.34..  !Bu+ 1.638  PDG	
	vcoef=1.638*103.34356*(fac1*fac2)**3./    
     &      (12*pi**3.*xmB*(xmB2)**2.)                       !Br=vcoef*|A|^2
	DAtot=vcoef*(FFGS(xw)**2.)*(AAA*BBB+CCC*DDD)*xlamb**2. !notice xlamb^2     
	return
      end
	!mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm Form factor 
      double precision function FFGS(xw)
      implicit double precision(a-z)
	COMPLEX::Gs1,funw,Gsp,cp,Gs2p,c2p,Gs3p,c3p

	s=xw*xw
	sqrs=Sqrt(s)
	sqfs=Sqrt(-0.07792+s)
	lgs=Log(3.58243*(sqfs+sqrs))
	!----------------
	Gs1=(2.0618635468*sqrs)/(cmplx(0,0.057941)*sqfs  ! <-- for rho770
     &    +2.09875*sqrs-cmplx(0,0.7436)*sqfs*s
     &    -4.1379*s**1.5+(-0.03688646+0.473394*s)*sqfs*lgs)
	!----------------
      Gsp=(7.768275*sqrs)/(cmplx(0,0.07385)*sqfs+7.8153*sqrs - 
     &    cmplx(0,0.947758)*sqfs*s-4.8566763*s**1.5+ 
     &    (-0.047+0.60336143*s)*sqfs*lgs)	             ! <-- for rho1450

	FFGS=abs(Gs1)*sqrt(2.)  !notice the sqrt(2.)
	return
	end
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
	
czzzzzzzzzzzzzzzzzzzzzzzzzzzzz WAVE Funs zz Begin zzzzzz	
      subroutine DAh(x,vah)  ! h here  for Kaon    
	implicit double precision(a-h,o-z)  
 	common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h  	
	double precision vah(3) 
	vt =(2.*x-1);  vt2=vt*vt
			
	vah(1)=fh*6.*x*(1.-x)*(1. +a1h*3.*vt +a2h*(5.*vt2-1.)*1.5           	
     &       +a4h*15./8.*(1.-14.0*vt2+21.*vt**4.) )*0.204124    	    
      vah(2)=fh*(1.+(30.*et3h-2.5*rhoh**2.)*(3.*vt2-1.)/2.
     &      -3.*(et3h*w3h+rhoh**2.*9.0*(1.+6.*a2h)/20. )            
     &       *(3.-30.*vt2+35.*vt**4.)/8. )*0.204124   	
      vah(3)=fh*(-vt)*(1. +6.*(5.*et3h-0.5*et3h*w3h-7.*rhoh**2./20.
     &       -0.6*rhoh**2.*a2h)*(1.-10.*x+10.*x**2.))*0.204124

	return
      end
	!zzzzzzzzzzzzz-------------
      subroutine DAR(y,vpR)  ! R rho770 here 
	implicit double precision(a-h,o-z)  
	common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta 
	double precision vpR(3) !phi0z=vpR(1); phisz=vpR(2); phitz=vpR(3)
	vy =(2.*y-1);  vy2=vy*vy	
     
	vpR(1)=6.*y*(1.-y)*(1.-a1R*(3.*vy)+a2R*1.5*(5.*vy2-1.))
     &                                               *0.204124
!	vpR(2)=fTR/fR*3.*(-vy)*0.204124
!	vpR(3)=fTR/fR*3.*(vy2)*0.204124 	
	vpR(2)=0.852*3.*0.204124*(-vy)*(1.+a4R*(1.-10.*y+10.*y*y))  !fTR/fR=0.852
	vpR(3)=0.852*3.*0.204124*(vy2+xeta*1.5*vy2*(5.*vy2-1.))          	

	return                                   
      end	           
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzz WAVE Funs zz END zzzzzz   

c############################# Decay Amplitudes ## Begin #
      subroutine Figab(xB,z,x3,bB,b,b3,et,v1ab)      ! Fig1-ab
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1ab(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=z
      be=xB*z
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-a
      t1a=MAX(tmp,1./b,1./bB)
      call hfun(al,be,b,bB,vfun)
      har=vfun(1)*1.77493*(z*(1.-z))**0.3
      hai=vfun(2)*1.77493*(z*(1.-z))**0.3

      al=xB   !-et+
      be=xB*z
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-b
      t1b=MAX(tmp,1./bB,1./b)
      call hfun(al,be,bB,b,vfun)
      hbr=vfun(1)*1.77493*(xB*(1.-xB))**0.3
      hbi=vfun(2)*1.77493*(xB*(1.-xB))**0.3

      sud=suda(xB,bB,xmB)+suda(z,b,xmB)+suda((1.-z),b,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp(5./3.*ss(bB,t1a)+2.*ss(b,t1a))
      su2=ssu*dexp(5./3.*ss(bB,t1b)+2.*ss(b,t1b))

      phiBx=phib(xB,bB)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------
      ctmp=4./3*8*(-1.+et)*fh*phiBx*Pi*xmB**4*bB*b
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=phi0z*(1+z)-Sqrt(et)*(phisz+phitz)*(-1+2*z)
      factb=et*phi0z+2*Sqrt(et)*phisz

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------Fsp---------
      ctmp=4./3*16*fh*phiBx*Pi*r0*xmB**4*bB*b
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=Sqrt(et)*(-(phitz*z)+phisz*(2+z))+phi0z*(1+et*(-1+2*z))
      factb=2*Sqrt(et)*phisz*(1+et-xB)+phi0z*xB

      Fspar=coea*facta*har
      Fspai=coea*facta*hai
      Fspbr=coeb*factb*hbr
      Fspbi=coeb*factb*hbi
      !------------Fsp------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(1)/3+wc(2)
      wcay=wc(3)/3+wc(4)+wc(9)/3+wc(10)
      wcaz=wc(5)/3+wc(6)+wc(7)/3+wc(8)
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(1)/3+wc(2)
      wcby=wc(3)/3+wc(4)+wc(9)/3+wc(10)
      wcbz=wc(5)/3+wc(6)+wc(7)/3+wc(8)
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FSPar-vti*FSPai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FSPai+vti*FSPar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FSPbr-vti*FSPbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FSPbi+vti*FSPbr)

	v1ab(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1ab(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb  
      return
      end
      !------------------------
      subroutine Figcd(xB,z,x3,bB,b,b3,et,v1cd)      ! Fig1-cd
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1cd(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=xB*z
      be=(-1+et+x3-et*x3+xB)*z
      !be=(-1.+x3+xB)*z
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-c
      t1a=MAX(tmp,1./bB,1./b3)
      call hfun(al,be,bB,b3,vfun)
      har=vfun(1)
      hai=vfun(2)

      al=xB*z
      be=((-1+et)*x3+xB)*z
      !be=(-x3+xB)*z
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-d
      t1b=MAX(tmp,1./bB,1./b3)
      call hfun(al,be,bB,b3,vfun)
      hbr=vfun(1)
      hbi=vfun(2)

      sud=suda(xB,bB,xmB)+suda(z,bB,xmB)+suda((1.-z),bB,xmB)
     &    +suda(x3,b3,xmB)+suda((1.-x3),b3,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp((5./3.+2.)*ss(bB,t1a)+2.*ss(b3,t1a))
      su2=ssu*dexp((5./3.+2.)*ss(bB,t1b)+2.*ss(b3,t1b))

      phiBx=phib(xB,bB)
      call DAh(x3,vah) !for h
      phAx =vah(1)
      phPx =vah(2)
      phTx =vah(3)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------
      ctmp=4./3*16*Sqrt(2./3)*(-1.+et)*phAx*phiBx*Pi*xmB**4*bB*b3
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=phi0z*(1-x3-xB+et*(-1+x3-z))+Sqrt(et)*(-phisz+phitz)*z
      factb=(-1+et)*phi0z*x3+phi0z*xB+(-phi0z+Sqrt(et)*(phisz+phitz))
     &      *z

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------Flr---------
      ctmp=4./3*16*Sqrt(2./3)*phiBx*Pi*r0*xmB**4*bB*b3
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=phi0z*phPx*(1-x3-xB+et*(-1+x3-z))+phi0z*phTx*(1-x3-xB+et
     &      *(-1+x3+z))+Sqrt(et)*(phTx*(phisz*(1+et*(-1+x3)-x3-xB-z)
     &      +phitz*(-1+et+x3-et*x3+xB-z))+phPx*(phisz*(1+et*(-1+x3)
     &      -x3-xB+z)+phitz*(-1+et+x3-et*x3+xB+z)))
      factb=(-1+et)*(phi0z+Sqrt(et)*(phisz-phitz))*(phPx-phTx)*x3
     &      +(phi0z+Sqrt(et)*(phisz-phitz))*(phPx-phTx)*xB+Sqrt(et)
     &      *(Sqrt(et)*phi0z-phisz-phitz)*(phPx+phTx)*z

      Flrar=coea*facta*har
      Flrai=coea*facta*hai
      Flrbr=coeb*factb*hbr
      Flrbi=coeb*factb*hbi
      !------------Flr------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(1)
      wcay=wc(3)+wc(9)
      wcaz=wc(5)+wc(7)
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(1)
      wcby=wc(3)+wc(9)
      wcbz=wc(5)+wc(7)
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FLRar-vti*FLRai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FLRai+vti*FLRar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FLRbr-vti*FLRbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FLRbi+vti*FLRbr)

	v1cd(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1cd(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb   
      return
      end
      !------------------------
      subroutine Figef(xB,z,x3,bB,b,b3,et,v1ef)      ! Fig1-ef
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1ef(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=-1.+z
      be=-((et*(-1+x3)-x3)*(-1+z))
      !be=x3*(-1.+z)
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-e
      t1a=MAX(tmp,1./b,1./b3)
      call hfun(al,be,b,b3,vfun)
      har=vfun(1)*1.77493*(z*(1.-z))**0.3
      hai=vfun(2)*1.77493*(z*(1.-z))**0.3

      al=et*(-1+x3)-x3
      be=-((et*(-1+x3)-x3)*(-1+z))
      !al=-x3
      !be=x3*(-1.+z)
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-f
      t1b=MAX(tmp,1./b3,1./b)
      call hfun(al,be,b3,b,vfun)
      hbr=vfun(1)*1.77493*(x3*(1.-x3))**0.3
      hbi=vfun(2)*1.77493*(x3*(1.-x3))**0.3

      sud=suda(z,b,xmB)+suda((1.-z),b,xmB)
     &    +suda(x3,b3,xmB)+suda((1.-x3),b3,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp(2.*ss(b,t1a)+2.*ss(b3,t1a))
      su2=ssu*dexp(2.*ss(b,t1b)+2.*ss(b3,t1b))

      phiBx=phib(xB,bB)
      call DAh(x3,vah) !for h
      phAx =vah(1)
      phPx =vah(2)
      phTx =vah(3)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------
      ctmp=4./3*8*fB*Pi*xmB**4*b3*b
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=(-1.+et)*phAx*phi0z*(-1+z)+2*Sqrt(et)*phPx*r0*(phisz*(-2
     &      +z)-phitz*z)
      factb=-(et**2*phAx*phi0z*(-1+x3))-2*et**1.5*phisz*(phPx+phTx)
     &      *r0*(-1+x3)-phAx*phi0z*x3+et*phAx*phi0z*(-1+2*x3)+2
     &      *Sqrt(et)*phisz*r0*(phTx*(-1+x3)+phPx*(1+x3))

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------Fsp---------
      ctmp=4./3*16*fB*Pi*xmB**4*b3*b 
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=2*phi0z*phPx*r0*(1+et*(-1+z))-(-1+et)*Sqrt(et)*phAx
     &      *(phisz+phitz)*(-1+z)
      factb=2*(-1+et)*Sqrt(et)*phAx*phisz-phi0z*r0*(2*et*phTx+(-1+et)
     &      *(phPx-phTx)*x3)

      Fspar=coea*facta*har
      Fspai=coea*facta*hai
      Fspbr=coeb*factb*hbr
      Fspbi=coeb*factb*hbi
      !------------Fsp------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(1)/3+wc(2)
      wcay=wc(3)/3+wc(4)+wc(9)/3+wc(10)
      wcaz=wc(5)/3+wc(6)+wc(7)/3+wc(8)
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(1)/3+wc(2)
      wcby=wc(3)/3+wc(4)+wc(9)/3+wc(10)
      wcbz=wc(5)/3+wc(6)+wc(7)/3+wc(8)
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FSPar-vti*FSPai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FSPai+vti*FSPar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FSPbr-vti*FSPbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FSPbi+vti*FSPbr)

	v1ef(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1ef(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb  
      return
      end
      !------------------------
      subroutine Figgh(xB,z,x3,bB,b,b3,et,v1gh)      ! Fig1-gh
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1gh(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=-((et*(-1+x3)-x3)*(-1+z))
      be=1+(-1+et+x3-et*x3+xB)*z
      !al=x3*(-1.+z)
      !be=1.+(-1.+x3+xB)*z
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-g
      t1a=MAX(tmp,1./b,1./bB)
      call hfun(al,be,b,bB,vfun)
      har=vfun(1)
      hai=vfun(2)

      al=-((et*(-1+x3)-x3)*(-1+z))
      be=-((et*(-1+x3)-x3+xB)*(-1+z))
      !al=x3*(-1.+z)
      !be=(x3-xB)*(-1.+z)
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-h
      t1b=MAX(tmp,1./b,1./bB)
      call hfun(al,be,b,bB,vfun)
      hbr=vfun(1)
      hbi=vfun(2)

      sud=suda(xB,bB,xmB)+suda(z,b,xmB)+suda((1.-z),b,xmB)
     &    +suda(x3,b,xmB)+suda((1.-x3),b,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp((5./3.)*ss(bB,t1a)+4.*ss(b,t1a))
      su2=ssu*dexp((5./3.)*ss(bB,t1b)+4.*ss(b,t1b))

      phiBx=phib(xB,bB)
      call DAh(x3,vah) !for h
      phAx =vah(1)
      phPx =vah(2)
      phTx =vah(3)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------
      ctmp=4./3*16*Sqrt(2./3)*phiBx*Pi*xmB**4*b*bB
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=et**1.5*(phisz-phitz)*(phPx+phTx)*r0*(-1+x3)+phAx*phi0z
     &      *(x3+xB)+et**2*phAx*phi0z*(-1+x3-z)+et*phAx*phi0z*(1-2*x3
     &      -xB+z)+Sqrt(et)*r0*(-(phisz*phPx*(3+x3+xB-z))+phitz*phPx
     &      *(-1+x3+xB+z)-phTx*(phitz-phitz*(x3+xB-z)+phisz*(-1+x3+xB
     &      +z)))
      factb=(-1+et)**2*phAx*phi0z*(-1+z)+Sqrt(et)*r0*(phisz*phPx*(1
     &      +et+x3-et*x3-xB-z)+phitz*phPx*(-1+et+x3-et*x3-xB+z)+phTx
     &      *(phisz*(1+et*(-1+x3)-x3+xB-z)+phitz*(-1+et*(-1+x3)-x3+xB
     &      +z)))

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------Flr---------
      ctmp=4./3*16*Sqrt(2./3)*phiBx*Pi*xmB**4*b*bB
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=-((-1+et)*Sqrt(et)*phAx*(phisz-phitz)*(1+z))-phi0z*r0
     &      *(phTx*(-2+x3-et*x3+xB-et*z)+phPx*(-2+x3+xB+et*(2-x3+z)))
      factb=(-1+et)*Sqrt(et)*phAx*(phisz-phitz)*(-1+z)-phi0z*r0*((-1
     &      +et)*(phPx+phTx)*x3+(phPx+phTx)*xB+et*phTx*(-2+z)-et*phPx
     &      *z)

      Flrar=coea*facta*har
      Flrai=coea*facta*hai
      Flrbr=coeb*factb*hbr
      Flrbi=coeb*factb*hbi
      !------------Flr------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(1)
      wcay=wc(3)+wc(9)
      wcaz=wc(5)+wc(7)
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(1)
      wcby=wc(3)+wc(9)
      wcbz=wc(5)+wc(7)
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FLRar-vti*FLRai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FLRai+vti*FLRar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FLRbr-vti*FLRbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FLRbi+vti*FLRbr)

	v1gh(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1gh(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb 
      return
      end
      !------------------------
      subroutine Figmn(xB,z,x3,bB,b,b3,et,v1mn)      ! Fig1-mn
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1mn(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=x3-et*x3
      be=-((-1+et)*x3*xB)
      !al=x3
      !be=x3*xB
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-m
      t1a=MAX(tmp,1./b3,1./bB)
      call hfun(al,be,b3,bB,vfun)
      har=vfun(1)*1.77493*(x3*(1.-x3))**0.3
      hai=vfun(2)*1.77493*(x3*(1.-x3))**0.3

      al=xB-et*xB
      be=-((-1+et)*x3*xB)
      !al=xB
      !be=x3*xB
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-n
      t1b=MAX(tmp,1./bB,1./b3)
      call hfun(al,be,bB,b3,vfun)
      hbr=vfun(1)*1.77493*(xB*(1.-xB))**0.3
      hbi=vfun(2)*1.77493*(xB*(1.-xB))**0.3

      sud=suda(xB,bB,xmB)+suda(x3,b3,xmB)+suda((1.-x3),b3,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp(5./3.*ss(bB,t1a)+2.*ss(b3,t1a))
      su2=ssu*dexp(5./3.*ss(bB,t1b)+2.*ss(b3,t1b))

      phiBx=phib(xB,bB)
      call DAh(x3,vah) !for h
      phAx =vah(1)
      phPx =vah(2)
      phTx =vah(3)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------fR*
      ctmp=4./3*8*phiBx*Pi*xmB**4*bB*b3  
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=-(phTx*r0*(1+et+2*(-1+et)*x3))-(-1+et)*(phPx*r0*(-1+2*x3)
     &      +phAx*(-1+(-1+et)*x3))
      factb=-2*phPx*r0*(1+et*(-1+xB))-(-1+et)*et*phAx*xB

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------FLL---------
	!FLL=FLR
      FLRar=FLLar
      FLRai=FLLai
      FLRbr=FLLbr
      FLRbi=FLLbi
      !------------FLL------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(1)+wc(2)/3                 !notice pi0=(uu-dd)/sqrt2
      wcayA=wc(3)+wc(4)/3+wc(9)+wc(10)/3 !for uu_bar component
      wcazA=wc(5)+wc(6)/3+wc(7)+wc(8)/3  !for uu_bar component
      wcay=wcayA -( wc(3)+wc(4)/3-0.5*wc(9)-0.5*wc(10)/3 )  ! for dd_bar
      wcaz=wcazA -( wc(5)+wc(6)/3-0.5*wc(7)-0.5*wc(8)/3 )   ! for dd_bar
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(1)+wc(2)/3
      wcbyB=wc(3)+wc(4)/3+wc(9)+wc(10)/3
      wcbzB=wc(5)+wc(6)/3+wc(7)+wc(8)/3
      wcby=wcbyB -( wc(3)+wc(4)/3-0.5*wc(9)-0.5*wc(10)/3 )
      wcbz=wcbzB -( wc(5)+wc(6)/3-0.5*wc(7)-0.5*wc(8)/3 )
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FLRar-vti*FLRai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FLRai+vti*FLRar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FLRbr-vti*FLRbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FLRbi+vti*FLRbr)

	v1mn(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1mn(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb 
      return
      end
      !------------------------
      subroutine Figop(xB,z,x3,bB,b,b3,et,v1op)      ! Fig1-op
      implicit double precision(a-h,o-z)
      common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB,ckA,ckl,ckr,cke
      common/conh/xmh,xm0h,fh,a1h,a2h,a4h,rhoh,et3h,w3h
      common/conR/xmR,xGaR,fR,fTR,a1R,a2R,a4R,xeta,eta
      double precision wc(10),v1op(2),vfun(2),vah(3),vpR(3)
      r0=xm0h/xmB

      al=-((-1+et)*x3*xB)
      be=-((et*(-1+x3)-x3)*(-1+xB+z))
      !al=x3*xB
      !be=x3*(-1.+xB+z)
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-o
      t1a=MAX(tmp,1./bB,1./b)
      call hfun(al,be,bB,b,vfun)
      har=vfun(1)
      hai=vfun(2)

      al=-((-1+et)*x3*xB)
      be=-((-1+et)*x3*(xB-z))
      !al=x3*xB
      !be=x3*(xB-z)
      tmp=MAX(xmB*dsqrt(dabs(al)),xmB*dsqrt(dabs(be))) !fig-p
      t1b=MAX(tmp,1./bB,1./b)
      call hfun(al,be,bB,b,vfun)
      hbr=vfun(1)
      hbi=vfun(2)

      sud=suda(xB,bB,xmB)+suda(z,b,xmB)+suda((1.-z),b,xmB)
     &    +suda(x3,bB,xmB)+suda((1.-x3),bB,xmB)
      ssu=dexp(-sud)
      if (ssu .gt. 1.)  ssu=1.
      su1=ssu*dexp((5./3.+2.)*ss(bB,t1a)+2.*ss(b,t1a))
      su2=ssu*dexp((5./3.+2.)*ss(bB,t1b)+2.*ss(b,t1b))

      phiBx=phib(xB,bB)
      call DAh(x3,vah) !for h
      phAx =vah(1)
      phPx =vah(2)
      phTx =vah(3)
      call DAR(z,vpR)  !for R
      phi0z =vpR(1)
      phisz =vpR(2)
      phitz =vpR(3)

      !------------FLL---------
      ctmp=4./3*16*Sqrt(2./3)*phi0z*phiBx*Pi*xmB**4*bB*b
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=(-1+et)**2*phAx*(-1+xB+z)+r0*(-((-1+et)*(phPx-phTx)*x3)
     &      +et*phTx*(-2+xB+z)+et*phPx*(xB+z))
      factb=(-1+et)*((-1+et)*phAx+(phPx+phTx)*r0)*x3+((-1+et)*phAx+et
     &      *(-phPx+phTx)*r0)*(xB-z)

      FLLar=coea*facta*har
      FLLai=coea*facta*hai
      FLLbr=coeb*factb*hbr
      FLLbi=coeb*factb*hbi
      !------------FLL------end
      !------------Fsp---------
      ctmp=4./3*16*Sqrt(2./3)*phi0z*phiBx*Pi*xmB**4*bB*b
      coea=ctmp*alphas(t1a)*su1
      coeb=ctmp*alphas(t1b)*su2

      facta=-((-1+et)*phAx*(-1+et*(-1+x3)-x3+xB+z))+r0*(-((-1+et)
     &      *(phPx+phTx)*x3)-et*phTx*(-2+xB+z)+et*phPx*(xB+z))
      factb=(-1+et)*(phPx-phTx)*r0*x3-((-1+et)**2*phAx+et*(phPx+phTx)
     &      *r0)*xB+((-1+et)**2*phAx+et*(phPx+phTx)*r0)*z

      Fspar=coea*facta*har
      Fspai=coea*facta*hai
      Fspbr=coeb*factb*hbr
      Fspbi=coeb*factb*hbi
      !------------Fsp------end
      !------------wc---------
      if(t1a.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1a**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1a)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcax=wc(2)
      wcayA=wc(4)+wc(10)   !for uu_bar
      wcazA=wc(6)+wc(8)
      wcay=wcayA -( wc(4)-0.5*wc(10) )  !for dd_bar
      wcaz=wcazA -( wc(6)-0.5*wc(8) )
      !------------
      if(t1b.gt.4.8) then  !xmw=80.4 vmbq=4.8
      aa=dlog(80.4**2./xlam**2.)/dlog(t1b**2./xlam**2.)
      call wilsw(aa,wc)
      else
      aa=alphas(t1b)/alphas(vmbq)
      call wilsb(aa,wc)
      end if
      wcbx=wc(2)
      wcbyB=wc(4)+wc(10)   !for uu_bar
      wcbzB=wc(6)+wc(8)
      wcby=wcbyB -( wc(4)-0.5*wc(10) ) !for dd_bar
      wcbz=wcbzB -( wc(6)-0.5*wc(8) )
      !------------wc------end
      vur=ckA*(ckl**4.)*ckr
      vui=ckA*(ckl**4.)*cke
      vtr=-ckA*(ckl**2.)
      vti=0.

      f1ar=wcax*(vur*FLLar-vui*FLLai)-wcay*(vtr*FLLar-vti*FLLai)
     &    -wcaz*(vtr*FSPar-vti*FSPai)
      f1ai=wcax*(vur*FLLai+vui*FLLar)-wcay*(vtr*FLLai+vti*FLLar)
     &    -wcaz*(vtr*FSPai+vti*FSPar)
      f1br=wcbx*(vur*FLLbr-vui*FLLbi)-wcby*(vtr*FLLbr-vti*FLLbi)
     &    -wcbz*(vtr*FSPbr-vti*FSPbi)
      f1bi=wcbx*(vur*FLLbi+vui*FLLbr)-wcby*(vtr*FLLbi+vti*FLLbr)
     &    -wcbz*(vtr*FSPbi+vti*FSPbr)

	v1op(1)=(f1ar+f1br)/sqrt(2.) !*xlamb 
	v1op(2)=(f1ai+f1bi)/sqrt(2.) !*xlamb  
      return
      end
      !------------------------
c############################## Decay Amplitudes ## END #


!********************************************************! 
!*            ---- W.F. Wang 2018-07-06----             *!    
!*        --- Hard, VEGAS & BESSEL FUNCTIONS ---        *!
!********************************************************!   
     	!xxxxxxxxxxxxxxxxxx H FUN. xxxxxxxxxx Begin xxxx  
      subroutine hfun(al,be,b_1,b_2,vfun)		
	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
	double precision vfun(2)
	alp=xmB*dsqrt(dabs(al))
	bet=xmB*dsqrt(dabs(be))
	if(be .ge. 0.00000001)then  
	  facxr=fbsk0(1,bet*b_2,ier)
	  facxi=0.0
	elseif(be .ge. -0.00000001)then
	   vfun(1)=0.0
	   vfun(2)=0.0
	   goto 99
	else
	  facxr=-pi/2.*fbsn0(bet*b_2)
	  facxi=pi/2.*fbsj0(bet*b_2,ier)
	endif   
	if(al .ge. 0.00000001)then
	 goto 10
	elseif(al .ge. -0.00000001)then
	   vfun(1)=0.0
	   vfun(2)=0.0
	   goto 99
	else
	 goto 20
	endif
10	facyi=0.0
	if(b_1 .ge. b_2)then
	facyr=fbsi0(1,alp*b_2,ier)*fbsk0(1,alp*b_1,ier)
	else
	facyr=fbsi0(1,alp*b_1,ier)*fbsk0(1,alp*b_2,ier)
	endif
	goto 88
20	facyi=pi/2.*fbsj0(alp*b_2,ier)*fbsj0(alp*b_1,ier)
	if(b_1 .ge. b_2)then
	facyr=-pi/2.*fbsj0(alp*b_2,ier)*fbsn0(alp*b_1)
	else
	facyr=-pi/2.*fbsj0(alp*b_1,ier)*fbsn0(alp*b_2)
	endif
88	vfun(1)=facxr*facyr-facxi*facyi      !  get real part
	vfun(2)=facxi*facyr+facxr*facyi      !  get imag part
99	continue
 	return
      end
	!xxxxxxxxxxxxxxx H Fun. xxxxxxxxxxxxxxxxxx END xxxx
	
	!--------  N_B for B meson  --------
	double precision function funNB(x,wgt)    
    	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
      dimension x(1)
	y=x(1)
   	funNB=y*y*(1.-y)**2.*dexp(-0.5*(y*xmB/wb)**2.)  
	return
      end
	!------------------------
	double precision function phib(x,b)  ! B-meson
      implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
      fa=xnb*(x*(1.-x))**2.0
	phib=fa*dexp(-0.5*(x*xmB/wb)**2.-0.5*(wb*b)**2.)
	return
      end
	!+++++++++++++++++ alphas +++++++++++++++++!
	double precision function alphas(u)
	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
	xbeta0=(33.-2.*4.)/3.
	xbeta1=(306.-38.*4.)/3.
	 ct=dlog(u**2./(xlamb**2.))    
	alphas=4.*pi/xbeta0/ct !*(1.-xbeta1/xbeta0/xbeta0*dlog(ct)/ct)
	return	 
	End
	!---------------
      double precision function alphas0(u)
	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
	xbeta0=(33.-2.*5.)/3.
	xbeta1=(306.-38.*5.)/3.
	 ct=dlog(u**2./(xlam**2.))
	alphas0=4.*pi/xbeta0/ct !*(1.-xbeta1/xbeta0/xbeta0*dlog(ct)/ct)
	return
	End
	!xxxxxxxxxxxxxx    Sudakov fators	xxxxxxxxxxxxxx
      double precision function suda(x,b,qm)
	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
	!---------------- 
	euler=0.577215665
      a1=1.333333333333333
	beta1=2.083333333333333            !(33.0-2.0*xnf)/12.0
	xnf=4.0
	a2=67./9.-pi**2./3.-10.*xnf/27.+8./3.*beta1*dlog(dexp(euler)/2.)
	beta2=0.0
	!----------------
      qhat=dlog(x*qm/sqrt(2.)/xlamb)
      bhat=dlog(1./b/xlamb)
      if (qhat .gt. bhat) then
        suda=a1*qhat*dlog(qhat/bhat)/2./beta1
     &    -a1*(qhat-bhat)/2./beta1
     &    +a2*(qhat/bhat-1.)/4./beta1**2.
     &    -(a2/4./beta1**2.-a1*dlog(dexp(2.*euler-1.)/2.)/4./beta1)
     &             *dlog(qhat/bhat)
      else
        suda=0.
      endif
      return
      end
	!------------------------
      double precision function ss(b,t)
	implicit double precision(a-h,o-z)
	common/con1/xlamb,xlam,xmw,pi,xmB,vmbq,fB,wB,xNB
      if (t .gt. 1./b) then
        s2=0.48*dlog(-dlog(t/xlamb)/dlog(b*xlamb))	!1./beta1=0.48
        ss=s2/2.0 
      else
        ss=0.
      endif
      return
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c************** Wilson Coefs. at LO  2015.01.22 **************
      subroutine wilsw(aa,c)
      implicit double precision(a-h,o-z)
      double precision  c(10)
c     aa=alphas(t)/alphas(m_w)
      if(aa.lt.1.)aa=1.

        amw=80.41
        Pi=3.14159
        xt=173.8**2./amw**2.
        beta5=(33.-2.*5.)/12.
        asw=Pi/beta5/LOG(amw**2/0.193**2.)
        b0=(xt/(1.-xt)+xt/(xt-1.)**2.*log(xt))/4.
        c0=xt*((xt-6.)/(xt-1.)+(3.*xt+2.)/(xt-1.)**2.*log(xt))/8.
        d0=-4./9.*log(xt)+(-19.*xt**3.+25.*xt**2.)/(xt-1.)**3./36.
     &    +xt**2.*(5.*xt**2.-2.*xt-6.)/(xt-1.)**4.*log(xt)/18.

        e0=-2./3.*(log(xt))+xt*(xt**2.+11.*xt-18.)/(xt-1.)**3./12.
     &    +xt**2.*(4.*xt**2.-16.*xt+15.)/(xt-1.)**4.*log(xt)/6.
        c7w=(4.*c0+d0)/132./6./pi
        c9w=(4.*c0+d0+(10.*b0-4.*c0)/0.23124)/132./6./pi
        c4=asw/8./pi*e0
        c6=c4
        c3=-asw/24./pi*e0+(2.*b0+c0)/0.23124/132./6./pi
        c5=-asw/24./pi*e0

        c(1)= 0.5*aa**(-6./23.) - 0.5*aa**(12./23.)
        c(2)= 0.5*aa**(-6./23.) + 0.5*aa**(12./23.)

        c(3)= 0.0510*aa**(-0.4086) - 0.0714*aa**(-6./23.) +
     -   0.0054*aa**(-0.1456) -0.1403*aa**0.4230 -
     -   0.0113*aa**0.8994 + 1./6.*aa**(12./23.) +
     h   c9w*(-0.0325*aa**(-0.4086) + 0.0357*aa**(-6./23.)-
     -   0.0016*aa**(-0.1456) +0.2342*aa**0.4230 -
     -   0.25*aa**(12./23.) +0.0141*aa**0.8994) +
     h   c7w*(-0.0063*aa**(-0.4086)  + 0.01629*aa**(-0.1456)-
     -   0.0185*aa**0.4230 +0.0085*aa**0.8994 )+
     -  c4*(0.3287*aa**(-0.4086)  +   0.0424*aa**(-0.1456) -
     -   0.3263*aa**0.4230  -   0.0448*aa**0.8994) +
     -  c3*(0.2868*aa**(-0.4086) +  0.0491*aa**(-0.1456) +
     -   0.6579*aa**0.4230  +   0.0061*aa**0.8994)+
     -  c5*(-0.0629*aa**(-0.4086) + 0.1629*aa**(-0.1456) -
     -   0.1846*aa**0.4230 +  0.0846*aa**0.8994) +
     -  c6*(0.0447*aa**(-0.4086) -  0.0063*aa**(-0.1456) -
     -   0.2610*aa**0.4230+  0.2226*aa**0.8994)

        c(4)=  0.0984*aa**(-0.4086) - 0.0714*aa**(-6./23.) +
     -   0.0026*aa**(-0.1456) + 0.1214*aa**0.4230 -
     -   1./6.*aa**(12./23.)+0.0156*aa**0.8994 +
     h   c9w*(-0.0627*aa**(-0.4086) + 0.0357*aa**(-6./23.) -
     -   0.0008*aa**(-0.1456) - 0.2027*aa**0.4230 +
     -   0.25*aa**(12./23.) -0.0196*aa**0.8994)  +
     h   c7w*(-0.0122*aa**(-0.4086) +0.0079*aa**(-0.1456)  +
     -   0.0160*aa**0.4230  - 0.0117*aa**0.8994 )+
     -  c6*(0.0864*aa**(-0.4086)  - 0.0031*aa**(-0.1456) +
     -   0.2259*aa**0.4230  - 0.3092*aa**0.8994) +
     -  c5*(-0.1215*aa**(-0.4086)  +  0.0793*aa**(-0.1456) +
     -   0.1597*aa**0.4230  -  0.1175*aa**0.8994)+
     -  c3*(0.5539*aa**(-0.4086) + 0.0239*aa**(-0.1456) -
     -   0.5693*aa**0.4230  - 0.0085*aa**0.8994) +
     -  c4*(0.6348*aa**(-0.4086)  +  0.0206*aa**(-0.1456) +
     -   0.2823*aa**0.4230+  0.0623*aa**0.8994)

        c(5)=  -0.0397*aa**(-0.4086)  + 0.0304*aa**(-0.1456) +
     -   0.0117*aa**0.4230 - 0.0025*aa**0.8994 +
     h  c9w*(0.0253*aa**(-0.4086)  - 0.0089*aa**(-0.1456) -
     -   0.0195*aa**0.4230  +0.0031*aa**0.8994)  +
     h   c7w*(0.0049*aa**(-0.4086)  + 0.0917*aa**(-0.1456) -
     -   0.1*aa**(-3./23.) +0.0015*aa**0.4230 +0.0019*aa**0.8994)+
     g  c4*(-0.2559*aa**(-0.4086)+ 0.2385*aa**(-0.1456) +
     -   0.0271*aa**0.4230  -  0.0098*aa**0.8994)+
     -  c3*(-0.2233*aa**(-0.4086) +  0.2767*aa**(-0.1456) -
     -   0.0547*aa**0.4230  + 0.0013*aa**0.8994)+
     -  c5*(0.0490*aa**(-0.4086)+ 0.9171*aa**(-0.1456) +
     -   0.0154*aa**0.4230 +  0.0185*aa**0.8994) +
     -  c6*(-0.0348*aa**(-0.4086) - 0.0357*aa**(-0.1456) +
     -   0.0217*aa**0.4230 + 0.0488*aa**0.8994)

        c(6)=0.0335*aa**(-0.4086)  - 0.0112*aa**(-0.1456) +
     -   0.0239*aa**0.4230 - 0.0462*aa**0.8994 +
     h  c9w*(-0.0213*aa**(-0.4086) + 0.0033*aa**(-0.1456) -
     -   0.0399*aa**0.4230  +0.0579*aa**0.8994)  +
     h  c7w*(-0.0041*aa**(-0.4086)  - 0.0337*aa**(-0.1456) +
     -   aa**(-3./23.)/30. +0.0031*aa**0.4230 +
     -   0.0347*aa**0.8994 -aa**(24./23.)/30.)+
     g  c4*(0.2160*aa**(-0.4086)  -  0.0877*aa**(-0.1456) +
     -   0.0555*aa**0.4230 -  0.1839*aa**0.8994) +
     -  c3*(0.1885*aa**(-0.4086) - 0.1017*aa**(-0.1456) -
     -   0.1120*aa**0.4230  + 0.0251*aa**0.8994)+
     -  c5*(-0.0414*aa**(-0.4086) - 0.3370*aa**(-0.1456) +
     -   0.0314*aa**0.4230 +  0.3469*aa**0.8994) +
     -  c6*(0.0294*aa**(-0.4086)  + 0.0131*aa**(-0.1456) +
     -   0.0444*aa**0.4230  +   0.9131*aa**0.8994)

        c(7)=c7w*aa**(-3./23.)
        c(8)=c7w*(-aa**(-3./23.) +aa**(24./23.))/3.
        c(9)= c9w*( 0.5*aa**(-6./23.)  + 0.5*aa**(12./23.) )
        c(10)= c9w*(0.5*aa**(-6./23.) - 0.5*aa**(12./23.) )

      return
      end
c****************************************************************
      subroutine wilsb(aa,cc)
      implicit double precision(a-h,o-z)
	double precision  cc(10)
c     aa=alphas(t)/alphas(m_b)
      if(aa.lt.1.)aa=1.

        C1= -0.27033535
        C2=  1.11879098
        C3=  0.0126053
        C4= -0.02694983
        C5=  0.00847361
        C6= -0.03259826
        C7=  0.00108502
        C8=  0.00039601
        C9= -0.00894514
        C10= 0.00216143

        cc(1)= c2*( aa**(-6./25.) -aa**(12./25.) ) /2.
     -   + c1*(  aa**(-6./25.) +  aa**(12./25.)) /2.
        cc(2)= c2*( aa**(-6./25.) +  aa**(12./25.))/2. +
     -    c1*( aa**(-6./25.)  - aa**(12./25.)  )/2.

        cc(3)= c4*(0.3606*aa**(-0.3469) +  0.03166*aa**(-0.1317) -
     -  0.3626*aa**0.4201  - 0.0297*aa**0.8451) +
     - c10*(0.0149*aa**(-0.3469)  - 0.0020*aa**(-0.1317) -
     -  0.4981*aa**0.4201 + 0.5*aa**(12./25.) - 0.0148*aa**0.8451) +
     -  c2*(0.0651*aa**(-0.3469) -  0.0833*aa**(-6./25.) +
     -  0.0046*aa**(-0.1317) - 0.2265*aa**0.4201 + 0.25*aa**(12./25.)-
     u  0.0099*aa**0.8451) +
     9  c3*(0.3308*aa**(-0.3469)  + 0.0356*aa**(-0.1317) +
     g  0.6337*aa**0.4201  -  0.0001*aa**0.8451) +
     f  c1*(0.0502*aa**(-0.3469) - 0.0833*aa**(-6./25.) +
     y   0.0066*aa**(-0.1317)+0.2717*aa**0.4201-0.25*aa**(12./25.)+
     y   0.0049*aa**0.8451) +
     y  c9*(-0.0149*aa**(-0.3469) + 0.0020*aa**(-0.1317) +
     y    0.4981*aa**0.4201 - 0.5*aa**(12./25.) + 0.0148*aa**0.8451)+
     y  c5*(-0.0598*aa**(-0.3469) + 0.1371*aa**(-0.1317) -
     y   0.1473*aa**0.4201 +   0.0700*aa**0.8451) +
     y  c6*(0.0377*aa**(-0.3469) -   0.0045*aa**(-0.1317) -
     y   0.2210*aa**0.4201  + 0.18775*aa**0.8451) +
     y  c7*(-0.0150*aa**(-0.3469) + 0.0343*aa**(-0.1317)  -
     y   0.0368*aa**0.4201  +  0.0175*aa**0.8451) +
     y  c8*(0.009*aa**(-0.3469) -  0.0011*aa**(-0.1317) -
     y   0.0553*aa**0.4201  + 0.0469*aa**0.8451 )

        cc(4)=c6*(0.0640*aa**(-0.3469) - 0.0021*aa**(-0.1317) +
     y   0.2018*aa**0.4201  - 0.2637*aa**0.8451
     y  ) + c5*(-0.10156*aa**(-0.3469)  +  0.06538*aa**(-0.1317) +
     y   0.1345*aa**0.4201  -  0.09836*aa**0.8451) +
     y  c9*(-0.02528*aa**(-0.3469)  + 0.0009*aa**(-0.1317) -
     y   0.4549*aa**0.4201 + 0.5*aa**(12./25.) - 0.0207*aa**0.8451) +
     y  c1*(0.08515*aa**(-0.3469) -  0.0833*aa**(-6./25.) +
     y  0.0031*aa**(-0.1317)-0.24809*aa**0.4201+0.25*aa**(12./25.)-
     y  0.00688*aa**0.8451) +
     y  c3*(0.5615*aa**(-0.3469)  +  0.01699*aa**(-0.1317) -
     y    0.5787*aa**0.4201 + 0.0002*aa**0.8451) +
     y  c2*(0.1104*aa**(-0.3469)-0.0833*aa**(-6./25.)
     g      +0.0022*aa**(-0.1317)
     y  +0.2068*aa**0.4201 - 0.25*aa**(12./25.) + 0.0139*aa**0.8451) +
     y   c10*(0.0253*aa**(-0.3469) - 0.0009*aa**(-0.1317) +
     y   0.4549*aa**0.4201 - 0.5*aa**(12./25.) + 0.0207*aa**0.8451) +
     y  c4*(0.6121*aa**(-0.3469) +  0.0151*aa**(-0.1317) +
     y    0.3311*aa**0.4201  +  0.0417*aa**0.8451) +
     y c8*(0.0160*aa**(-0.3469)  - 0.0005*aa**(-0.1317) +
     y    0.0505*aa**0.4201  -  0.0659*aa**0.8451) +
     y  c7*(-0.0254*aa**(-0.3469) +
     y   0.0163*aa**(-0.1317) + 0.0336*aa**0.4201 - 0.0246*aa**0.8451)

        cc(5)=c4*(-0.2291*aa**(-0.3469) +  0.2167*aa**(-0.1317) +
     y  0.0192*aa**0.4201  - 0.0067*aa**0.8451) +
     y  c10*(-0.0095*aa**(-0.3469) -  0.0136*aa**(-0.1317) +
     y  0.0264*aa**0.4201  -  0.0034*aa**0.8451) +
     y c2*(-0.0413*aa**(-0.3469) + 0.0316*aa**(-0.1317) +
     y    0.0120*aa**0.4201  - 0.0022*aa**0.8451
     y  ) + c3*(-0.2102*aa**(-0.3469) +  0.2438*aa**(-0.1317) -
     y  0.0336*aa**0.4201 ) +
     y c1*(-0.0319*aa**(-0.3469)  + 0.0452*aa**(-0.1317) -
     y  0.0144*aa**0.4201  + 0.0011*aa**0.8451) +
     y   c9*(0.0095*aa**(-0.3469)  + 0.0136*aa**(-0.1317) -
     y  0.0264*aa**0.4201  +  0.0034*aa**0.8451) +
     y  c5*(0.0380*aa**(-0.3469) +  0.9382*aa**(-0.1317) +
     y   0.0078*aa**0.4201 + 0.0159*aa**0.8451) +
     y  c6*(-0.0240*aa**(-0.3469) - 0.0305*aa**(-0.1317) +
     y    0.0117*aa**0.4201  + 0.0427*aa**0.8451) +
     y c8*(-0.0060*aa**(-0.3469)  - 0.0076*aa**(-0.1317) +
     y  0.0029*aa**0.4201  + 0.0107*aa**0.8451  )+
     y c7*(0.0095*aa**(-0.3469)  +0.2346*aa**(-0.1317) -
     y  0.25*aa**(-3./25.) + 0.0020*aa**0.4201 +
     y  0.0040*aa**0.8451 )

        cc(6)= c4*(0.1825*aa**(-0.3469) -0.0784*aa**(-0.1317) +
     -  0.0449*aa**0.4201  - 0.14894*aa**0.8451)
     -   + c10*(0.0075*aa**(-0.3469) + 0.0049*aa**(-0.1317) +
     -   0.0617*aa**0.4201  - 0.07412*aa**0.8451
     -  ) + c2*(0.0329*aa**(-0.3469) -0.0114*aa**(-0.1317) +
     -  0.0280*aa**0.4201  -  0.0495*aa**0.8451) +
     - c3*(0.1674*aa**(-0.3469) -0.0882*aa**(-0.1317) -
     -  0.0784*aa**0.4201 -0.0007*aa**0.8451) +
     -  c1*(0.0254*aa**(-0.3469) - 0.0163*aa**(-0.1317) -
     -  0.0336*aa**0.4201  +  0.0246*aa**0.8451) +
     - c9*(-0.0075*aa**(-0.3469)  - 0.0049*aa**(-0.1317) -
     -  0.0617*aa**0.4201  + 0.07412*aa**0.8451
     - ) + c5*(-0.0303*aa**(-0.3469) - 0.3395*aa**(-0.1317) +
     -  0.0182*aa**0.4201 + 0.35157*aa**0.8451
     -  ) + c6*(0.0191*aa**(-0.3469) +  0.0110*aa**(-0.1317) +
     -  0.0274*aa**0.4201 + 0.94253*aa**0.8451
     -  ) + c8*(0.0048*aa**(-0.3469) + 0.0028*aa**(-0.1317) +
     -  0.0068*aa**0.4201 +0.2356*aa**0.8451 - 0.25*aa**(24./25.)) +
     -     c7*(-0.0076*aa**(-0.3469) -  0.0849*aa**(-0.1317) +
     -  0.0833*aa**(-3./25.) + 0.0046*aa**0.4201 +
     -  0.0879*aa**0.8451 - 0.0833*aa**(24./25.))

        cc(7)= c7*aa**(-3./25.)
        cc(8)=c7*(-aa**(-3./25.)+aa**(24./25.))/3. +c8*aa**(24./25.)

        cc(9)= c10*(0.5*aa**(-6./25.)-0.5*aa**(12./25.))+
     -     c9*( 0.5*aa**(-6./25.)+0.5*aa**(12./25.) )

        cc(10) =c9*( 0.5*aa**(-6./25.)-0.5*aa**(12./25.)) +
     -     c10*(0.5*aa**(-6./25.) + 0.5*aa**(12./25.))

      return
      end
cccccccccccccccccccccccccc - Vegas - cccccccccccccccccccccc    	  
      BLOCK DATA
C
C   MAKES DEFAULT PARAMETER ASSIGNMENTS FOR VEGAS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(20),XU(20),ACC
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,20)
      COMMON/BVEG3/ALPH,NDMX,MDS
      COMMON/RNSD/DSEED
      DATA NCALL/50000/,ITMX/10/,NPRN/-1/,ACC/-1./,
     1     XL/20*0.0000001/,
     2     XU/20*0.9999999/,     ! xu(1,2,3,4,5,6) for x_{1,2,3,1',2',3'}, respectively
     3     ALPH/1.5/,NDMX/50/,MDS/1/,NDEV/6/,
     4     NDO/1/,XI/1000*1./,IT/0/,SI,SWGT,SCHI/3*0./,
     5     DSEED/12345678./
      END
C
      SUBROUTINE VEGAS(NDIM,FXN,AVGI,SD,CHI2A)
C
C   SUBROUTINE PERFORMS NDIM-DIMENSIONAL MONTE CARLO INTEG'N
C      - BY G.P. LEPAGE    SEPT 1976/(REV)AUG 1979
C      - ALGORITHM DESCRIBED IN J COMP PHYS 27,192(1978)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(20),XU(20),ACC
      COMMON/BVEG2/IT,NDO,SI,SWGT,SCHI,XI(50,20)
      COMMON/BVEG3/ALPH,NDMX,MDS
      COMMON/BVEG4/CALLS,TI,TSI
      DIMENSION D(50,20),DI(50,20),XIN(50),R(50),DX(20),IA(20),
     1          KG(20),DT(20),X(20)
      DIMENSION RAND(20)
      DATA ONE/1./
      SQRT(A)=DSQRT(A)
      ALOG(A)=DLOG(A)
      ABS(A)=DABS(A)
C
      NDO=1
      DO 1 J=1,NDIM
1     XI(1,J)=ONE
C
      ENTRY VEGAS1(NDIM,FXN,AVGI,SD,CHI2A)
C         - INITIALIZES CUMULATIVE VARIABLES, BUT NOT GRID
      IT=0
      SI=0.
      SWGT=SI
      SCHI=SI
C
      ENTRY VEGAS2(NDIM,FXN,AVGI,SD,CHI2A)
C         - NO INITIALIZATION
      ND=NDMX
      NG=1
      IF(MDS.EQ.0) GO TO 2
      NG=(NCALL/2.)**(1./NDIM)
      MDS=1
      IF((2*NG-NDMX).LT.0) GO TO 2
      MDS=-1
      NPG=NG/NDMX+1
      ND=NG/NPG
      NG=NPG*ND
2     K=NG**NDIM
      NPG=NCALL/K
      IF(NPG.LT.2) NPG=2
      CALLS=NPG*K
      DXG=ONE/NG
      DV2G=(CALLS*DXG**NDIM)**2/NPG/NPG/(NPG-ONE)
      XND=ND
      NDM=ND-1
      DXG=DXG*XND
      XJAC=ONE/CALLS
      DO 3 J=1,NDIM
      DX(J)=XU(J)-XL(J)
3     XJAC=XJAC*DX(J)
C
C   REBIN, PRESERVING BIN DENSITY
      IF(ND.EQ.NDO) GO TO 8
      RC=NDO/XND
      DO 7 J=1, NDIM
      K=0
      XN=0.
      DR=XN
      I=K
4     K=K+1
      DR=DR+ONE
      XO=XN
      XN=XI(K,J)
5     IF(RC.GT.DR) GO TO 4
      I=I+1
      DR=DR-RC
      XIN(I)=XN-(XN-XO)*DR
      IF(I.LT.NDM) GO TO 5
      DO 6 I=1,NDM
6     XI(I,J)=XIN(I)
7     XI(ND,J)=ONE
      NDO=ND
C
8     IF(NPRN.GE.0) WRITE(NDEV,200) NDIM,CALLS,IT,ITMX,ACC,NPRN,
     1                    ALPH,MDS,ND,(XL(J),XU(J),J=1,NDIM)
C
      ENTRY VEGAS3(NDIM,FXN,AVGI,SD,CHI2A)
C         - MAIN INTEGRATION LOOP
9     IT=IT+1
      TI=0.
      TSI=TI
      DO 10 J=1,NDIM
      KG(J)=1
      DO 10 I=1,ND
      D(I,J)=TI
10    DI(I,J)=TI
C
11    FB=0.
      F2B=FB
      K=0
12    K=K+1
      CALL RANDA(NDIM,RAND)
      WGT=XJAC
      DO 15 J=1,NDIM
      XN=(KG(J)-RAND(J))*DXG+ONE
      IA(J)=XN
      IF(IA(J).GT.1) GO TO 13
      XO=XI(IA(J),J)
      RC=(XN-IA(J))*XO
      GO TO 14
13    XO=XI(IA(J),J)-XI(IA(J)-1,J)
      RC=XI(IA(J)-1,J)+(XN-IA(J))*XO
14    X(J)=XL(J)+RC*DX(J)
15    WGT=WGT*XO*XND
C
      F=WGT
      F=F*FXN(X,WGT)
      F2=F*F
      FB=FB+F
      F2B=F2B+F2
      DO 16 J=1,NDIM
      DI(IA(J),J)=DI(IA(J),J)+F
16    IF(MDS.GE.0) D(IA(J),J)=D(IA(J),J)+F2
      IF(K.LT.NPG) GO TO 12
C
      F2B=SQRT(F2B*NPG)
      F2B=(F2B-FB)*(F2B+FB)
      TI=TI+FB
      TSI=TSI+F2B
      IF(MDS.GE.0) GO TO 18
      DO 17 J=1,NDIM
17    D(IA(J),J)=D(IA(J),J)+F2B
18    K=NDIM
19    KG(K)=MOD(KG(K),NG)+1
      IF(KG(K).NE.1) GO TO 11
      K=K-1
      IF(K.GT.0) GO TO 19
C
C   COMPUTE FINAL RESULTS FOR THIS ITERATION
      TSI=TSI*DV2G
      TI2=TI*TI
      WGT=ONE/TSI
      SI=SI+TI*WGT
      SWGT=SWGT+WGT
      SCHI=SCHI+TI2*WGT
      AVGI=SI/SWGT
      CHI2A=(SCHI-SI*AVGI)/(IT-.9999)
      SD=SQRT(ONE/SWGT)
C
      IF(NPRN.LT.0) GO TO 21
      TSI=SQRT(TSI)
      WRITE(NDEV,201) IT,TI,TSI,AVGI,SD,CHI2A
      IF(NPRN.EQ.0) GO TO 21
      DO 20 J=1,NDIM
20    WRITE(NDEV,202) J,(XI(I,J),DI(I,J),I=1+NPRN/2,ND,NPRN)
C
C   REFINE GRID
21    DO 23 J=1,NDIM
      XO=D(1,J)
      XN=D(2,J)
      D(1,J)=(XO+XN)/2.
      DT(J)=D(1,J)
      DO 22 I=2,NDM
      D(I,J)=XO+XN
      XO=XN
      XN=D(I+1,J)
      D(I,J)=(D(I,J)+XN)/3.
22    DT(J)=DT(J)+D(I,J)
      D(ND,J)=(XN+XO)/2.
23    DT(J)=DT(J)+D(ND,J)
C
      DO 28 J=1,NDIM
      RC=0.
      DO 24 I=1,ND
      R(I)=0.
      IF(D(I,J).LE.0.) GO TO 24
      XO=DT(J)/D(I,J)
      R(I)=((XO-ONE)/XO/ALOG(XO))**ALPH
24    RC=RC+R(I)
      RC=RC/XND
      K=0
      XN=0.
      DR=XN
      I=K
25    K=K+1
      DR=DR+R(K)
      XO=XN
      XN=XI(K,J)
26    IF(RC.GT.DR) GO TO 25
      I=I+1
      DR=DR-RC
      XIN(I)=XN-(XN-XO)*DR/R(K)
      IF(I.LT.NDM) GO TO 26
      DO 27 I=1,NDM
27    XI(I,J)=XIN(I)
28    XI(ND,J)=ONE
C
      IF(IT.LT.ITMX.AND.ACC*ABS(AVGI).LT.SD) GO TO 9
200   FORMAT(/35H INPUT PARAMETERS FOR VEGAS:  NDIM=I3,8H  NCALL=F8.0
     1  /28X,5H  IT=I5,7H  ITMX=I5/28X,6H  ACC=G9.3
     2  /28X,7H  NPRN=I3,7H  ALPH=F5.2/28X,6H  MDS=I3,6H   ND=I4
     3  /28X,10H  (XL,XU)=(T40,2H( G12.6,3H , G12.6,2H )))
201   FORMAT(///21H INTEGRATION BY VEGAS//14H ITERATION NO.,I3,
     1  14H:   INTEGRAL =G14.8/24X,10HSTD DEV  =G10.4/
     2  34H ACCUMULATED RESULTS:   INTEGRAL =G14.8/
     3  24X,10HSTD DEV  =G10.4/24X,17HCHI**2 PER IT'N =G10.4)
202   FORMAT(/15H DATA FOR AXIS I2/25H    X       DELTA I       ,
     1  24H   X       DELTA I      ,18H   X       DELTA I,
     2  /(1H F7.6,1X,G11.4,5X,F7.6,1X,G11.4,5X,F7.6,1X,G11.4))
      RETURN
      END
C
      SUBROUTINE RANDA(NR,R)
C                                  SPECIFICATIONS FOR ARGUMENTS
      DOUBLE PRECISION DSEED,R
      DIMENSION R(NR)
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER            I
      DOUBLE PRECISION   D2P31M,D2P31
C                                  D2P31M=(2**31) - 1
C                                  D2P31 =(2**31)(OR AN ADJUSTED VALUE)
      COMMON/RNSD/DSEED
      DATA               D2P31M/2147483647.D0/
      DATA               D2P31/2147483711.D0/
C                                  FIRST EXECUTABLE STATEMENT
      DO 5 I=1,NR
         DSEED = DMOD(16807.D0*DSEED,D2P31M)
    5 R(I) = DSEED / D2P31
      RETURN
      END	   
cxxxxxxxxxxxxxxxxxxxxxxx Bessel functions xxxxxxxxxxxxxxxxxxxxxxx
      double precision function fbsn0(arg)
      implicit double precision(a-h,o-z)
      call fbsyn(arg,0.0d0,1,yn,ier)
      fbsn0=yn
      return
      end
	!--------------------------------------------
      double precision function fbsi0 (iopt,arg,ier)
      implicit double precision(a-h,o-z)
      dimension p(15),pp(8),q(5),qq(7)
      data xmax /91.200d0/
      data dexp40/ .235385266837019985d+18/
      data xsmall/8.67361738d-19/
      data xinf /1.7d+38/
      data const / .666666666666666667d-01/
      data p /-.524878303872486911d-17,-.159822176963152107d-13,
     *        -.268434353985728351d-10,-.305172129516732270d-07,
     *        -.251726344203303272d-04,-.154539719077661595d-01,
     *        -.709353218849083214d+01,-.241251875635161846d+04,
     *        -.595456062371357883d+06,-.103130633806846395d+09,
     *        -.119127423465868193d+11,-.849250748990540193d+12,
     *        -.329400775241734956d+14,-.550503528993272799d+15,
     *        -.223355758483236763d+16/
      data q /-.372775563866849905d+04, .651584945870435991d+07,
     *        -.656265446265324359d+10, .376041775067259003d+13,
     *        -.970879130649400689d+15/
      data pp/-.398437500000000000d+00, .292053801142487550d+01,
     *        -.247084602559363172d+01, .479148550995084342d+00,
     *        -.373847878179342018d-02,-.268014935222659065d-02,
     *         .991686607024382073d-04,-.218771393882892240d-05/
      data qq/-.314466876245542596d+02, .855395457249381542d+02,
     *        -.602279794663117256d+02, .139825866959438771d+02,
     *        -.111517488948088167d+01, .325476719383372820d-01,
     *        -.551943584858776969d-03/
      ier = 0
      xcon = 85.0d0
      if (iopt.ne.1.and.iopt.ne.2) go to 30
      x = dabs(arg)
      if (x.lt.xsmall) go to 20
      if ((iopt.eq.1).and.(x.gt.xmax)) go to 25
      if (x.ge.15.0d0) go to 10
      xx = x*x
      sump = p(1)
      do 5 i=2,15
         sump = sump*xx+p(i)
    5 continue
      xx = xx-225.d0
      sumq = ((((xx+q(1))*xx+q(2))*xx+q(3))*xx+q(4))*xx+q(5)
      fbsi0 = sump/sumq
      if (iopt.eq.2) fbsi0 = fbsi0*dexp(-x)
      go to 9005
   10 xx = 1.0d0/x-const
      sump = ((((((pp(1)*xx+pp(2))*xx+pp(3))*xx+pp(4))*xx+pp(5))*xx+
     *pp(6))*xx+pp(7))*xx+pp(8)
      sumq = ((((((xx+qq(1))*xx+qq(2))*xx+qq(3))*xx+qq(4))*xx+qq(5))*xx+
     *qq(6))*xx+qq(7)
      fbsi0 = (sump/sumq-pp(1))/dsqrt(x)
      if (iopt.eq.2) go to 9005
      if (x.gt.xcon) go to 15
      fbsi0 = fbsi0*dexp(x)
      go to 9005
   15 a = fbsi0*dexp(x-40.0d0)
      fbsi0 = a*dexp40
      go to 9005
   20 fbsi0 = 1.0d0
      go to 9005
   25 fbsi0 = xinf
      ier = 130
      go to 9000
   30 fbsi0 = xinf
      ier = 129
 9000 continue
c     stop
 9005 return
      end
c----------------------------------------------------------ok
      double precision function fbsi1 (iopt,arg,ier)
      implicit double precision(a-h,o-z)
      dimension p(15),pp(8),q(5),qq(6)
      data xcon/85.0d0/
      data xmax1 /.18014398509481984d17/
      data xmax  /91.200d0/
      data dexp40/ .235385266837019985d+18/
      data xsmall/8.67361738d-19/
      data xinf  /1.7d+38/
      data const / .666666666666666667d-01/
      data xmins/1.058791184d-22/
      data p /-.197052918025351399d-18,-.652455155831519027d-15,
     *        -.119287889036032387d-11,-.148319049359946477d-08,
     *        -.134668298276351529d-05,-.917464432878175013d-03,
     *        -.472070908273101624d+00,-.182259466316573159d+03,
     *        -.518940919823080174d+05,-.105885507247693471d+08,
     *        -.148282676066123661d+10,-.133574376822754930d+12,
     *        -.698767796480100898d+13,-.177320378407915915d+15,
     *        -.145771802781434636d+16/
      data q /-.400768646799041899d+04, .748105803566550690d+07,
     *        -.800595189986197647d+10, .485447142582736225d+13,
     *        -.132181683073214422d+16/
      data pp/-.604371590561376000d-01, .457481229019334590d+00,
     *        -.428437669033048064d+00, .973560001508866121d-01,
     *        -.324577239744655683d-02,-.363952647121217953d-03,
     *         .162586618674408364d-04,-.363475784046082235d-06/
      data qq/-.388065867215565934d+01, .325937148890369963d+01,
     *        -.850174764632179244d+00, .742120108131865300d-01,
     *        -.228356244894925126d-02, .375104331119228246d-04/
      data pbar  / .398437500000000000d+00/
      ier = 0
      if (iopt.ne.1.and.iopt.ne.2) go to 30
      x = dabs(arg)
      if (x.lt.xsmall) go to 20
      if ((iopt.eq.1).and.(x.gt.xmax)) go to 25
      if (x.ge.15.0d0) go to 10
      xx = x*x
      sump = p(1)
      do 5 j=2,15
         sump = sump*xx+p(j)
    5 continue
      xx = xx-225.0d0
      sumq = ((((xx+q(1))*xx+q(2))*xx+q(3))*xx+q(4))*xx+q(5)
      fbsi1 = (sump/sumq)*x
      if (iopt.eq.2) fbsi1 = fbsi1*dexp(-x)
      go to 35
   10 xx = 1.0d0/x-const
      sump = ((((((pp(1)*xx+pp(2))*xx+pp(3))*xx+pp(4))*xx+pp(5))*xx+
     *pp(6))*xx+pp(7))*xx+pp(8)
      sumq = (((((xx+qq(1))*xx+qq(2))*xx+qq(3))*xx+qq(4))*xx+qq(5))*xx+
     *qq(6)
      fbsi1 = (sump/sumq+pbar)/dsqrt(x)
      if (iopt.eq.2) go to 35
      if (x.gt.xcon) go to 15
      fbsi1 = fbsi1*dexp(x)
      go to 35
   15 a = fbsi1*dexp(x-40.0d0)
      fbsi1 = a*dexp40
      go to 35
   20 if (x*xmax1.lt.xmins) x = 0.0d0
      fbsi1 = 0.5d0*x
      go to 35
   25 fbsi1 = xinf
      ier = 130
      go to 35
   30 ier = 129
      fbsi1 = xinf
   35 if (arg.lt.0.0d0) fbsi1 = -fbsi1
 9000 continue
c     stop
 9005 return
      end
c----------------------------------------------------------ok
      double precision function fbsk0(iopt,arg,ier)
      implicit double precision(a-h,o-z)
      dimension ismall(2),f(4),g(3),p(6),pp(10),q(2),qq(10)
      equivalence        (xsmall,ismall(1))
      data const / .115931515658412449d+00/
      data xinf/1.7d+38/
      data xmax  /87.40d0/
      data xsmall/2.117582368d-22/
      data p / .585992214128261000d-03, .131660525649895718d+00,
     *         .119994637249107141d+02, .468509012019348321d+03,
     *         .591690598522705122d+04, .247081527203995526d+04/
      data q /-.249944189728323036d+03, .213127143038491203d+05/
      data f /-.164144528372990641d+01,-.296016578929588438d+03,
     *        -.177337846849529859d+05,-.403203407611454822d+06/
      data g /-.250649724458779927d+03, .298657131630540254d+05,
     *        -.161281363044581940d+07/
      data pp/ .113949805573847782d+03, .368325899573402679d+04,
     *         .310754089806843923d+05, .105770689480340219d+06,
     *         .173988679025656862d+06, .150976463532899145d+06,
     *         .715570627837640375d+05, .183215258701835377d+05,
     *         .234447387641993150d+04, .116002494250760355d+03/
      data qq/ .200134430649492424d+03, .443296288897464087d+04,
     *         .314746557502952788d+05, .974188297622680758d+05,
     *         .151446446735201578d+06, .126898395879775987d+06,
     *         .588246167858570276d+05, .148472283718023609d+05,
     *         .188218908409827137d+04, .925565991773048397d+02/
      ier = 0
      if (iopt.ne.1.and.iopt.ne.2) go to 40
      x = arg
      if (x.le.0.0d0) go to 40
      go to (5,10), iopt
    5 if (x.gt.xmax) go to 45
      if (x.le.1.0d0) go to 15
      go to 25
   10 if (x.gt.1.0d0) go to 25
   15 temp = dlog(x)
      if (x.lt.xsmall) go to 20
      xx = x*x
      sump = ((((p(1)*xx+p(2))*xx+p(3))*xx+p(4))*xx+p(5))*xx+p(6)
      sumq = (xx+q(1))*xx+q(2)
      sumf = ((f(1)*xx+f(2))*xx+f(3))*xx+f(4)
      sumg = ((xx+g(1))*xx+g(2))*xx+g(3)
      fbsk0 = sump/sumq-xx*sumf*temp/sumg-temp
      if (iopt.eq.2) fbsk0 = fbsk0*dexp(x)
      go to 9005
   20 fbsk0 = const-temp
      go to 9005
   25 xx = 1.0d0/x
      sump = pp(1)
      do 30 i=2,10
         sump = sump*xx+pp(i)
   30 continue
      sumq = xx
      do 35 i=1,9
         sumq = (sumq+qq(i))*xx
   35 continue
      sumq = sumq+qq(10)
      fbsk0 = sump/sumq/dsqrt(x)
      if (iopt.eq.1) fbsk0 = fbsk0*dexp(-x)
      go to 9005
   40 ier = 129
      fbsk0 = xinf
      go to 9000
   45 ier = 130
      fbsk0 = 0.0d0
 9000 continue
c     stop
 9005 return
      end
c-------------------------------------------------------ok
      double precision function fbsk1 (iopt,arg,ier)
      implicit double precision(a-h,o-z)
      dimension f(5),g(3),p(5),pp(11),q(3),qq(9)
      data xmax/86.7d0/
      data xinf /1.7d+38/
      data xmins/1.058791184d-22/
      data xsmall/1.734723476d-18/
      data xmax1/.1801439850948198d+17/
      data p / .481270704568784423d+00, .999913735674293098d+02,
     *         .718853826040847984d+04, .177333240351470156d+06,
     *         .719389200654205860d+06/
      data q /-.281439157545387258d+03, .372642986720676978d+05,
     *        -.221493748782433045d+07/
      data f /-.227955908269550024d+00,-.531039133351802751d+02,
     *        -.450516237634360869d+04,-.147580692054142224d+06,
     *        -.135311614927854213d+07/
      data g /-.305071515787875957d+03, .431176532113510799d+05,
     *        -.270623229855708426d+07/
      data pp/ .642577458591731388d-01, .755845846311760307d+01,
     *         .131826099185699413d+03, .810942561465374021d+03,
     *         .231237422091688715d+04, .345406755855445843d+04,
     *         .285906576979102882d+04, .133194864331832220d+04,
     *         .341229534868013129d+03, .441371761142304139d+02,
     *         .221967924968745489d+01/
      data qq/ .360010693068615188d+02, .330310200887653908d+03,
     *         .120826923160023486d+04, .211810004871719438d+04,
     *         .194484407889180061d+04, .969291657268026486d+03,
     *         .259512236555790513d+03, .345522284527589128d+02,
     *         .177104780326010866d+01/
      ier = 0
      if (iopt.ne.1.and.iopt.ne.2) go to 40
      x = arg
      go to (5,10), iopt
    5 if (x.gt.xmax) go to 45
      if (x.le.1.0d0) go to 15
      go to 25
   10 if (x.gt.1.0d0) go to 25
   15 if (x.lt.xsmall) go to 20
      xx = x*x
      sump = ((((p(1)*xx+p(2))*xx+p(3))*xx+p(4))*xx+p(5))*xx+q(3)
      sumq = ((xx+q(1))*xx+q(2))*xx+q(3)
      sumf = (((f(1)*xx+f(2))*xx+f(3))*xx+f(4))*xx+f(5)
      sumg = ((xx+g(1))*xx+g(2))*xx+g(3)
      fbsk1 = (xx*dlog(x)*sumf/sumg+sump/sumq)/x
      if (iopt.eq.2) fbsk1 = fbsk1*dexp(x)
      go to 9005
   20 if (x*xmax1.lt.xmins) go to 40
      fbsk1 = 1.0d0/x
      go to 9005
   25 xx = 1.0d0/x
      sump = pp(1)
      do 30 i=2,11
         sump = sump*xx+pp(i)
   30 continue
      sumq = xx
      do 35 i=1,8
         sumq = (sumq+qq(i))*xx
   35 continue
      sumq = sumq+qq(9)
      fbsk1 = sump/sumq/dsqrt(x)
      if (iopt.eq.1) fbsk1 = fbsk1*dexp(-x)
      go to 9005
   40 ier = 129
      fbsk1 = xinf
      go to 9000
   45 ier = 130
      fbsk1 = 0.0d0
 9000 continue
c     stop
 9005 return
      end
c----------------------------------------------------ok
      double precision function fbsj0(arg,ier)
      integer            ier
      double precision   arg
      double precision   ax,down,pi2,prod,prp,prq,
     1                   p0,p1,qrp,qrq,q0,q1,r0,r1,twopi1,twopi2,u,up,
     2                   w,xden,xmax,xnum,xsmall,x01,x02,x11,
     3                   x12,z,zsq
      dimension          prp(7),prq(8),p0(6),p1(6),qrp(5),qrq(7),
     1                   q0(6),q1(6)
      data pi2   / .636619772367581343d+00/
      data xmax  / .429496729600000000d+10/
      data xsmall/2.328306437d-10/
      data               x01  /    0.2404825557695772d 01/
      data               x02  /    0.1043754397719454d-15/
      data               x11  /    0.5520078110286310d 01/
      data               x12  /    0.8088597146146419d-16/
      data               twopi1/   0.6283185005187988d 01/
      data               twopi2/   0.3019915981956752d-06/
      data prp   / .663029979048337941d+07,-.621407004235401206d+09,
     *         .272825078786059427d+11,-.412986685009908667d+12,
     *        -.121170361645935283d+00, .103442228154431889d+03,
     *        -.366298146551070864d+05/
      data qrp   / .456126962242199381d+06, .139850973722634333d+09,
     *         .263281983008596486d+11, .238837879963322903d+13,
     *         .936140223923377105d+03/
      data prq   / .441767070253250875d+04, .117250462797571036d+05,
     *         .103419106415837267d+05,-.728797024644646189d+04,
     *        -.122540781613789895d+05,-.183193979693920840d+04,
     *         .485917033559164992d+02, .743211966806242457d+03/
      data qrq   / .333073107746490711d+03,-.294587665455093373d+04,
     *         .186809900083591883d+05,-.840550625911695622d+05,
     *         .245991022625863089d+06,-.357834780261523010d+06,
     *        -.252580762408015550d+02/
      data p0    / .153982653262391147d+06, .201613528304998364d+07,
     *         .841304145655043921d+07, .123323847681763814d+08,
     *         .539348508386943832d+07, .248527192895740401d+04/
      data q0    / .156001727694003094d+06, .202506680157013401d+07,
     *         .842644905062979733d+07, .123383102278632496d+08,
     *         .539348508386943832d+07, .261570073692083968d+04/
      data p1    /-.226263064193370411d+03,-.236595617077910818d+04,
     *        -.823906631348560656d+04,-.103814169874846409d+05,
     *        -.398461735759522246d+04,-.488719939584126153d+01/
      data q1    / .157048919151539552d+05, .156021320667929165d+06,
     *         .533291363421689717d+06, .666745423931982698d+06,
     *         .255015510886094237d+06, .408771467398349922d+03/
      ier = 0
      ax = dabs(arg)
      if (ax.gt.xmax) go to 35
      if (ax.gt.xsmall) go to 5
      fbsj0 = 1.0d0
      go to 9005
    5 if (ax.gt.8.0d0) go to 25
      if (ax.gt.4.0d0) go to 15
      zsq = ax*ax
      xnum = (prp(5)*zsq+prp(6))*zsq+prp(7)
      xden = zsq+qrp(5)
      do 10 i=1,4
         xnum = xnum*zsq+prp(i)
         xden = xden*zsq+qrp(i)
   10 continue
      prod = ((ax-x01)-x02)*(ax+x01)
      fbsj0 = prod*xnum/xden
      go to 9005
   15 zsq = 1.0d0-(ax*ax)/64.0d0
      xnum = prq(7)*zsq+prq(8)
      xden = zsq+qrq(7)
      do 20 i=1,6
         xnum = xnum*zsq+prq(i)
         xden = xden*zsq+qrq(i)
   20 continue
      prod = (ax+x11)*((ax-x11)-x12)
      fbsj0 = prod*xnum/xden
      go to 9005
   25 z = 8.0d0/ax
      w = ax/twopi1
      w = dint(w) + 0.125d0
      u = (ax-w*twopi1)-w*twopi2
      zsq = z*z
      xnum = p0(6)
      xden = zsq+q0(6)
      up = p1(6)
      down = zsq+q1(6)
      do 30 i=1,5
         xnum = xnum*zsq+p0(i)
         xden = xden*zsq+q0(i)
         up = up*zsq+p1(i)
         down = down*zsq+q1(i)
   30 continue
      r0 = xnum/xden
      r1 = up/down
      fbsj0 = dsqrt(pi2/ax)*(r0*dcos(u)-z*r1*dsin(u))
      go to 9005
   35 fbsj0 = 0.0d0
      ier = 129
 9000 continue
 9005 return
      end
c-----------------------------------------------------------------------
      subroutine fbsyn (arg,order,n,yn,ier)
      integer            n,ier
      double precision   arg,order,yn(n)
      integer            i,j,k,l,m,nmax,nn   !itwpi1(2),itwpi2(2),
      integer            iend,mend,m1,m2,m3,iendp2,mendm3
      double precision   co,cum,cump,c2,c2p,delx,den,denp,
     1                   dlnhaf,dlnx,dlnxnu,fact(34),foolu,f1(17),
     2                   f2(16),f3(17),gam,gnu,gnu1,gnu2,gp(8),gq(8),
     3                   half,p(12,7),pi2,q(11,7),
     5                   twopi1,twopi2,t1,t2,t3(6),v,vcos,vpi,vsin,
     4                   rtpi2,rtxnu,s,scrat(6),sum,sump,t,twopi,
	6                   w(19),x,xc,xin,xinf,xk,xm,xnu,xnum,xnump, 
     7                   z,zh,zmax,zmin,zstop,con1,con2
      double precision   twopia,twopib,zmax2
      logical            nflag
      data pi2   / .636619772367581343d+00/
      data rtpi2 / .797884560802865356d+00/
      data dlnhaf/-.693147180559945310d+00/
      data twopi1/0.62831850051879880d 01/
      data twopi2/0.30199159819567520d-06/
      data twopia/0.62812500000000000d 01/
      data twopib/0.19350051879882810d-02/
      data twopi / .628318530717958647d+01/
      data half  / .500000000000000000d+00/
      data zmax2 / 4294967296.0d0/
      data zmin  / .294d-38/
      data xinf  / 1.7d+38/
      data gnu1  / .500488281250000000d-01/
      data gnu2  / .996093750000000000d+00/
      data zmax  / .429496729600000000d+10/
      data co    / .886226925452758014d+00/
      data gp    /-.160243404626643007d+01, .148472957685686770d+02,
     *        -.280750941362157079d+03,-.417571326141158257d+03,
     *         .594221994599142643d+03,-.273925800230766587d+05,
     *        -.778326730220496473d+05,-.577828296022400365d+04/
      data gq    /-.178681796751286802d+06,-.268406415244627174d+02,
     *         .214414834924538890d+03,-.224026496383528068d+03,
     *        -.459374814623425838d+04, .145077504178693444d+05,
     *         .329089709308774549d+05,-.114818917860301046d+06/
      data p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),
     *     p(8,1),p(9,1),p(10,1),p(11,1),p(12,1)
     *           / .966220363059263586d+02,-.135166484200055635d+04,
     *         .864912264696622203d+04,-.265430829974670190d+05,
     *        -.258277503715144970d+04, .253958451563109432d+06,
     *        -.741666962661522701d+06,-.185303404753455436d+06,
     *         .249424577877157393d+07,-.328627552380477017d+07,
     *        -.588915321642776927d+07, .456049615140191762d+06/
      data p(1,2),p(2,2),p(3,2),p(4,2),p(5,2),p(6,2),p(7,2),
     *     p(8,2),p(9,2),p(10,2)
     *           / .827836377496675271d+02,-.126469763244904028d+04,
     *         .580689617323250046d+04, .212643684799830998d+05,
     *        -.229280692127117227d+06, .172961254163067550d+06,
     *         .265359659939117365d+07,-.205519713780387999d+07,
     *        -.834208771374013572d+07, .228340012506985874d+06/
      data p(1,3),p(2,3),p(3,3),p(4,3),p(5,3),p(6,3),p(7,3),
     *     p(8,3),p(9,3),p(10,3)
     *           /-.655464718390655217d+02, .164166175577358534d+04,
     *        -.106542820607880074d+05,-.319526554910620985d+05,
     *         .488095894287561287d+06,-.108912380388916879d+06,
     *        -.626824801218970995d+07, .138102270261905953d+07,
     *         .184668308157951082d+08, .100645521111378579d+07/
      data p(1,4),p(2,4),p(3,4),p(4,4),p(5,4),p(6,4),p(7,4),
     *     p(8,4),p(9,4),p(10,4)
     *           / .696963123814590427d+02,-.133020689714923170d+04,
     *         .452303293370654331d+04, .471702732719340836d+05,
     *        -.260746560822329284d+06,-.494295070906759632d+06,
     *         .333646692280230407d+07, .169164170133818384d+07,
     *        -.948024850300566443d+07,-.136168612925693158d+07/
      data p(1,5),p(2,5),p(3,5),p(4,5),p(5,5),p(6,5),p(7,5),
     *     p(8,5),p(9,5),p(10,5)
     *           /-.887596516599601438d+02, .145281203139013753d+04,
     *        -.250210995757854514d+04,-.600773443105863251d+05,
     *         .213021520919733862d+06, .770892892749782111d+06,
     *        -.295732940068198094d+07,-.311681822772432647d+07,
     *         .869420380824761184d+07, .211470077845670982d+07/
      data p(1,6),p(2,6),p(3,6),p(4,6),p(5,6),p(6,6),p(7,6),
     *     p(8,6),p(9,6),p(10,6)
     *           / .105141263149205264d+03,-.161405265322962398d+04,
     *         .154355389910810923d+04, .728265946613332423d+05,
     *        -.208095986994592776d+06,-.104246097148516522d+07,
     *         .312777482643423898d+07, .474231617027130897d+07,
     *        -.951912134902371432d+07,-.338528003709295307d+07/
      data p(1,7),p(2,7),p(3,7),p(4,7),p(5,7),p(6,7),p(7,7),
     *     p(8,7),p(9,7),p(10,7),p(11,7)
     *           / .336889850414961424d+02,-.422443710176094398d+03,
     *        -.203831469281181514d+04, .331030626094266419d+05,
     *         .604252878602211945d+05,-.933840931715176502d+06,
     *        -.113618427457866563d+07, .959057577413757707d+07,
     *         .903854881191026607d+07,-.236832975247582365d+08,
     *        -.104340643565508042d+08/
      data q(1,1),q(2,1),q(3,1),q(4,1),q(5,1),q(6,1),q(7,1),
     *     q(8,1),q(9,1),q(10,1),q(11,1)
     *           /-.101241970477141194d+02, .216005595832218891d+02,
     *        -.102815919664190652d+03, .114516150295087211d+04,
     *        -.433090914101065738d+03, .460732511990055463d+03,
     *        -.851274820032285445d+05,-.768523228710292092d+05,
     *         .205560195496649931d+06, .364590804564888233d+07,
     *         .516729324640858565d+07/
      data q(1,2),q(2,2),q(3,2),q(4,2),q(5,2),q(6,2),q(7,2),
     *     q(8,2),q(9,2)
     *           / .225647845389957788d+01, .711593716504863289d+02,
     *        -.141959532218571208d+03,-.609226713172482790d+03,
     *        -.331287173625895019d+05,-.179090061095354737d+06,
     *        -.134472863756067396d+07,-.392640605688557498d+07,
     *        -.134787513180809696d+08/
      data q(1,3),q(2,3),q(3,3),q(4,3),q(5,3),q(6,3),q(7,3),
     *     q(8,3),q(9,3)
     *           /-.939317088650692602d+01, .523806373239931005d+02,
     *        -.121057364184354667d+04,-.258682317593112574d+04,
     *        -.844920979601552086d+05,-.292487285966100229d+06,
     *        -.301034666110382428d+07,-.713311890076026399d+07,
     *        -.387847838647490074d+08/
      data q(1,4),q(2,4),q(3,4),q(4,4),q(5,4),q(6,4),q(7,4),
     *     q(8,4),q(9,4)
     *           /-.437886178371543549d+01, .108803383282850214d+03,
     *        -.396170247546345465d+03, .482799277343449138d+04,
     *        -.260956190444154348d+05, .795561217047195187d+05,
     *        -.121501380828845716d+07,-.560365677885192256d+06,
     *        -.244594499396832557d+08/
      data q(1,5),q(2,5),q(3,5),q(4,5),q(5,5),q(6,5),q(7,5),
     *     q(8,5),q(9,5)
     *           /-.581690214590269560d+01, .121643903355730538d+03,
     *        -.519952651037368415d+03, .634791999150319729d+04,
     *        -.291287737773905651d+05, .162626855091396249d+06,
     *        -.124038130344249191d+07, .121060159492453767d+07,
     *        -.270394900976098049d+08/
      data q(1,6),q(2,6),q(3,6),q(4,6),q(5,6),q(6,6),q(7,6),
     *     q(8,6),q(9,6)
     *           /-.788702380969431594d+01, .130637206116758225d+03,
     *        -.771892722230546718d+03, .715457420628860810d+04,
     *        -.431226753881803730d+05, .199060268477405625d+06,
     *        -.169463010448782541d+07, .196130269020172899d+07,
     *        -.353328964340194695d+08/
      data q(1,7),q(2,7),q(3,7),q(4,7),q(5,7),q(6,7),q(7,7),
     *     q(8,7),q(9,7)
     *           / .359926636077034681d+01, .144952858594130714d+03,
     *         .998518492990023872d+03, .121000936522091880d+05,
     *         .998406746198100760d+05, .637436088274496556d+06,
     *         .488124457435495168d+07, .158873416545813497d+08,
     *         .952711263569843984d+08/
      data f1    / .526129484784613088d-05, .313533730646244067d-05,
     *         .468966778832002960d-04, .387395004571051596d-04,
     *        -.433600138305899192d-03, .233918540353778055d-02,
     *        -.840830427176014116d-03,-.210619747888291847d-01,
     *         .747920302081092203d-01,-.603883559512617772d-02,
     *        -.417981294786837086d+00, .668645110914854594d+00,
     *         .638141340718330557d+00,-.143970347809096185d+01,
     *        -.134576433585482691d+00, .318309886183790672d+00,
     *         .258709136016583669d-05/
      data f2    / .572867669670037890d-06, .229397311396353895d-04,
     *         .209873686010719157d-04, .143717132163052397d-03,
     *        -.510284654371587433d-03, .213837479821310760d-02,
     *         .104955004042030943d-02,-.209403929898152524d-01,
     *         .521523272704097449d-01, .732116630026647386d-01,
     *        -.427655556436905675d+00, .274382993322142905d+00,
     *         .130911899735173940d+01,-.696436028549460850d+00,
     *        -.157079632679489662d+01, .367466905196615962d+00/
      data f3    /-.484943071073704969d-05, .969366827458976827d-05,
     *        -.193509783140010526d-04, .387299755087458891d-04,
     *        -.766940410140475115d-04, .156095491967178727d-03,
     *        -.292709443974278280d-03, .669702782215668059d-03,
     *        -.908025366654622614d-03, .355044302302707344d-02,
     *        -.849830319153386275d-04, .236341941622100586d-01,
     *         .259667398807643121d-01, .131092848703934761d+00,
     *         .134576433585482691d+00, .318309886183790672d+00,
     *         .242608782760487999d-05/
      data fact  / .100000000000000000d+01, .100000000000000000d+01,
     *         .500000000000000000d+00, .166666666666666667d+00,
     *         .416666666666666667d-01, .833333333333333333d-02,
     *         .138888888888888889d-02, .198412698412698413d-03,
     *         .248015873015873016d-04, .275573192239858906d-05,
     *         .275573192239858906d-06, .250521083854417187d-07,
     *         .208767569878680990d-08, .160590438368216146d-09,
     *         .114707455977297247d-10, .764716373181981647d-12,
     *         .477947733238738529d-13, .281145725434552076d-14,
     *         .156192069685862264d-15, .822063524662432970d-17,
     *         .411031762331216485d-18, .195729410633912612d-19,
     *         .889679139245057327d-21, .386817017063068404d-22,
     *         .161173757109611835d-23, .644695028438447338d-25,
     *         .247959626322479746d-26, .918368986379554613d-28,
     *         .327988923706983791d-29, .113099628864477169d-30,
     *         .376998762881590564d-32, .121612504155351795d-33,
     *         .380039075485474359d-35, .115163356207719503d-36/
      zstop = 20.5d0
      con1 = 25.0d0
      con2 = 35.0d0
      iend = 17
      mend = 16
      m1 = 11
      m2 = 8
      m3 = 4
      iendp2 = iend+2
      mendm3 = mend-3
      z = arg
      ier = 0
      v = order
      nflag = (n.eq.1)
      iopt = 2
      if (nflag) iopt = 1
      if (z.ge.zmin) go to 5
      ier = 129
    5 if ((v.ge.0.0d0).and.(v.lt.1.0d0)) go to 10
      ier = 129
   10 if (n.gt.0) go to 15
      ier = 130
   15 if (ier.eq.0) go to 25
      j = n
      if (ier.eq.130) j = 1
      do 20 i=1,j
         yn(i) = -xinf
   20 continue
      go to 9000
   25 if (v.lt.zmin) v = 0.0d0
      if (z.lt..03125d0) go to 70
      if (z.ge.zstop) go to 145
      j = idint((z+3.5d0)/3.0d0)
      x = dble(float(j))*3.0d0-2.0d0
      m = 12
      if (j.gt.1) m = 10
      gnu = v+1.0d0
      xnum = p(1,j)*v+p(2,j)
      xnump = xnum+p(1,j)
      den = v+q(1,j)
      denp = den+1.0d0
      do 30 i=3,m
         xnum = xnum*v+p(i,j)
         den = den*v+q(i-1,j)
         xnump = xnump*gnu+p(i,j)
         denp = denp*gnu+q(i-1,j)
   30 continue
      if (j.ne.7) go to 35
      xnum = xnum*v+p(11,j)
      xnump = xnump*gnu+p(11,j)
   35 delx = .0625d0
      if (z.gt.4.0d0) delx = delx+delx
      t = z-x
      if (t.lt.0.0d0) delx = -delx
      m = mend
      if (z.gt.4.0d0) m = mendm3
      xnu = 1.0d0
      if (x.ne.1.0d0) xnu = (x**v)*x
      w(1) = xnump/denp*xnu
      w(2) = xnum/den*xnu
      xnu = 2.0d0*(gnu-half)
      if (z.ge..5d0) go to 40
      t = z+z-1.5d0
      if (z.ge..25d0) go to 40
      t = 4.0d0*z-2.0d0
      if (z.lt..125d0) t = 8.0d0*z-2.5d0
   40 i = idint(t/delx)+1
      c2 = 0.0d0
      c2p = 0.0d0
      cum = w(1)
      cump = w(2)
      nmax = m+3
      do 65 nn=1,i
         if (nn.lt.i) go to 45
         delx = z-x
         if (delx.eq.0.0d0) go to 65
   45    if (z.ge.1.0d0) go to 50
         if (nn.eq.9) delx = dmax1(delx,-.03125d0)
         if (nn.eq.17) delx = dmax1(delx,-.015625d0)
         if (nn.eq.25) delx = dmax1(delx,-.0078125d0)
   50    xin = 1.0d0/x
         w(3) = xnu*xin*w(2)-w(1)
         xk = 1.0d0
         s = xnu
         do 55 k=4,nmax
            s = s-1.0d0
            w(k) = (s*w(k-1)-xk*w(k-3))*xin-w(k-2)
            xk = xk+1.0d0
   55    continue
         sum = fact(m+2)*w(m+2)*delx
         sump = fact(m+2)*w(m+3)*delx
         j = m+1
         do 60 k=1,m
            sum = (sum+fact(j)*w(j))*delx
            sump = (sump+fact(j)*w(j+1))*delx
            j = j-1
   60    continue
         c2 = c2+sum
         t = cum+c2
         foolu = cum-t
         c2 = c2+foolu
         cum = t
         c2p = c2p+sump
         t = cump+c2p
         foolu = cump-t
         c2p = c2p+foolu
         cump = t
         x = x+delx
         w(1) = cum
         w(2) = cump
   65 continue
      xnu = (z**v)*z
      yn(1) = w(2)/xnu
      if (nflag) go to 9005
      yn(2) = w(1)/xnu
      go to 170
   70 zh = z*0.5d0
      x = zh*zh
      if ((v.lt.gnu1).or.(v.gt.gnu2)) go to 95
      x = -x
      gnu = v
      vpi = v/pi2
      if (v.gt.0.5d0) vpi = (1.0d0-v)/pi2
      vpi = vpi+vpi
      vsin = dsin(vpi)
      vcos = dcos(vpi)
      if (v.gt.0.5d0) vcos = -vcos
      t1 = zh**gnu
      t2 = 1.0d0/t1
      do 90 i=1,2
         xc = t1
         do 85 j=1,2
            t = gnu-0.5d0
            if (j.eq.2) t = t+2.0d0
            if (t.gt.0.5d0) t = t-1.0d0
            xnum = gp(1)*t
            den = t
            do 75 l=2,8
               xnum = (xnum+gp(l))*t
               den = (den+gq(l))*t
   75       continue
            den = den+gq(1)
            gam = co+xnum/den
            if (j.lt.i) gam = gam*gnu
            if (j.eq.2) gam = gam/(gnu+1.0d0)
            if (i*j.eq.4) gam = gam/(gnu+2.0d0)
            sum = fact(5)
            xk = 4.0d0
            do 80 k=1,4
               sum = sum*x/(gnu+xk)+fact(5-k)
               xk = xk-1.0d0
   80       continue
            scrat(j) = xc*sum/gam
            gnu = -gnu
            xc = t2
   85    continue
         yn(i) = (vcos*scrat(1)-scrat(2))/vsin
         if (nflag) go to 9005
         if (i.eq.2) go to 170
         vcos = -vcos
         vsin = -vsin
         gnu = v+1.0d0
         t1 = t1*zh
         t2 = t2/zh
   90 continue
      go to 170
   95 gnu = v
      if (v.gt.gnu2) gnu = (v-half)-0.5d0
      t1 = f1(17)
      t2 = 0.0d0
      t3(1) = f3(17)
      do 100 i=1,16
         t1 = t1*gnu+f1(i)
         t2 = t2*gnu+f2(i)
         t3(1) = t3(1)*gnu+f3(i)
  100 continue
      t1 = t1/(1.0d0-gnu)
      t2 = t2/(1.0d0-gnu*gnu)
      t3(1) = t3(1)/(1.0d0+gnu)
      dlnx = 2.0d0*dlog(zh)
      dlnxnu = gnu*dlnx
      if (dlnxnu.lt.dlnhaf) go to 110
      xnu = fact(iendp2)
      do 105 i=1,iend
         xnu = (xnu*dlnxnu)+fact(iendp2-i)
  105 continue
      xnu = dlnx*xnu
      scrat(1) = t2+xnu*t1
      xnu = gnu*xnu+1.0d0
      go to 115
  110 xnu = zh**(gnu+gnu)
      scrat(1) = (xnu*t1-t3(1))/gnu
  115 xk = 1.0d0
      do 120 i=2,6
         t = xk-gnu
         t3(i) = t3(i-1)/(t*xk)
         scrat(i) = (scrat(i-1)*(1.0d0-gnu/xk)/t-2.0d0*t3(i))/(xk+gnu)
         xk = xk+1.0d0
  120 continue
      xk = 5.0d0+gnu
      sum = scrat(5)
      sump = xk*scrat(6)+t3(6)
      if (x.lt.1.0d-30) x = 0.0d0
      do 125 i=1,4
         k = 5-i
         xk = xk-1.0d0
         sump = (-sump)*x+xk*scrat(k+1)+t3(k+1)
         sum = sum*(-x)+scrat(k)
  125 continue
      sump = xnu*t1/zh-sump*zh
      rtxnu = dsqrt(xnu)
      xnu = gnu/zh
      sum = sum/rtxnu
      if (gnu.lt.0.0d0) go to 130
      yn(1) = sum
      if (nflag) go to 9005
      yn(2) = xnu*yn(1)-sump/rtxnu
      go to 170
  130 yn(1) = xnu*sum-sump/rtxnu
      if (nflag) go to 9005
      xnu = xnu+2.0d0/z
      if (xnu.ge.-xinf/yn(1)) go to 135
      yn(2) = xnu*yn(1)-sum
      go to 170
  135 do 140 i=2,n
         yn(i) = -xinf
  140 continue
      ier = 131
      go to 9000
  145 if (z.le.zmax) go to 155
      ier = 132
      do 150 i=1,n
         yn(i) = 0.0d0
  150 continue
      go to 9000
  155 xc = rtpi2/dsqrt(z)
      xin = .015625d0/(z*z)
      m = 16
      if (z.ge.con1) m = m1
      if (z.ge.con2) m = m2
      if (z.ge.130.0d0) m = m3
      xm = 4.0d0*dble(float(m))
      t = dint(z/twopi)
      if (t.le.zmax2) go to 158
      x = (((z-t*twopia)-t*twopib)-t*twopi2)-(v+0.5d0)/pi2
      go to 159
  158 x = ((z-t*twopi1)-t*twopi2)-(v+0.5d0)/pi2
  159 vsin = dsin(x)
      vcos = dcos(x)
      gnu = v+v
      do 165 i=1,2
         s = ((xm-1.0d0)-gnu)*((xm-1.0d0)+gnu)*xin*0.5d0
         t = (gnu-(xm-3.0d0))*(gnu+(xm-3.0d0))
         xnum = s*t*fact(2*m+1)
         t1 = (gnu-(xm+1.0d0))*(gnu+(xm+1.0d0))
         den = s*t1*fact(2*m+2)
         xk = xm
         k = m+m
         t1 = t
         do 160 j=2,m
            xk = xk-4.0d0
            s = ((xk-1.0d0)-gnu)*((xk-1.0d0)+gnu)
            t = (gnu-(xk-3.0d0))*(gnu+(xk-3.0d0))
            xnum = (xnum+fact(k-1))*s*t*xin
            den = (den+fact(k))*s*t1*xin
            k = k-2
            t1 = t
  160    continue
         xnum = xnum+1.0d0
         den = (den+1.0d0)*(gnu*gnu-1.0d0)/(8.0d0*z)
         yn(i) = xc*(xnum*vsin+den*vcos)
         if (nflag) go to 9005
         t = vsin
         vsin = -vcos
         vcos = t
         gnu = gnu+2.0d0
  165 continue
  170 if (n.le.2) go to 9005
      gnu = v+v+2.0d0
      nn = n
      do 185 j=3,nn
         if (dabs(yn(j-1)).lt.1.0d0) go to 180
         if (dabs(gnu/z).lt.dabs(xinf/yn(j-1))) go to 180
         do 175 l=j,n
            yn(l) = -xinf
  175    continue
         ier = 131
         go to 9000
  180    yn(j) = gnu*yn(j-1)/z-yn(j-2)
         gnu = gnu+2.0d0
  185 continue
      go to 9005
 9000 continue
 9005 return	 
      End	  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC - END
