!========================存在奇异点的积分===================================
subroutine Ymean1(alpha0,beta0,k1,i1,j1)
use MODGLOB
use IMSL
USE INTEGRAL
implicit none
real*8 alpha0,beta0,es,sum0,sum01,sum02,sum03,width1,width2,width3,sum11,sum12,sum13,sum1,width11,width12,width13,&
            & sum3,sum00,sum10,width
INTEGER,PARAMETER :: N=81
real*8 f(12),x1(41),x2(20),x3(20),fun1(41),fun2(20),fun3(20),jf1(11),jf11(11)
real*8 f3(12),x31(41),x32(20),x33(20),fun31(41),fun32(20),fun33(20)
real*8,external:: func,fun
real datas(n)
integer ii,jj,k1,i1,j1
es=1.0d-6
f=(/es,1.0d-5,1.0d-4,1.0d-3,1.0d-2,0.1d0,0.9d0,1.0-1.0d-2,1.0-1.0d-3,1.0-1.0d-4,1.0-1.0d-5,1.0-es/)

do ii=1,12
  call enlg3(f0,yy(1:81,k1,i1,j1),81,f(ii),f3(ii))
enddo
sum0=es**alpha0/alpha0+es**beta0/beta0
sum1=es**alpha0/alpha0*yy(1,k1,i1,j1)+es**beta0/beta0*yy(81,k1,i1,j1)
sum00=0.0
sum10=0.0
if(ii==6) then
    width1=(f(ii+1)-f(ii))/40.d0
    width11=(f3(ii+1)-f3(ii))/40.d0
    do jj=1,41
    x1(jj)=f(ii)+width1*(jj-1)
    x31(jj)=f3(ii)+width11*(jj-1)
    fun1(jj)=x1(jj)**(alpha0-1.d0)*(1.d0-x1(jj))**(beta0-1.d0)
    fun31(jj)=x31(jj)*x1(jj)**(alpha0-1.d0)*(1.d0-x1(jj))**(beta0-1.d0)
    enddo
    call GenerateData(DATAS,WIDTH,FUNC)
    jf1(6)=Trape_Integral(fun1,width1)
    jf11(6)=Trape_Integral(fun31,width1)
    sum01=sum00+jf1(6)
    sum11=sum10+jf11(6)
else if(ii<=5) then
    width2=(f(ii+1)-f(ii))/20.d0
    width12=(f3(ii+1)-f3(ii))/20.d0
    do jj=1,20
    x2(jj)=f(ii)+width2*(jj-1)
    x32(jj)=f3(ii)+width12*(jj-1)
    fun2(jj)=x2(jj)**(alpha0-1.0)*(1.d0-x2(jj))**(beta0-1.d0)
    fun32(jj)=x32(jj)*x2(jj)**(alpha0-1.0)*(1.d0-x2(jj))**(beta0-1.d0)
    enddo
    do ii=1,5
    jf1(ii)=Trape_Integral(fun2,width2)
    jf11(ii)=Trape_Integral(fun32,width2)
    sum02=sum00+jf1(ii)
    sum12=sum10+jf11(ii)
    enddo
else
    width3=(f(ii+1)-f(ii))/20.d0
    width13=(f3(ii+1)-f3(ii))/20.d0
    do jj=1,20
    x3(jj)=f(ii)+width3*(jj-1)
    x33(jj)=f3(ii)+width13*(jj-1)
    fun3(jj)=x3(jj)**(alpha0-1.0)*(1.d0-x3(jj))**(beta0-1.d0)
    fun33(jj)=x33(jj)*x3(jj)**(alpha0-1.0)*(1.d0-x3(jj))**(beta0-1.d0)
    enddo
    do ii=6,11
    jf1(ii)=Trape_Integral(fun3,width3)
    jf11(ii)=Trape_Integral(fun33,width3)
    sum01=sum00+jf1(ii)
    sum11=sum10+jf11(ii)
    enddo 
endif
sum0=sum0+sum01+sum02+sum03
sum1=sum1+sum11+sum12+sum13
YYM(k1,i1,j1)=sum1/sum0
return
end subroutine
!=======================无奇异点的积分====================================
subroutine Ymean2(k1,i1,j1)
use MODGLOB
use IMSL
USE INTEGRAL
implicit none
real*8  sum,width1,width2,width
INTEGER,PARAMETER :: N=81
real datas(N)
real*8 x(11),fun(11),y1(11),p(11),jf(11)
real*8 ,external:: func
integer ii,jj,k1,i1,j1

sum=0.d0
do ii=1,80
  width=(f0(ii+1)-f0(ii))/10.d0
  width1=(yy(ii+1,k1,i1,j1)-yy(ii,k1,i1,j1))/10.d0
  width2=(pdf(i1,j1,ii+1)-pdf(i1,j1,ii))/10.d0
  do jj=1,11
    x(jj)=f0(ii)+width*(jj-1)
    y1(jj)=yy(ii,k1,i1,j1)+width1*(jj-1)
    p(jj)=pdf(i1,j1,ii)+width2*(jj-1)
    fun(jj)=y1(jj)*p(jj)
  enddo
  call GenerateData(datas,width,func)
  JF(jj)=Trape_Integral(fun,width)
  sum=sum+JF(jj)
enddo
yym(k1,i1,j1)=sum
return
end subroutine