[alfa,beta,ri,rj,rg,rh,integr]=qmomo(alfa,beta,ri,rj,rg,rh,integr);

function [alfa,beta,ri,rj,rg,rh,integr]=qmomo(alfa,beta,ri,rj,rg,rh,integr);
%***BEGIN PROLOGUE  QMOMO
%***PURPOSE  This routine computes modified Chebyshev moments.  The K-th
%            modified Chebyshev moment is defined as the integral over
%            (-1,1) of W(X)*T(K,X), where T(K,X) is the Chebyshev
%            polynomial of degree K.
%***LIBRARY   SLATEC (QUADPACK)
%***CATEGORY  H2A2A1, C3A2
%***TYPE      SINGLE PRECISION (QMOMO-S, DQMOMO-D)
%***KEYWORDS  MODIFIED CHEBYSHEV MOMENTS, QUADPACK, QUADRATURE
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%        MODIFIED CHEBYSHEV MOMENTS
%        STANDARD FORTRAN SUBROUTINE
%        REAL VERSION
%
%        PARAMETERS
%           ALFA   - Real
%                    Parameter in the weight function W(X), ALFA.GT.(-1)
%
%           BETA   - Real
%                    Parameter in the weight function W(X), BETA.GT.(-1)
%
%           RI     - Real
%                    Vector of dimension 25
%                    RI(K) is the integral over (-1,1) of
%                    (1+X)**ALFA*T(K-1,X), K = 1, ..., 25.
%
%           RJ     - Real
%                    Vector of dimension 25
%                    RJ(K) is the integral over (-1,1) of
%                    (1-X)**BETA*T(K-1,X), K = 1, ..., 25.
%
%           RG     - Real
%                    Vector of dimension 25
%                    RG(K) is the integral over (-1,1) of
%                    (1+X)**ALFA*LOG((1+X)/2)*T(K-1,X), K = 1, ..., 25.
%
%           RH     - Real
%                    Vector of dimension 25
%                    RH(K) is the integral over (-1,1) of
%                    (1-X)**BETA*LOG((1-X)/2)*T(K-1,X), K = 1, ..., 25.
%
%           INTEGR - Integer
%                    Input parameter indicating the modified
%                    Moments to be computed
%                    INTEGR = 1 compute RI, RJ
%                           = 2 compute RI, RJ, RG
%                           = 3 compute RI, RJ, RH
%                           = 4 compute RI, RJ, RG, RH
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   810101  DATE WRITTEN
%   891009  Removed unreferenced statement label.  (WRB)
%   891009  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  QMOMO
%
persistent alfp1 alfp2 an anm1 betp1 betp2 i im1 ralf rbet ; 

if isempty(alfp1), alfp1=0; end;
if isempty(alfp2), alfp2=0; end;
if isempty(an), an=0; end;
if isempty(anm1), anm1=0; end;
if isempty(betp1), betp1=0; end;
if isempty(betp2), betp2=0; end;
if isempty(ralf), ralf=0; end;
if isempty(rbet), rbet=0; end;
if isempty(i), i=0; end;
if isempty(im1), im1=0; end;
%
%
%
%***FIRST EXECUTABLE STATEMENT  QMOMO
alfp1 = alfa + 0.1e+01;
betp1 = beta + 0.1e+01;
alfp2 = alfa + 0.2e+01;
betp2 = beta + 0.2e+01;
ralf = 0.2e+01.^alfp1;
rbet = 0.2e+01.^betp1;
%
%           COMPUTE RI, RJ USING A FORWARD RECURRENCE RELATION.
%
ri(1) = ralf./alfp1;
rj(1) = rbet./betp1;
ri(2) = ri(1).*alfa./alfp2;
rj(2) = rj(1).*beta./betp2;
an = 0.2e+01;
anm1 = 0.1e+01;
for 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.1e+01;
end; i = fix(25+1);
while (1);
if( integr~=1 )
if( integr~=3 )
%
%           COMPUTE RG USING A FORWARD RECURRENCE RELATION.
%
rg(1) = -ri(1)./alfp1;
rg(2) = -(ralf+ralf)./(alfp2.*alfp2) - rg(1);
an = 0.2e+01;
anm1 = 0.1e+01;
im1 = 2;
for i = 3 : 25;
rg(i) = -(an.*(an-alfp2).*rg(im1)-an.*ri(im1)+anm1.*ri(i))./(anm1.*(an+alfp1));
anm1 = an;
an = an + 0.1e+01;
im1 = fix(i);
end; i = fix(25+1);
if( integr==2 )
break;
end;
end;
%
%           COMPUTE RH USING A FORWARD RECURRENCE RELATION.
%
rh(1) = -rj(1)./betp1;
rh(2) = -(rbet+rbet)./(betp2.*betp2) - rh(1);
an = 0.2e+01;
anm1 = 0.1e+01;
im1 = 2;
for i = 3 : 25;
rh(i) = -(an.*(an-betp2).*rh(im1)-an.*rj(im1)+anm1.*rj(i))./(anm1.*(an+betp1));
anm1 = an;
an = an + 0.1e+01;
im1 = fix(i);
end; i = fix(25+1);
for i = 2 : 2: 25 ;
rh(i) = -rh(i);
end; i = fix(25 +1);
end;
break;
end;
for i = 2 : 2: 25 ;
rj(i) = -rj(i);
end; i = fix(25 +1);
end %subroutine qmomo
%DECK QNC79