subroutine dqmomo(alfa,beta,ri,rj,rg,rh,integr)
c***begin prologue  dqmomo
c***date written   820101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a2a1,c3a2
c***keywords  modified chebyshev moments
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
c***purpose  this routine computes modified chebsyshev moments. the k-th
c            modified chebyshev moment is defined as the integral over
c            (-1,1) of w(x)*t(k,x), where t(k,x) is the chebyshev
c            polynomial of degree k.
c***description
c
c        modified chebyshev moments
c        standard fortran subroutine
c        double precision version
c
c        parameters
c           alfa   - double precision
c                    parameter in the weight function w(x), alfa.gt.(-1)
c
c           beta   - double precision
c                    parameter in the weight function w(x), beta.gt.(-1)
c
c           ri     - double precision
c                    vector of dimension 25
c                    ri(k) is the integral over (-1,1) of
c                    (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
c
c           rj     - double precision
c                    vector of dimension 25
c                    rj(k) is the integral over (-1,1) of
c                    (1-x)**beta*t(k-1,x), k = 1, ..., 25.
c
c           rg     - double precision
c                    vector of dimension 25
c                    rg(k) is the integral over (-1,1) of
c                    (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           rh     - double precision
c                    vector of dimension 25
c                    rh(k) is the integral over (-1,1) of
c                    (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           integr - integer
c                    input parameter indicating the modified
c                    moments to be computed
c                    integr = 1 compute ri, rj
c                           = 2 compute ri, rj, rg
c                           = 3 compute ri, rj, rh
c                           = 4 compute ri, rj, rg, rh
c
c***references  (none)
c***routines called  (none)
c***end prologue  dqmomo
c
      double precision alfa,alfp1,alfp2,an,anm1,beta,betp1,betp2,ralf,
     *  rbet,rg,rh,ri,rj
      integer i,im1,integr
c
      dimension rg(25),rh(25),ri(25),rj(25)
c
c
c***first executable statement  dqmomo
      alfp1 = alfa+0.1d+01
      betp1 = beta+0.1d+01
      alfp2 = alfa+0.2d+01
      betp2 = beta+0.2d+01
      ralf = 0.2d+01**alfp1
      rbet = 0.2d+01**betp1
c
c           compute ri, rj using a forward recurrence relation.
c
      ri(1) = ralf/alfp1
      rj(1) = rbet/betp1
      ri(2) = ri(1)*alfa/alfp2
      rj(2) = rj(1)*beta/betp2
      an = 0.2d+01
      anm1 = 0.1d+01
      do 20 i=3,25
        ri(i) = -(ralf+an*(an-alfp2)*ri(i-1))/(anm1*(an+alfp1))
        rj(i) = -(rbet+an*(an-betp2)*rj(i-1))/(anm1*(an+betp1))
        anm1 = an
        an = an+0.1d+01
   20 continue
      if(integr.eq.1) go to 70
      if(integr.eq.3) go to 40
c
c           compute rg using a forward recurrence relation.
c
      rg(1) = -ri(1)/alfp1
      rg(2) = -(ralf+ralf)/(alfp2*alfp2)-rg(1)
      an = 0.2d+01
      anm1 = 0.1d+01
      im1 = 2
      do 30 i=3,25
        rg(i) = -(an*(an-alfp2)*rg(im1)-an*ri(im1)+anm1*ri(i))/
     *  (anm1*(an+alfp1))
        anm1 = an
        an = an+0.1d+01
        im1 = i
   30 continue
      if(integr.eq.2) go to 70
c
c           compute rh using a forward recurrence relation.
c
   40 rh(1) = -rj(1)/betp1
      rh(2) = -(rbet+rbet)/(betp2*betp2)-rh(1)
      an = 0.2d+01
      anm1 = 0.1d+01
      im1 = 2
      do 50 i=3,25
        rh(i) = -(an*(an-betp2)*rh(im1)-an*rj(im1)+
     *  anm1*rj(i))/(anm1*(an+betp1))
        anm1 = an
        an = an+0.1d+01
        im1 = i
   50 continue
      do 60 i=2,25,2
        rh(i) = -rh(i)
   60 continue
   70 do 80 i=2,25,2
        rj(i) = -rj(i)
   80 continue
   90 return
      end