C  Copyright (c) 2003-2010 University of Florida
C
C  This program is free software; you can redistribute it and/or modify
C  it under the terms of the GNU General Public License as published by
C  the Free Software Foundation; either version 2 of the License, or
C  (at your option) any later version.

C  This program is distributed in the hope that it will be useful,
C  but WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C  GNU General Public License for more details.

C  The GNU General Public License is included in this distribution
C  in the file COPYRIGHT.
         SUBROUTINE  OED__RYS_3_ROOTS_WEIGHTS
     +
     +                    ( NT,NTGQP,
     +                      TVAL,
     +
     +                             RTS,
     +                             WTS )
     +
C------------------------------------------------------------------------
C  OPERATION   : OED__RYS_3_ROOTS_WEIGHTS
C  MODULE      : ONE ELECTRON INTEGRALS DIRECT
C  MODULE-ID   : OED
C  SUBROUTINES : none
C  DESCRIPTION : This operation returns Rys polynomial roots and weights
C                in case the number of roots and weights required
C                is = 3. All T's are treated at once so the complete
C                set of roots and weights is returned.
C
C                For the moment taken essentially unchanged from the
C                GAMESS package (routine RTS123, but removing their
C                'spaghetti' code from the 70's of unreadable
C                internested IFs and GOTOs!).
C
C                One interesting aspect of the GAMESS routines is that
C                their code returns scaled roots, i.e. their roots
C                do not ly between the range 0 and 1. To get to the
C                proper roots as needed for our package, we simply
C                set:
C
C                   root (our) = root (gamess) / (1 + root (games))
C
C
C                  Input:
C
C                    NT           =  # of T-exponents
C                    NTGQP        =  # of roots times # of T-exponents
C                                    (= 3 * NT)
C                    TVAL         =  the set of NT T-exponents defining
C                                    the Rys weight functions
C
C                  Output:
C
C                    RTS          =  all NTGQP quadrature roots
C                    WTS          =  all NTGQP quadrature weights
C
C
C  AUTHOR      : Norbert Flocke
C------------------------------------------------------------------------
C
C
C             ...include files and declare variables.
C
C
         IMPLICIT   NONE

         INTEGER    M,N
         INTEGER    NT,NTGQP
         INTEGER    TCASE

         INTEGER    JUMP3 (1:48)

         DOUBLE PRECISION   A1,A2
         DOUBLE PRECISION   E,T,X,Y,Z
         DOUBLE PRECISION   F1,F2
         DOUBLE PRECISION   PI4
         DOUBLE PRECISION   R1,R2,R3
         DOUBLE PRECISION   R13,R23,R33
         DOUBLE PRECISION   T1,T2,T3
         DOUBLE PRECISION   T3MAX
         DOUBLE PRECISION   W1,W2,W3
         DOUBLE PRECISION   W23,W33
         DOUBLE PRECISION   THIRD,HALF,ONE,TWO,FOUR,HALFP7,HALFP12

         DOUBLE PRECISION   TVAL (1:NT)
         DOUBLE PRECISION   RTS  (1:NTGQP)
         DOUBLE PRECISION   WTS  (1:NTGQP)

         PARAMETER  (R13 = 1.90163509193487D-01)
         PARAMETER  (R23 = 1.78449274854325D+00)
         PARAMETER  (R33 = 5.52534374226326D+00)

         PARAMETER  (W23 = 1.77231492083829D-01)
         PARAMETER  (W33 = 5.11156880411248D-03)

         PARAMETER  (PI4 = 7.85398163397448D-01)

         PARAMETER  (T3MAX = 48.D0)

         PARAMETER  (THIRD   = 0.333333333333333D0)
         PARAMETER  (HALF    = 0.5D0)
         PARAMETER  (ONE     = 1.D0)
         PARAMETER  (TWO     = 2.D0)
         PARAMETER  (FOUR    = 4.D0)
         PARAMETER  (HALFP7  = 7.5D0)
         PARAMETER  (HALFP12 = 12.5D0)

         DATA  JUMP3  /1,2,2,3,3,4,4,4,4,4,
     +                 5,5,5,5,5,6,6,6,6,6,
     +                 7,7,7,7,7,7,7,7,7,7,
     +                 7,7,7,8,8,8,8,8,8,8,
     +                 8,8,8,8,8,8,8,9/
C
C
C------------------------------------------------------------------------
C
C
C                 ********************************
C             ... *  # of roots and weights = 3  *
C                 ********************************
C
C
      M = 1

      DO 300 N = 1,NT

         T = TVAL (N)

         IF (T. LE. 3.0D-07) THEN
C
C
C             ...T-range: T essentially 0
C
C
             R1 = 6.03769246832797D-02 - 9.28875764357368D-03 * T
             R2 = 7.76823355931043D-01 - 1.19511285527878D-01 * T
             R3 = 6.66279971938567D+00 - 1.02504611068957D+00 * T

             WTS (M)   = 4.67913934572691D-01 - 5.64876917232519D-02 * T
             WTS (M+1) = 3.60761573048137D-01 - 1.49077186455208D-01 * T
             WTS (M+2) = 1.71324492379169D-01 - 1.27768455150979D-01 * T

             RTS (M)   = R1 / (ONE + R1)
             RTS (M+1) = R2 / (ONE + R2)
             RTS (M+2) = R3 / (ONE + R3)

             M = M + 3

             GOTO 300
         END IF

         TCASE = INT ( MIN (T+ONE,T3MAX))

         GOTO (3100,3200,3300,3400,
     +         3500,3600,3700,3800,3900) JUMP3 (TCASE)
C
C
C             ...T-range: 0 < T < 1
C
C
 3100    E = DEXP (-T)

         F2 = ((((((((-7.60911486098850D-08  * T
     +                +1.09552870123182D-06) * T
     +                -1.03463270693454D-05) * T
     +                +8.16324851790106D-05) * T
     +                -5.55526624875562D-04) * T
     +                +3.20512054753924D-03) * T
     +                -1.51515139838540D-02) * T
     +                +5.55555554649585D-02) * T
     +                -1.42857142854412D-01) * T
     +                +1.99999999999986D-01

         F1 = ((T + T) * F2 + E) * THIRD

         W1 = (T + T) * F1 + E

         R1 = ((((((-5.10186691538870D-10  * T
     +              +2.40134415703450D-08) * T
     +              -5.01081057744427D-07) * T
     +              +7.58291285499256D-06) * T
     +              -9.55085533670919D-05) * T
     +              +1.02893039315878D-03) * T
     +              -9.28875764374337D-03) * T
     +              +6.03769246832810D-02

         R2 = ((((((-1.29646524960555D-08  * T
     +              +7.74602292865683D-08) * T
     +              +1.56022811158727D-06) * T
     +              -1.58051990661661D-05) * T
     +              -3.30447806384059D-04) * T
     +              +9.74266885190267D-03) * T
     +              -1.19511285526388D-01) * T
     +              +7.76823355931033D-01

         R3 = ((((((-9.28536484109606D-09  * T
     +              -3.02786290067014D-07) * T
     +              -2.50734477064200D-06) * T
     +              -7.32728109752881D-06) * T
     +              +2.44217481700129D-04) * T
     +              +4.94758452357327D-02) * T
     +              -1.02504611065774D+00) * T
     +              +6.66279971938553D+00

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 1 =< T < 3
C
C
 3200    E = DEXP (-T)
         X = T - TWO

         F2 = ((((((((((-1.48044231072140D-10  * X
     +                  +1.78157031325097D-09) * X
     +                  -1.92514145088973D-08) * X
     +                  +1.92804632038796D-07) * X
     +                  -1.73806555021045D-06) * X
     +                  +1.39195169625425D-05) * X
     +                  -9.74574633246452D-05) * X
     +                  +5.83701488646511D-04) * X
     +                  -2.89955494844975D-03) * X
     +                  +1.13847001113810D-02) * X
     +                  -3.23446977320647D-02) * X
     +                  +5.29428148329709D-02

         F1 = ((T + T) * F2 + E) * THIRD

         W1 = (T + T) * F1 + E

         R1 = (((((((( 1.44687969563318D-12  * X
     +                +4.85300143926755D-12) * X
     +                -6.55098264095516D-10) * X
     +                +1.56592951656828D-08) * X
     +                -2.60122498274734D-07) * X
     +                +3.86118485517386D-06) * X
     +                -5.13430986707889D-05) * X
     +                +6.03194524398109D-04) * X
     +                -6.11219349825090D-03) * X
     +                +4.52578254679079D-02

         R2 = ((((((( 6.95964248788138D-10  * X
     +               -5.35281831445517D-09) * X
     +               -6.745205954533D-08)   * X
     +               +1.502366784525D-06)   * X
     +               +9.923326947376D-07)   * X
     +               -3.89147469249594D-04) * X
     +               +7.51549330892401D-03) * X
     +               -8.48778120363400D-02) * X
     +               +5.73928229597613D-01

         R3 = ((((((((-2.81496588401439D-10  * X
     +                +3.61058041895031D-09) * X
     +                +4.53631789436255D-08) * X
     +                -1.40971837780847D-07) * X
     +                -6.05865557561067D-06) * X
     +                -5.15964042227127D-05) * X
     +                +3.34761560498171D-05) * X
     +                +5.04871005319119D-02) * X
     +                -8.24708946991557D-01) * X
     +                +4.81234667357205D+00

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 3 =< T < 5
C
C
 3300    E = DEXP (-T)
         X = T - FOUR

         F2 = ((((((((((-2.36788772599074D-11  * X
     +                  +2.89147476459092D-10) * X
     +                  -3.18111322308846D-09) * X
     +                  +3.25336816562485D-08) * X
     +                  -3.00873821471489D-07) * X
     +                  +2.48749160874431D-06) * X
     +                  -1.81353179793672D-05) * X
     +                  +1.14504948737066D-04) * X
     +                  -6.10614987696677D-04) * X
     +                  +2.64584212770942D-03) * X
     +                  -8.66415899015349D-03) * X
     +                  +1.75257821619922D-02

         F1 = ((T + T) * F2 + E) * THIRD

         W1 = (T + T) * F1 + E

         R1 = ((((((( 1.44265709189601D-11  * X
     +               -4.66622033006074D-10) * X
     +               +7.649155832025D-09)   * X
     +               -1.229940017368D-07)   * X
     +               +2.026002142457D-06)   * X
     +               -2.87048671521677D-05) * X
     +               +3.70326938096287D-04) * X
     +               -4.21006346373634D-03) * X
     +               +3.50898470729044D-02

         R2 = ((((((((-2.65526039155651D-11  * X
     +                +1.97549041402552D-10) * X
     +                +2.15971131403034D-09) * X
     +                -7.95045680685193D-08) * X
     +                +5.15021914287057D-07) * X
     +                +1.11788717230514D-05) * X
     +                -3.33739312603632D-04) * X
     +                +5.30601428208358D-03) * X
     +                -5.93483267268959D-02) * X
     +                +4.31180523260239D-01

         R3 = ((((((((-3.92833750584041D-10  * X
     +                -4.16423229782280D-09) * X
     +                +4.42413039572867D-08) * X
     +                +6.40574545989551D-07) * X
     +                -3.05512456576552D-06) * X
     +                -1.05296443527943D-04) * X
     +                -6.14120969315617D-04) * X
     +                +4.89665802767005D-02) * X
     +                -6.24498381002855D-01) * X
     +                +3.36412312243724D+00

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 5 =< T < 10
C
C
 3400    E = DEXP (-T)
         X = ONE / T
         Y = T - HALFP7
         Z = HALF * X

         W1 = (((((( 4.6897511375022D-01  * X
     +              -6.9955602298985D-01) * X
     +              +5.3689283271887D-01) * X
     +              -3.2883030418398D-01) * X
     +              +2.4645596956002D-01) * X
     +              -4.9984072848436D-01) * X
     +              -3.1501078774085D-06) * E + DSQRT (PI4*X)

         F1 = (W1 - E) * Z
         F2 = (F1 + F1 + F1 - E) * Z

         R1 = ((((((((((( 5.74429401360115D-16  * Y
     +                   +7.11884203790984D-16) * Y
     +                   -6.736701449826D-14)   * Y
     +                   -6.264613873998D-13)   * Y
     +                   +1.315418927040D-11)   * Y
     +                   -4.23879635610964D-11) * Y
     +                   +1.39032379769474D-09) * Y
     +                   -4.65449552856856D-08) * Y
     +                   +7.34609900170759D-07) * Y
     +                   -1.08656008854077D-05) * Y
     +                   +1.77930381549953D-04) * Y
     +                   -2.39864911618015D-03) * Y
     +                   +2.39112249488821D-02

         R2 = ((((((((((( 1.13464096209120D-14  * Y
     +                   +6.99375313934242D-15) * Y
     +                   -8.595618132088D-13)   * Y
     +                   -5.293620408757D-12)   * Y
     +                   -2.492175211635D-11)   * Y
     +                   +2.73681574882729D-09) * Y
     +                   -1.06656985608482D-08) * Y
     +                   -4.40252529648056D-07) * Y
     +                   +9.68100917793911D-06) * Y
     +                   -1.68211091755327D-04) * Y
     +                   +2.69443611274173D-03) * Y
     +                   -3.23845035189063D-02) * Y
     +                   +2.75969447451882D-01

         R3 = (((((((((((( 6.66339416996191D-15  * Y
     +                    +1.84955640200794D-13) * Y
     +                    -1.985141104444D-12)   * Y
     +                    -2.309293727603D-11)   * Y
     +                    +3.917984522103D-10)   * Y
     +                    +1.663165279876D-09)   * Y
     +                    -6.205591993923D-08)   * Y
     +                    +8.769581622041D-09)   * Y
     +                    +8.97224398620038D-06) * Y
     +                    -3.14232666170796D-05) * Y
     +                    -1.83917335649633D-03) * Y
     +                    +3.51246831672571D-02) * Y
     +                    -3.22335051270860D-01) * Y
     +                    +1.73582831755430D+00

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 10 =< T < 15
C
C
 3500    E = DEXP (-T)
         X = ONE / T
         Y = T - HALFP12
         Z = HALF * X

         W1 = (((-1.8784686463512D-01  * X
     +           +2.2991849164985D-01) * X
     +           -4.9893752514047D-01) * X
     +           -2.1916512131607D-05) * E + DSQRT (PI4*X)

         F1 = (W1 - E) * Z
         F2 = (F1 + F1 + F1 - E) * Z

         R1 = ((((((((((( 4.42133001283090D-16  * Y
     +                   -2.77189767070441D-15) * Y
     +                   -4.084026087887D-14)   * Y
     +                   +5.379885121517D-13)   * Y
     +                   +1.882093066702D-12)   * Y
     +                   -8.67286219861085D-11) * Y
     +                   +7.11372337079797D-10) * Y
     +                   -3.55578027040563D-09) * Y
     +                   +1.29454702851936D-07) * Y
     +                   -4.14222202791434D-06) * Y
     +                   +8.04427643593792D-05) * Y
     +                   -1.18587782909876D-03) * Y
     +                   +1.53435577063174D-02

         R2 = ((((((((((( 6.85146742119357D-15  * Y
     +                   -1.08257654410279D-14) * Y
     +                   -8.579165965128D-13)   * Y
     +                   +6.642452485783D-12)   * Y
     +                   +4.798806828724D-11)   * Y
     +                   -1.13413908163831D-09) * Y
     +                   +7.08558457182751D-09) * Y
     +                   -5.59678576054633D-08) * Y
     +                   +2.51020389884249D-06) * Y
     +                   -6.63678914608681D-05) * Y
     +                   +1.11888323089714D-03) * Y
     +                   -1.45361636398178D-02) * Y
     +                   +1.65077877454402D-01

         R3 = (((((((((((( 3.20622388697743D-15  * Y
     +                    -2.73458804864628D-14) * Y
     +                    -3.157134329361D-13)   * Y
     +                    +8.654129268056D-12)   * Y
     +                    -5.625235879301D-11)   * Y
     +                    -7.718080513708D-10)   * Y
     +                    +2.064664199164D-08)   * Y
     +                    -1.567725007761D-07)   * Y
     +                    -1.57938204115055D-06) * Y
     +                    +6.27436306915967D-05) * Y
     +                    -1.01308723606946D-03) * Y
     +                    +1.13901881430697D-02) * Y
     +                    -1.01449652899450D-01) * Y
     +                    +7.77203937334739D-01

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 15 =< T < 20
C
C
 3600    E = DEXP (-T)
         X = ONE / T
         Z = HALF * X

         W1 = (( 1.9623264149430D-01  * X
     +          -4.9695241464490D-01) * X
     +          -6.0156581186481D-05) * E + DSQRT (PI4*X)

         F1 = (W1 - E) * Z
         F2 = (F1 + F1 + F1 - E) * Z

         R1 = ((((((-2.43270989903742D-06  * T
     +              +3.57901398988359D-04) * T
     +              -2.34112415981143D-02) * T
     +              +7.81425144913975D-01) * T
     +              -1.73209218219175D+01) * T
     +              +2.43517435690398D+02) * T
     +              + (-1.97611541576986D+04  * X
     +                 +9.82441363463929D+03) * X
     +              -2.07970687843258D+03) * E + R13 / (T - R13)

         R2 = (((((-2.62627010965435D-04  * T
     +             +3.49187925428138D-02) * T
     +             -3.09337618731880D+00) * T
     +             +1.07037141010778D+02) * T
     +             -2.36659637247087D+03) * T
     +             + ((-2.91669113681020D+06  * X
     +                 +1.41129505262758D+06) * X
     +                 -2.91532335433779D+05) * X
     +             +3.35202872835409D+04) * E + R23 / (T - R23)

         R3 = ((((( 9.31856404738601D-05  * T
     +             -2.87029400759565D-02) * T
     +             -7.83503697918455D-01) * T
     +             -1.84338896480695D+01) * T
     +             +4.04996712650414D+02) * T
     +             + (-1.89829509315154D+05  * X
     +                +5.11498390849158D+04) * X
     +             -6.88145821789955D+03) * E + R33 / (T - R33)

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 20 =< T < 33
C
C
 3700    E = DEXP (-T)
         X = ONE / T
         Z = HALF * X

         W1 = (( 1.9623264149430D-01  * X
     +          -4.9695241464490D-01) * X
     +          -6.0156581186481D-05) * E + DSQRT (PI4*X)

         F1 = (W1 - E) * Z
         F2 = (F1 + F1 + F1 - E) * Z

         R1 = ((((-4.97561537069643D-04  * T
     +            -5.00929599665316D-02) * T
     +            +1.31099142238996D+00) * T
     +            -1.88336409225481D+01) * T
     +            -6.60344754467191D+02  * X
     +            +1.64931462413877D+02) * E + R13 / (T - R13)

         R2 = ((((-4.48218898474906D-03  * T
     +            -5.17373211334924D-01) * T
     +            +1.13691058739678D+01) * T
     +            -1.65426392885291D+02) * T
     +            -6.30909125686731D+03  * X
     +            +1.52231757709236D+03) * E + R23 / (T - R23)

         R3 = ((((-1.38368602394293D-02  * T
     +            -1.77293428863008D+00) * T
     +            +1.73639054044562D+01) * T
     +            -3.57615122086961D+02) * T
     +            -1.45734701095912D+04  * X
     +            +2.69831813951849D+03) * E + R33 / (T - R33)

         T1 = R1 / (R1 + ONE)
         T2 = R2 / (R2 + ONE)
         T3 = R3 / (R3 + ONE)

         A1 = F1 - T1 * W1
         A2 = F2 - T1 * F1

         W2 = (T3 * A1 - A2) / ((T3 - T2) * (T2 - T1))
         W3 = (A2 - T2 * A1) / ((T3 - T2) * (T3 - T1))

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = T1
         RTS (M+1) = T2
         RTS (M+2) = T3

         M = M + 3

         GOTO 300
C
C
C             ...T-range: 33 =< T < 47
C
C
 3800    E = DEXP (-T)

         W1 = DSQRT (PI4/T)

         W2 = (( 6.15072615497811D+01  * T
     +          -2.91980647450269D+03) * T
     +          +3.80794303087338D+04) * E + W23 * W1

         W3 = ((( 1.52258947224714D-01  * T
     +           -8.30661900042651D+00) * T
     +           +1.92977367967984D+02) * T
     +           -1.67787926005344D+03) * E + W33 * W1

         R1 = ((-7.39058467995275D+00  * T
     +          +3.21318352526305D+02) * T
     +          -3.99433696473658D+03) * E + R13 / (T - R13)

         R2 = ((-7.38726243906513D+01  * T
     +          +3.13569966333873D+03) * T
     +          -3.86862867311321D+04) * E + R23 / (T - R23)

         R3 = ((-2.63750565461336D+02  * T
     +          +1.04412168692352D+04) * T
     +          -1.28094577915394D+05) * E + R33 / (T - R33)

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = R1 / (ONE + R1)
         RTS (M+1) = R2 / (ONE + R2)
         RTS (M+2) = R3 / (ONE + R3)

         M = M + 3

         GOTO 300
C
C
C             ...T-range: T >= 47
C
C
 3900    W1 = DSQRT (PI4/T)
         W2 = W23 * W1
         W3 = W33 * W1

C         R1 = R13 / (T - R13)
C         R2 = R23 / (T - R23)
C         R3 = R33 / (T - R33)
C         RTS (M)   = R1 / (ONE + R1)
C         RTS (M+1) = R2 / (ONE + R2)
C         RTS (M+2) = R3 / (ONE + R3)

         WTS (M)   = W1 - W2 - W3
         WTS (M+1) = W2
         WTS (M+2) = W3

         RTS (M)   = R13 / T
         RTS (M+1) = R23 / T
         RTS (M+2) = R33 / T

         M = M + 3

  300 CONTINUE
C
C
C             ...ready!
C
C
      RETURN
      END
