C 17 JAN 06 - MWS - REPLACE ROOT6 WITH DISCRETIZED STIELTJES CODE
C 19 SEP 05 - GDF - DIMENSION ROOT6 FAMILY FOR 13 ROOTS
C  7 MAR 05 - JMS - ROOTSA: EXTEND LARGE X CODE TO NROOTS=13
C 13 FEB 05 - MWS - PAD COMMON BLOCK ROOT
C 10 NOV 04 - JMS - USE ROOTSA FOR ASYMPTOTICALLY LARGE X VALUES
C 20 NOV 97 - MWS - RYSFUN: TEST ON ASYMPTOTIC NUMERICAL ACCURACY
C 15 NOV 94 - DK  - RYSSMT: AVOID SQRT OF A NEGATIVE NUMBER IN FAC
C  7 JAN 92 - TLW - MAKE WRITE PARALLEL; ADD COMMON PAR
C 28 SEP 90 - MWS - RYSFUN: DIMINISH E USED FOR LARGE X ASYMPTOTE
C                   RYSSMT: AVOID DIVISION BY ZERO WHEN FAC UNDERFLOWS
C 21 SEP 89 - STE - RYSNOD: CALL EXIT CHANGED TO CALL ABRT
C 10 APR 88 - MWS - ADD ROOT6,RYSFUN,RYSSMT,RYSNOD
C 14 FEB 88 - MWS - INCREASE RYS ROOTS TO 9
C 13 OCT 86 - JAB - REMOVE BLANK FROM CONSTANT IN ROOT5
C 21 OCT 85 - STE - USE GENERIC FUNCTIONS
C  2 APR 84 - STE - ISOLATE ROOT ROUTINES
C
C*MODULE RYSPOL  *DECK RT123
      SUBROUTINE RT123
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C             *****   VERSION FEBRUARY 13,1975   *****
C
      COMMON /ROOT  / X,U(13),W(13),NROOTS
C
      EQUIVALENCE (U(1),RT1),(U(2),RT2),(U(3),RT3),(U(4),RT4),(U(5),RT5)
      EQUIVALENCE (W(1),WW1),(W(2),WW2),(W(3),WW3),(W(4),WW4),(W(5),WW5)
C
      DATA R12,PIE4/2.75255128608411D-01, 7.85398163397448D-01/
      DATA R22,W22/ 2.72474487139158D+00, 9.17517095361369D-02/
      DATA R13/     1.90163509193487D-01/
      DATA R23,W23/ 1.78449274854325D+00, 1.77231492083829D-01/
      DATA R33,W33/ 5.52534374226326D+00, 5.11156880411248D-03/
C
C     X LARGE (ASYMPTOTIC SOLUTION)              NROOTS=1,2,3
C
CA    IF((NROOTS.EQ.1 .AND. X.GT.33.0D+00) .OR.
CA   *   (NROOTS.EQ.2 .AND. X.GT.40.0D+00) .OR.
CA   *   (NROOTS.EQ.3 .AND. X.GT.47.0D+00)) THEN
CA       CALL ROOTSA
CA       RETURN
CA    ENDIF
C
      IF (X .GT. 5.0D+00) GO TO 400
      IF (X .GT. 1.0D+00) GO TO 280
      IF (X .GT. 3.0D-07) GO TO 180
C     X IS APPROXIMATELY ZERO.         NROOTS=1,2, OR 3
      IF (NROOTS-2) 120,140,160
  120 RT1 = 0.5D+00 -X/5.0D+00
      WW1 = 1.0D+00 -X/3.0D+00
      RETURN
  140 RT1 = 1.30693606237085D-01 -2.90430236082028D-02 *X
      RT2 = 2.86930639376291D+00 -6.37623643058102D-01 *X
      WW1 = 6.52145154862545D-01 -1.22713621927067D-01 *X
      WW2 = 3.47854845137453D-01 -2.10619711404725D-01 *X
      RETURN
  160 RT1 = 6.03769246832797D-02 -9.28875764357368D-03 *X
      RT2 = 7.76823355931043D-01 -1.19511285527878D-01 *X
      RT3 = 6.66279971938567D+00 -1.02504611068957D+00 *X
      WW1 = 4.67913934572691D-01 -5.64876917232519D-02 *X
      WW2 = 3.60761573048137D-01 -1.49077186455208D-01 *X
      WW3 = 1.71324492379169D-01 -1.27768455150979D-01 *X
      RETURN
C     X = 0.0 TO 1.0                   NROOTS=1,2, OR 3
  180 IF (NROOTS .EQ. 3) GO TO 220
      F1 = ((((((((-8.36313918003957D-08*X+1.21222603512827D-06 )*X-
     +     1.15662609053481D-05 )*X+9.25197374512647D-05 )*X-
     +     6.40994113129432D-04 )*X+3.78787044215009D-03 )*X-
     +     1.85185172458485D-02 )*X+7.14285713298222D-02 )*X-
     +     1.99999999997023D-01 )*X+3.33333333333318D-01
      WW1 = (X+X)*F1+EXP(-X)
      IF (NROOTS .EQ. 2) GO TO 200
      RT1 = F1/(WW1-F1)
      RETURN
  200 RT1 = (((((((-2.35234358048491D-09*X+2.49173650389842D-08)*X-
     +     4.558315364581D-08)*X-2.447252174587D-06)*X+
     +     4.743292959463D-05)*X-5.33184749432408D-04 )*X+
     +     4.44654947116579D-03 )*X-2.90430236084697D-02 )*X+
     +     1.30693606237085D-01
      RT2 = (((((((-2.47404902329170D-08*X+2.36809910635906D-07)*X+
     +     1.835367736310D-06)*X-2.066168802076D-05)*X-
     +     1.345693393936D-04)*X-5.88154362858038D-05 )*X+
     +     5.32735082098139D-02 )*X-6.37623643056745D-01 )*X+
     +     2.86930639376289D+00
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  220 RT1 = ((((((-5.10186691538870D-10*X+2.40134415703450D-08)*X-
     +     5.01081057744427D-07 )*X+7.58291285499256D-06 )*X-
     +     9.55085533670919D-05 )*X+1.02893039315878D-03 )*X-
     +     9.28875764374337D-03 )*X+6.03769246832810D-02
      RT2 = ((((((-1.29646524960555D-08*X+7.74602292865683D-08)*X+
     +     1.56022811158727D-06 )*X-1.58051990661661D-05 )*X-
     +     3.30447806384059D-04 )*X+9.74266885190267D-03 )*X-
     +     1.19511285526388D-01 )*X+7.76823355931033D-01
      RT3 = ((((((-9.28536484109606D-09*X-3.02786290067014D-07)*X-
     +     2.50734477064200D-06 )*X-7.32728109752881D-06 )*X+
     +     2.44217481700129D-04 )*X+4.94758452357327D-02 )*X-
     +     1.02504611065774D+00 )*X+6.66279971938553D+00
      F2 = ((((((((-7.60911486098850D-08*X+1.09552870123182D-06 )*X-
     +     1.03463270693454D-05 )*X+8.16324851790106D-05 )*X-
     +     5.55526624875562D-04 )*X+3.20512054753924D-03 )*X-
     +     1.51515139838540D-02 )*X+5.55555554649585D-02 )*X-
     +     1.42857142854412D-01 )*X+1.99999999999986D-01
  240 E = EXP(-X)
      F1 = ((X+X)*F2+E)/3.0D+00
      WW1 = (X+X)*F1+E
  260 T1 = RT1/(RT1+1.0D+00)
      T2 = RT2/(RT2+1.0D+00)
      T3 = RT3/(RT3+1.0D+00)
      A2 = F2-T1*F1
      A1 = F1-T1*WW1
      WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1))
      WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1))
      WW1 = WW1-WW2-WW3
      RETURN
  280 IF (X .GT. 3.0D+00) GO TO 340
C     X = 1.0 TO 3.0                   NROOTS=1,2, OR 3
      Y = X-2.0D+00
      IF (NROOTS .EQ. 3) GO TO 320
      F1 = ((((((((((-1.61702782425558D-10*Y+1.96215250865776D-09 )*Y-
     +     2.14234468198419D-08 )*Y+2.17216556336318D-07 )*Y-
     +     1.98850171329371D-06 )*Y+1.62429321438911D-05 )*Y-
     +     1.16740298039895D-04 )*Y+7.24888732052332D-04 )*Y-
     +     3.79490003707156D-03 )*Y+1.61723488664661D-02 )*Y-
     +     5.29428148329736D-02 )*Y+1.15702180856167D-01
      WW1 = (X+X)*F1+EXP(-X)
      IF (NROOTS .EQ. 2) GO TO 300
      RT1 = F1/(WW1-F1)
      RETURN
  300 RT1 = (((((((((-6.36859636616415D-12*Y+8.47417064776270D-11)*Y-
     +     5.152207846962D-10)*Y-3.846389873308D-10)*Y+
     +     8.472253388380D-08)*Y-1.85306035634293D-06 )*Y+
     +     2.47191693238413D-05 )*Y-2.49018321709815D-04 )*Y+
     +     2.19173220020161D-03 )*Y-1.63329339286794D-02 )*Y+
     +     8.68085688285261D-02
      RT2 = ((((((((( 1.45331350488343D-10*Y+2.07111465297976D-09)*Y-
     +     1.878920917404D-08)*Y-1.725838516261D-07)*Y+
     +     2.247389642339D-06)*Y+9.76783813082564D-06 )*Y-
     +     1.93160765581969D-04 )*Y-1.58064140671893D-03 )*Y+
     +     4.85928174507904D-02 )*Y-4.30761584997596D-01 )*Y+
     +     1.80400974537950D+00
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  320 RT1 = (((((((( 1.44687969563318D-12*Y+4.85300143926755D-12)*Y-
     +     6.55098264095516D-10 )*Y+1.56592951656828D-08 )*Y-
     +     2.60122498274734D-07 )*Y+3.86118485517386D-06 )*Y-
     +     5.13430986707889D-05 )*Y+6.03194524398109D-04 )*Y-
     +     6.11219349825090D-03 )*Y+4.52578254679079D-02
      RT2 = ((((((( 6.95964248788138D-10*Y-5.35281831445517D-09)*Y-
     +     6.745205954533D-08)*Y+1.502366784525D-06)*Y+
     +     9.923326947376D-07)*Y-3.89147469249594D-04 )*Y+
     +     7.51549330892401D-03 )*Y-8.48778120363400D-02 )*Y+
     +     5.73928229597613D-01
      RT3 = ((((((((-2.81496588401439D-10*Y+3.61058041895031D-09)*Y+
     +     4.53631789436255D-08 )*Y-1.40971837780847D-07 )*Y-
     +     6.05865557561067D-06 )*Y-5.15964042227127D-05 )*Y+
     +     3.34761560498171D-05 )*Y+5.04871005319119D-02 )*Y-
     +     8.24708946991557D-01 )*Y+4.81234667357205D+00
      F2 = ((((((((((-1.48044231072140D-10*Y+1.78157031325097D-09 )*Y-
     +     1.92514145088973D-08 )*Y+1.92804632038796D-07 )*Y-
     +     1.73806555021045D-06 )*Y+1.39195169625425D-05 )*Y-
     +     9.74574633246452D-05 )*Y+5.83701488646511D-04 )*Y-
     +     2.89955494844975D-03 )*Y+1.13847001113810D-02 )*Y-
     +     3.23446977320647D-02 )*Y+5.29428148329709D-02
      GO TO 240
C     X = 3.0 TO 5.0                   NROOTS =1,2, OR 3
  340 Y = X-4.0D+00
      IF (NROOTS .EQ. 3) GO TO 380
      F1 = ((((((((((-2.62453564772299D-11*Y+3.24031041623823D-10 )*Y-
     +     3.614965656163D-09)*Y+3.760256799971D-08)*Y-
     +     3.553558319675D-07)*Y+3.022556449731D-06)*Y-
     +     2.290098979647D-05)*Y+1.526537461148D-04)*Y-
     +     8.81947375894379D-04 )*Y+4.33207949514611D-03 )*Y-
     +     1.75257821619926D-02 )*Y+5.28406320615584D-02
      WW1 = (X+X)*F1+EXP(-X)
      IF (NROOTS .EQ. 2) GO TO 360
      RT1 = F1/(WW1-F1)
      RETURN
  360 RT1 = ((((((((-4.11560117487296D-12*Y+7.10910223886747D-11)*Y-
     +     1.73508862390291D-09 )*Y+5.93066856324744D-08 )*Y-
     +     9.76085576741771D-07 )*Y+1.08484384385679D-05 )*Y-
     +     1.12608004981982D-04 )*Y+1.16210907653515D-03 )*Y-
     +     9.89572595720351D-03 )*Y+6.12589701086408D-02
      RT2 = (((((((((-1.80555625241001D-10*Y+5.44072475994123D-10)*Y+
     +     1.603498045240D-08)*Y-1.497986283037D-07)*Y-
     +     7.017002532106D-07)*Y+1.85882653064034D-05 )*Y-
     +     2.04685420150802D-05 )*Y-2.49327728643089D-03 )*Y+
     +     3.56550690684281D-02 )*Y-2.60417417692375D-01 )*Y+
     +     1.12155283108289D+00
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  380 RT1 = ((((((( 1.44265709189601D-11*Y-4.66622033006074D-10)*Y+
     +     7.649155832025D-09)*Y-1.229940017368D-07)*Y+
     +     2.026002142457D-06)*Y-2.87048671521677D-05 )*Y+
     +     3.70326938096287D-04 )*Y-4.21006346373634D-03 )*Y+
     +     3.50898470729044D-02
      RT2 = ((((((((-2.65526039155651D-11*Y+1.97549041402552D-10)*Y+
     +     2.15971131403034D-09 )*Y-7.95045680685193D-08 )*Y+
     +     5.15021914287057D-07 )*Y+1.11788717230514D-05 )*Y-
     +     3.33739312603632D-04 )*Y+5.30601428208358D-03 )*Y-
     +     5.93483267268959D-02 )*Y+4.31180523260239D-01
      RT3 = ((((((((-3.92833750584041D-10*Y-4.16423229782280D-09)*Y+
     +     4.42413039572867D-08 )*Y+6.40574545989551D-07 )*Y-
     +     3.05512456576552D-06 )*Y-1.05296443527943D-04 )*Y-
     +     6.14120969315617D-04 )*Y+4.89665802767005D-02 )*Y-
     +     6.24498381002855D-01 )*Y+3.36412312243724D+00
      F2 = ((((((((((-2.36788772599074D-11*Y+2.89147476459092D-10 )*Y-
     +     3.18111322308846D-09 )*Y+3.25336816562485D-08 )*Y-
     +     3.00873821471489D-07 )*Y+2.48749160874431D-06 )*Y-
     +     1.81353179793672D-05 )*Y+1.14504948737066D-04 )*Y-
     +     6.10614987696677D-04 )*Y+2.64584212770942D-03 )*Y-
     +     8.66415899015349D-03 )*Y+1.75257821619922D-02
      GO TO 240
  400 IF (X .GT. 15.0D+00) GO TO 560
      E = EXP(-X)
      IF (X .GT. 10.0D+00) GO TO 480
C     X = 5.0 TO 10.0                  NROOTS =1,2, OR 3
      WW1 = (((((( 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 + SQRT(PIE4/X)
      F1 = (WW1-E)/(X+X)
      IF (NROOTS-2) 420,440,460
  420 RT1 = F1/(WW1-F1)
      RETURN
  440 Y = X-7.5D+00
      RT1 = (((((((((((((-1.43632730148572D-16*Y+2.38198922570405D-16)*
     +     Y+1.358319618800D-14)*Y-7.064522786879D-14)*Y-
     +     7.719300212748D-13)*Y+7.802544789997D-12)*Y+
     +     6.628721099436D-11)*Y-1.775564159743D-09)*Y+
     +     1.713828823990D-08)*Y-1.497500187053D-07)*Y+
     +     2.283485114279D-06)*Y-3.76953869614706D-05 )*Y+
     +     4.74791204651451D-04 )*Y-4.60448960876139D-03 )*Y+
     +     3.72458587837249D-02
      RT2 = (((((((((((( 2.48791622798900D-14*Y-1.36113510175724D-13)*Y-
     +     2.224334349799D-12)*Y+4.190559455515D-11)*Y-
     +     2.222722579924D-10)*Y-2.624183464275D-09)*Y+
     +     6.128153450169D-08)*Y-4.383376014528D-07)*Y-
     +     2.49952200232910D-06 )*Y+1.03236647888320D-04 )*Y-
     +     1.44614664924989D-03 )*Y+1.35094294917224D-02 )*Y-
     +     9.53478510453887D-02 )*Y+5.44765245686790D-01
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  460 F2 = (F1+F1+F1-E)/(X+X)
      Y = X-7.5D+00
      RT1 = ((((((((((( 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
      RT2 = ((((((((((( 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
      RT3 = (((((((((((( 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
      GO TO 260
C     X = 10.0 TO 15.0                 NROOTS=1,2, OR 3
  480 WW1 = (((-1.8784686463512D-01/X+2.2991849164985D-01)/X -
     +     4.9893752514047D-01)/X-2.1916512131607D-05)*E + SQRT(PIE4/X)
      F1 = (WW1-E)/(X+X)
      IF (NROOTS-2) 500,520,540
  500 RT1 = F1/(WW1-F1)
      RETURN
  520 RT1 = ((((-1.01041157064226D-05*X+1.19483054115173D-03)*X -
     +     6.73760231824074D-02)*X+1.25705571069895D+00)*X + (((-
     +     8.57609422987199D+03/X+5.91005939591842D+03)/X -
     +     1.70807677109425D+03)/X+2.64536689959503D+02)/X -
     +     2.38570496490846D+01)*E + R12/(X-R12)
      RT2 = ((( 3.39024225137123D-04*X-9.34976436343509D-02)*X -
     +     4.22216483306320D+00)*X + (((-2.08457050986847D+03/X -
     +     1.04999071905664D+03)/X+3.39891508992661D+02)/X -
     +     1.56184800325063D+02)/X+8.00839033297501D+00)*E + R22/(X-R22)
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  540 F2 = (F1+F1+F1-E)/(X+X)
      Y = X-12.5D+00
      RT1 = ((((((((((( 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
      RT2 = ((((((((((( 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
      RT3 = (((((((((((( 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
      GO TO 260
  560 IF (X .GT. 33.0D+00) GO TO 660
C     X = 15.0 TO 33.0                 NROOTS=1,2, OR 3
      E = EXP(-X)
      WW1 = (( 1.9623264149430D-01/X-4.9695241464490D-01)/X -
     +     6.0156581186481D-05)*E + SQRT(PIE4/X)
      F1 = (WW1-E)/(X+X)
      IF (NROOTS-2) 580,600,620
  580 RT1 = F1/(WW1-F1)
      RETURN
  600 RT1 = ((((-1.14906395546354D-06*X+1.76003409708332D-04)*X -
     +     1.71984023644904D-02)*X-1.37292644149838D-01)*X + (-
     +     4.75742064274859D+01/X+9.21005186542857D+00)/X -
     +     2.31080873898939D-02)*E + R12/(X-R12)
      RT2 = ((( 3.64921633404158D-04*X-9.71850973831558D-02)*X -
     +     4.02886174850252D+00)*X + (-1.35831002139173D+02/X -
     +     8.66891724287962D+01)/X+2.98011277766958D+00)*E + R22/(X-R22)
      WW2 = ((F1-WW1)*RT1+F1)*(1.0D+00+RT2)/(RT2-RT1)
      WW1 = WW1-WW2
      RETURN
  620 F2 = (F1+F1+F1-E)/(X+X)
      IF (X .GT. 20.0D+00) GO TO 640
      RT1 = ((((((-2.43270989903742D-06*X+3.57901398988359D-04)*X -
     +     2.34112415981143D-02)*X+7.81425144913975D-01)*X -
     +     1.73209218219175D+01)*X+2.43517435690398D+02)*X + (-
     +     1.97611541576986D+04/X+9.82441363463929D+03)/X -
     +     2.07970687843258D+03)*E + R13/(X-R13)
      RT2 = (((((-2.62627010965435D-04*X+3.49187925428138D-02)*X -
     +     3.09337618731880D+00)*X+1.07037141010778D+02)*X -
     +     2.36659637247087D+03)*X + ((-2.91669113681020D+06/X +
     +     1.41129505262758D+06)/X-2.91532335433779D+05)/X +
     +     3.35202872835409D+04)*E + R23/(X-R23)
      RT3 = ((((( 9.31856404738601D-05*X-2.87029400759565D-02)*X -
     +     7.83503697918455D-01)*X-1.84338896480695D+01)*X +
     +     4.04996712650414D+02)*X + (-1.89829509315154D+05/X +
     +     5.11498390849158D+04)/X-6.88145821789955D+03)*E + R33/(X-R33)
      GO TO 260
  640 RT1 = ((((-4.97561537069643D-04*X-5.00929599665316D-02)*X +
     +     1.31099142238996D+00)*X-1.88336409225481D+01)*X -
     +     6.60344754467191D+02 /X+1.64931462413877D+02)*E + R13/(X-R13)
      RT2 = ((((-4.48218898474906D-03*X-5.17373211334924D-01)*X +
     +     1.13691058739678D+01)*X-1.65426392885291D+02)*X -
     +     6.30909125686731D+03 /X+1.52231757709236D+03)*E + R23/(X-R23)
      RT3 = ((((-1.38368602394293D-02*X-1.77293428863008D+00)*X +
     +     1.73639054044562D+01)*X-3.57615122086961D+02)*X -
     +     1.45734701095912D+04 /X+2.69831813951849D+03)*E + R33/(X-R33)
      GO TO 260
C     X = 33.0 TO INFINITY             NROOTS=1,2, OR 3
  660 WW1 = SQRT(PIE4/X)
      IF (NROOTS-2) 680,700,720
  680 RT1 = 0.5D+00/(X-0.5D+00)
      RETURN
  700 IF (X .GT. 40.0D+00) GO TO 740
      E = EXP(-X)
      RT1 = (-8.78947307498880D-01*X+1.09243702330261D+01)*E + R12/(X-
     +     R12)
      RT2 = (-9.28903924275977D+00*X+8.10642367843811D+01)*E + R22/(X-
     +     R22)
      WW2 = ( 4.46857389308400D+00*X-7.79250653461045D+01)*E + W22*WW1
      WW1 = WW1-WW2
      RETURN
  720 IF (X .GT. 47.0D+00) GO TO 760
      E = EXP(-X)
      RT1 = ((-7.39058467995275D+00*X+3.21318352526305D+02)*X -
     +     3.99433696473658D+03)*E + R13/(X-R13)
      RT2 = ((-7.38726243906513D+01*X+3.13569966333873D+03)*X -
     +     3.86862867311321D+04)*E + R23/(X-R23)
      RT3 = ((-2.63750565461336D+02*X+1.04412168692352D+04)*X -
     +     1.28094577915394D+05)*E + R33/(X-R33)
      WW3 = ((( 1.52258947224714D-01*X-8.30661900042651D+00)*X +
     +     1.92977367967984D+02)*X-1.67787926005344D+03)*E + W33*WW1
      WW2 = (( 6.15072615497811D+01*X-2.91980647450269D+03)*X +
     +     3.80794303087338D+04)*E + W23*WW1
      WW1 = WW1-WW2-WW3
      RETURN
  740 RT1 = R12/(X-R12)
      RT2 = R22/(X-R22)
      WW2 = W22*WW1
      WW1 = WW1-WW2
      RETURN
  760 RT1 = R13/(X-R13)
      RT2 = R23/(X-R23)
      RT3 = R33/(X-R33)
      WW2 = W23*WW1
      WW3 = W33*WW1
      WW1 = WW1-WW2-WW3
      RETURN
      END
C*MODULE RYSPOL  *DECK ROOT4
      SUBROUTINE ROOT4
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C          *****   VERSION FEBRUARY 16,1975   *****
      COMMON /ROOT  / X,U(13),W(13),NROOTS
C
      EQUIVALENCE (U(1),RT1),(U(2),RT2),(U(3),RT3),(U(4),RT4),(U(5),RT5)
      EQUIVALENCE (W(1),WW1),(W(2),WW2),(W(3),WW3),(W(4),WW4),(W(5),WW5)
C
      DATA R14,PIE4/1.45303521503316D-01, 7.85398163397448D-01/
      DATA R24,W24/ 1.33909728812636D+00, 2.34479815323517D-01/
      DATA R34,W34/ 3.92696350135829D+00, 1.92704402415764D-02/
      DATA R44,W44/ 8.58863568901199D+00, 2.25229076750736D-04/
C
C     X LARGE (ASYMPTOTIC SOLUTION)              NROOTS=4
C
CA    IF(X.GT.53.0D+00) THEN
CA       CALL ROOTSA
CA       RETURN
CA    ENDIF
C
      IF (X .GT. 15.0D+00) GO TO 180
      IF (X .GT. 5.0D+00) GO TO 140
      IF (X .GT. 1.0D+00) GO TO 120
      IF (X .GT. 3.0D-07) GO TO 100
C     X IS APPROXIMATELY ZERO.                   NROOTS = 4
      RT1 = 3.48198973061471D-02 -4.09645850660395D-03 *X
      RT2 = 3.81567185080042D-01 -4.48902570656719D-02 *X
      RT3 = 1.73730726945891D+00 -2.04389090547327D-01 *X
      RT4 = 1.18463056481549D+01 -1.39368301742312D+00 *X
      WW1 = 3.62683783378362D-01 -3.13844305713928D-02 *X
      WW2 = 3.13706645877886D-01 -8.98046242557724D-02 *X
      WW3 = 2.22381034453372D-01 -1.29314370958973D-01 *X
      WW4 = 1.01228536290376D-01 -8.28299075414321D-02 *X
      RETURN
C
C     X=0.0 TO 1.0                               NROOTS = 4
  100 RT1 = ((((((-1.95309614628539D-10*X+5.19765728707592D-09)*X-
     +     1.01756452250573D-07 )*X+1.72365935872131D-06 )*X-
     +     2.61203523522184D-05 )*X+3.52921308769880D-04 )*X-
     +     4.09645850658433D-03 )*X+3.48198973061469D-02
      RT2 = (((((-1.89554881382342D-08*X+3.07583114342365D-07)*X+
     +     1.270981734393D-06)*X-1.417298563884D-04)*X+
     +     3.226979163176D-03)*X-4.48902570678178D-02 )*X+
     +     3.81567185080039D-01
      RT3 = (((((( 1.77280535300416D-09*X+3.36524958870615D-08)*X-
     +     2.58341529013893D-07 )*X-1.13644895662320D-05 )*X-
     +     7.91549618884063D-05 )*X+1.03825827346828D-02 )*X-
     +     2.04389090525137D-01 )*X+1.73730726945889D+00
      RT4 = (((((-5.61188882415248D-08*X-2.49480733072460D-07)*X+
     +     3.428685057114D-06)*X+1.679007454539D-04)*X+
     +     4.722855585715D-02)*X-1.39368301737828D+00 )*X+
     +     1.18463056481543D+01
      WW1 = ((((((-1.14649303201279D-08*X+1.88015570196787D-07)*X-
     +     2.33305875372323D-06 )*X+2.68880044371597D-05 )*X-
     +     2.94268428977387D-04 )*X+3.06548909776613D-03 )*X-
     +     3.13844305680096D-02 )*X+3.62683783378335D-01
      WW2 = ((((((((-4.11720483772634D-09*X+6.54963481852134D-08)*X-
     +     7.20045285129626D-07 )*X+6.93779646721723D-06 )*X-
     +     6.05367572016373D-05 )*X+4.74241566251899D-04 )*X-
     +     3.26956188125316D-03 )*X+1.91883866626681D-02 )*X-
     +     8.98046242565811D-02 )*X+3.13706645877886D-01
      WW3 = ((((((((-3.41688436990215D-08*X+5.07238960340773D-07)*X-
     +     5.01675628408220D-06 )*X+4.20363420922845D-05 )*X-
     +     3.08040221166823D-04 )*X+1.94431864731239D-03 )*X-
     +     1.02477820460278D-02 )*X+4.28670143840073D-02 )*X-
     +     1.29314370962569D-01 )*X+2.22381034453369D-01
      WW4 = ((((((((( 4.99660550769508D-09*X-7.94585963310120D-08)*X+
     +     8.359072409485D-07)*X-7.422369210610D-06)*X+
     +     5.763374308160D-05)*X-3.86645606718233D-04 )*X+
     +     2.18417516259781D-03 )*X-9.99791027771119D-03 )*X+
     +     3.48791097377370D-02 )*X-8.28299075413889D-02 )*X+
     +     1.01228536290376D-01
      RETURN
C
C     X= 1.0 TO 5.0                              NROOTS = 4
  120 Y = X-3.0D+00
      RT1 = (((((((((-1.48570633747284D-15*Y-1.33273068108777D-13)*Y+
     +     4.068543696670D-12)*Y-9.163164161821D-11)*Y+
     +     2.046819017845D-09)*Y-4.03076426299031D-08 )*Y+
     +     7.29407420660149D-07 )*Y-1.23118059980833D-05 )*Y+
     +     1.88796581246938D-04 )*Y-2.53262912046853D-03 )*Y+
     +     2.51198234505021D-02
      RT2 = ((((((((( 1.35830583483312D-13*Y-2.29772605964836D-12)*Y-
     +     3.821500128045D-12)*Y+6.844424214735D-10)*Y-
     +     1.048063352259D-08)*Y+1.50083186233363D-08 )*Y+
     +     3.48848942324454D-06 )*Y-1.08694174399193D-04 )*Y+
     +     2.08048885251999D-03 )*Y-2.91205805373793D-02 )*Y+
     +     2.72276489515713D-01
      RT3 = ((((((((( 5.02799392850289D-13*Y+1.07461812944084D-11)*Y-
     +     1.482277886411D-10)*Y-2.153585661215D-09)*Y+
     +     3.654087802817D-08)*Y+5.15929575830120D-07 )*Y-
     +     9.52388379435709D-06 )*Y-2.16552440036426D-04 )*Y+
     +     9.03551469568320D-03 )*Y-1.45505469175613D-01 )*Y+
     +     1.21449092319186D+00
      RT4 = (((((((((-1.08510370291979D-12*Y+6.41492397277798D-11)*Y+
     +     7.542387436125D-10)*Y-2.213111836647D-09)*Y-
     +     1.448228963549D-07)*Y-1.95670833237101D-06 )*Y-
     +     1.07481314670844D-05 )*Y+1.49335941252765D-04 )*Y+
     +     4.87791531990593D-02 )*Y-1.10559909038653D+00 )*Y+
     +     8.09502028611780D+00
      WW1 = ((((((((((-4.65801912689961D-14*Y+7.58669507106800D-13)*Y-
     +     1.186387548048D-11)*Y+1.862334710665D-10)*Y-
     +     2.799399389539D-09)*Y+4.148972684255D-08)*Y-
     +     5.933568079600D-07)*Y+8.168349266115D-06)*Y-
     +     1.08989176177409D-04 )*Y+1.41357961729531D-03 )*Y-
     +     1.87588361833659D-02 )*Y+2.89898651436026D-01
      WW2 = ((((((((((((-1.46345073267549D-14*Y+2.25644205432182D-13)*Y-
     +     3.116258693847D-12)*Y+4.321908756610D-11)*Y-
     +     5.673270062669D-10)*Y+7.006295962960D-09)*Y-
     +     8.120186517000D-08)*Y+8.775294645770D-07)*Y-
     +     8.77829235749024D-06 )*Y+8.04372147732379D-05 )*Y-
     +     6.64149238804153D-04 )*Y+4.81181506827225D-03 )*Y-
     +     2.88982669486183D-02 )*Y+1.56247249979288D-01
      WW3 = ((((((((((((( 9.06812118895365D-15*Y-1.40541322766087D-13)*
     +     Y+1.919270015269D-12)*Y-2.605135739010D-11)*Y+
     +     3.299685839012D-10)*Y-3.86354139348735D-09 )*Y+
     +     4.16265847927498D-08 )*Y-4.09462835471470D-07 )*Y+
     +     3.64018881086111D-06 )*Y-2.88665153269386D-05 )*Y+
     +     2.00515819789028D-04 )*Y-1.18791896897934D-03 )*Y+
     +     5.75223633388589D-03 )*Y-2.09400418772687D-02 )*Y+
     +     4.85368861938873D-02
      WW4 = ((((((((((((((-9.74835552342257D-16*Y+1.57857099317175D-14)*
     +     Y-2.249993780112D-13)*Y+3.173422008953D-12)*Y-
     +     4.161159459680D-11)*Y+5.021343560166D-10)*Y-
     +     5.545047534808D-09)*Y+5.554146993491D-08)*Y-
     +     4.99048696190133D-07 )*Y+3.96650392371311D-06 )*Y-
     +     2.73816413291214D-05 )*Y+1.60106988333186D-04 )*Y-
     +     7.64560567879592D-04 )*Y+2.81330044426892D-03 )*Y-
     +     7.16227030134947D-03 )*Y+9.66077262223353D-03
      RETURN
C
  140 IF (X .GT. 10.0D+00) GO TO 160
C     X=5.0 TO 10.0                              NROOTS = 4
      Y = X-7.5D+00
      RT1 = ((((((((( 4.64217329776215D-15*Y-6.27892383644164D-15)*Y+
     +     3.462236347446D-13)*Y-2.927229355350D-11)*Y+
     +     5.090355371676D-10)*Y-9.97272656345253D-09 )*Y+
     +     2.37835295639281D-07 )*Y-4.60301761310921D-06 )*Y+
     +     8.42824204233222D-05 )*Y-1.37983082233081D-03 )*Y+
     +     1.66630865869375D-02
      RT2 = ((((((((( 2.93981127919047D-14*Y+8.47635639065744D-13)*Y-
     +     1.446314544774D-11)*Y-6.149155555753D-12)*Y+
     +     8.484275604612D-10)*Y-6.10898827887652D-08 )*Y+
     +     2.39156093611106D-06 )*Y-5.35837089462592D-05 )*Y+
     +     1.00967602595557D-03 )*Y-1.57769317127372D-02 )*Y+
     +     1.74853819464285D-01
      RT3 = (((((((((( 2.93523563363000D-14*Y-6.40041776667020D-14)*Y-
     +     2.695740446312D-12)*Y+1.027082960169D-10)*Y-
     +     5.822038656780D-10)*Y-3.159991002539D-08)*Y+
     +     4.327249251331D-07)*Y+4.856768455119D-06)*Y-
     +     2.54617989427762D-04 )*Y+5.54843378106589D-03 )*Y-
     +     7.95013029486684D-02 )*Y+7.20206142703162D-01
      RT4 = (((((((((((-1.62212382394553D-14*Y+7.68943641360593D-13)*Y+
     +     5.764015756615D-12)*Y-1.380635298784D-10)*Y-
     +     1.476849808675D-09)*Y+1.84347052385605D-08 )*Y+
     +     3.34382940759405D-07 )*Y-1.39428366421645D-06 )*Y-
     +     7.50249313713996D-05 )*Y-6.26495899187507D-04 )*Y+
     +     4.69716410901162D-02 )*Y-6.66871297428209D-01 )*Y+
     +     4.11207530217806D+00
      WW1 = ((((((((((-1.65995045235997D-15*Y+6.91838935879598D-14)*Y-
     +     9.131223418888D-13)*Y+1.403341829454D-11)*Y-
     +     3.672235069444D-10)*Y+6.366962546990D-09)*Y-
     +     1.039220021671D-07)*Y+1.959098751715D-06)*Y-
     +     3.33474893152939D-05 )*Y+5.72164211151013D-04 )*Y-
     +     1.05583210553392D-02 )*Y+2.26696066029591D-01
      WW2 = ((((((((((((-3.57248951192047D-16*Y+6.25708409149331D-15)*Y-
     +     9.657033089714D-14)*Y+1.507864898748D-12)*Y-
     +     2.332522256110D-11)*Y+3.428545616603D-10)*Y-
     +     4.698730937661D-09)*Y+6.219977635130D-08)*Y-
     +     7.83008889613661D-07 )*Y+9.08621687041567D-06 )*Y-
     +     9.86368311253873D-05 )*Y+9.69632496710088D-04 )*Y-
     +     8.14594214284187D-03 )*Y+8.50218447733457D-02
      WW3 = ((((((((((((( 1.64742458534277D-16*Y-2.68512265928410D-15)*
     +     Y+3.788890667676D-14)*Y-5.508918529823D-13)*Y+
     +     7.555896810069D-12)*Y-9.69039768312637D-11 )*Y+
     +     1.16034263529672D-09 )*Y-1.28771698573873D-08 )*Y+
     +     1.31949431805798D-07 )*Y-1.23673915616005D-06 )*Y+
     +     1.04189803544936D-05 )*Y-7.79566003744742D-05 )*Y+
     +     5.03162624754434D-04 )*Y-2.55138844587555D-03 )*Y+
     +     1.13250730954014D-02
      WW4 = ((((((((((((((-1.55714130075679D-17*Y+2.57193722698891D-16)*
     +     Y-3.626606654097D-15)*Y+5.234734676175D-14)*Y-
     +     7.067105402134D-13)*Y+8.793512664890D-12)*Y-
     +     1.006088923498D-10)*Y+1.050565098393D-09)*Y-
     +     9.91517881772662D-09 )*Y+8.35835975882941D-08 )*Y-
     +     6.19785782240693D-07 )*Y+3.95841149373135D-06 )*Y-
     +     2.11366761402403D-05 )*Y+9.00474771229507D-05 )*Y-
     +     2.78777909813289D-04 )*Y+5.26543779837487D-04
      RETURN
C
C     X=10.0 TO 15.0                             NROOTS = 4
  160 Y = X-12.5D+00
      RT1 = ((((((((((( 4.94869622744119D-17*Y+8.03568805739160D-16)*Y-
     +     5.599125915431D-15)*Y-1.378685560217D-13)*Y+
     +     7.006511663249D-13)*Y+1.30391406991118D-11 )*Y+
     +     8.06987313467541D-11 )*Y-5.20644072732933D-09 )*Y+
     +     7.72794187755457D-08 )*Y-1.61512612564194D-06 )*Y+
     +     4.15083811185831D-05 )*Y-7.87855975560199D-04 )*Y+
     +     1.14189319050009D-02
      RT2 = ((((((((((( 4.89224285522336D-16*Y+1.06390248099712D-14)*Y-
     +     5.446260182933D-14)*Y-1.613630106295D-12)*Y+
     +     3.910179118937D-12)*Y+1.90712434258806D-10 )*Y+
     +     8.78470199094761D-10 )*Y-5.97332993206797D-08 )*Y+
     +     9.25750831481589D-07 )*Y-2.02362185197088D-05 )*Y+
     +     4.92341968336776D-04 )*Y-8.68438439874703D-03 )*Y+
     +     1.15825965127958D-01
      RT3 = (((((((((( 6.12419396208408D-14*Y+1.12328861406073D-13)*Y-
     +     9.051094103059D-12)*Y-4.781797525341D-11)*Y+
     +     1.660828868694D-09)*Y+4.499058798868D-10)*Y-
     +     2.519549641933D-07)*Y+4.977444040180D-06)*Y-
     +     1.25858350034589D-04 )*Y+2.70279176970044D-03 )*Y-
     +     3.99327850801083D-02 )*Y+4.33467200855434D-01
      RT4 = ((((((((((( 4.63414725924048D-14*Y-4.72757262693062D-14)*Y-
     +     1.001926833832D-11)*Y+6.074107718414D-11)*Y+
     +     1.576976911942D-09)*Y-2.01186401974027D-08 )*Y-
     +     1.84530195217118D-07 )*Y+5.02333087806827D-06 )*Y+
     +     9.66961790843006D-06 )*Y-1.58522208889528D-03 )*Y+
     +     2.80539673938339D-02 )*Y-2.78953904330072D-01 )*Y+
     +     1.82835655238235D+00
      WW4 = ((((((((((((( 2.90401781000996D-18*Y-4.63389683098251D-17)*
     +     Y+6.274018198326D-16)*Y-8.936002188168D-15)*Y+
     +     1.194719074934D-13)*Y-1.45501321259466D-12 )*Y+
     +     1.64090830181013D-11 )*Y-1.71987745310181D-10 )*Y+
     +     1.63738403295718D-09 )*Y-1.39237504892842D-08 )*Y+
     +     1.06527318142151D-07 )*Y-7.27634957230524D-07 )*Y+
     +     4.12159381310339D-06 )*Y-1.74648169719173D-05 )*Y+
     +     8.50290130067818D-05
      WW3 = ((((((((((((-4.19569145459480D-17*Y+5.94344180261644D-16)*Y-
     +     1.148797566469D-14)*Y+1.881303962576D-13)*Y-
     +     2.413554618391D-12)*Y+3.372127423047D-11)*Y-
     +     4.933988617784D-10)*Y+6.116545396281D-09)*Y-
     +     6.69965691739299D-08 )*Y+7.52380085447161D-07 )*Y-
     +     8.08708393262321D-06 )*Y+6.88603417296672D-05 )*Y-
     +     4.67067112993427D-04 )*Y+5.42313365864597D-03
      WW2 = ((((((((((-6.22272689880615D-15*Y+1.04126809657554D-13)*Y-
     +     6.842418230913D-13)*Y+1.576841731919D-11)*Y-
     +     4.203948834175D-10)*Y+6.287255934781D-09)*Y-
     +     8.307159819228D-08)*Y+1.356478091922D-06)*Y-
     +     2.08065576105639D-05 )*Y+2.52396730332340D-04 )*Y-
     +     2.94484050194539D-03 )*Y+6.01396183129168D-02
      WW1 = (((-1.8784686463512D-01/X+2.2991849164985D-01)/X -
     +     4.9893752514047D-01)/X-2.1916512131607D-05)*EXP(-X) +
     +     SQRT(PIE4/X)-WW4-WW3-WW2
      RETURN
C
  180 WW1 = SQRT(PIE4/X)
      IF (X .GT. 35.0D+00) GO TO 220
      IF (X .GT. 20.0D+00) GO TO 200
C     X=15.0 TO 20.0                             NROOTS = 4
      Y = X-17.5D+00
      RT1 = ((((((((((( 4.36701759531398D-17*Y-1.12860600219889D-16)*Y-
     +     6.149849164164D-15)*Y+5.820231579541D-14)*Y+
     +     4.396602872143D-13)*Y-1.24330365320172D-11 )*Y+
     +     6.71083474044549D-11 )*Y+2.43865205376067D-10 )*Y+
     +     1.67559587099969D-08 )*Y-9.32738632357572D-07 )*Y+
     +     2.39030487004977D-05 )*Y-4.68648206591515D-04 )*Y+
     +     8.34977776583956D-03
      RT2 = ((((((((((( 4.98913142288158D-16*Y-2.60732537093612D-16)*Y-
     +     7.775156445127D-14)*Y+5.766105220086D-13)*Y+
     +     6.432696729600D-12)*Y-1.39571683725792D-10 )*Y+
     +     5.95451479522191D-10 )*Y+2.42471442836205D-09 )*Y+
     +     2.47485710143120D-07 )*Y-1.14710398652091D-05 )*Y+
     +     2.71252453754519D-04 )*Y-4.96812745851408D-03 )*Y+
     +     8.26020602026780D-02
      RT3 = ((((((((((( 1.91498302509009D-15*Y+1.48840394311115D-14)*Y-
     +     4.316925145767D-13)*Y+1.186495793471D-12)*Y+
     +     4.615806713055D-11)*Y-5.54336148667141D-10 )*Y+
     +     3.48789978951367D-10 )*Y-2.79188977451042D-09 )*Y+
     +     2.09563208958551D-06 )*Y-6.76512715080324D-05 )*Y+
     +     1.32129867629062D-03 )*Y-2.05062147771513D-02 )*Y+
     +     2.88068671894324D-01
      RT4 = (((((((((((-5.43697691672942D-15*Y-1.12483395714468D-13)*Y+
     +     2.826607936174D-12)*Y-1.266734493280D-11)*Y-
     +     4.258722866437D-10)*Y+9.45486578503261D-09 )*Y-
     +     5.86635622821309D-08 )*Y-1.28835028104639D-06 )*Y+
     +     4.41413815691885D-05 )*Y-7.61738385590776D-04 )*Y+
     +     9.66090902985550D-03 )*Y-1.01410568057649D-01 )*Y+
     +     9.54714798156712D-01
      WW4 = ((((((((((((-7.56882223582704D-19*Y+7.53541779268175D-18)*Y-
     +     1.157318032236D-16)*Y+2.411195002314D-15)*Y-
     +     3.601794386996D-14)*Y+4.082150659615D-13)*Y-
     +     4.289542980767D-12)*Y+5.086829642731D-11)*Y-
     +     6.35435561050807D-10 )*Y+6.82309323251123D-09 )*Y-
     +     5.63374555753167D-08 )*Y+3.57005361100431D-07 )*Y-
     +     2.40050045173721D-06 )*Y+4.94171300536397D-05
      WW3 = (((((((((((-5.54451040921657D-17*Y+2.68748367250999D-16)*Y+
     +     1.349020069254D-14)*Y-2.507452792892D-13)*Y+
     +     1.944339743818D-12)*Y-1.29816917658823D-11 )*Y+
     +     3.49977768819641D-10 )*Y-8.67270669346398D-09 )*Y+
     +     1.31381116840118D-07 )*Y-1.36790720600822D-06 )*Y+
     +     1.19210697673160D-05 )*Y-1.42181943986587D-04 )*Y+
     +     4.12615396191829D-03
      WW2 = (((((((((((-1.86506057729700D-16*Y+1.16661114435809D-15)*Y+
     +     2.563712856363D-14)*Y-4.498350984631D-13)*Y+
     +     1.765194089338D-12)*Y+9.04483676345625D-12 )*Y+
     +     4.98930345609785D-10 )*Y-2.11964170928181D-08 )*Y+
     +     3.98295476005614D-07 )*Y-5.49390160829409D-06 )*Y+
     +     7.74065155353262D-05 )*Y-1.48201933009105D-03 )*Y+
     +     4.97836392625268D-02
      WW1 = (( 1.9623264149430D-01/X-4.9695241464490D-01)/X -
     +     6.0156581186481D-05)*EXP(-X)+WW1-WW2-WW3-WW4
      RETURN
C
C     X=20.0 TO 35.0                             NROOTS = 4
  200 E = EXP(-X)
      RT1 = ((((((-4.45711399441838D-05*X+1.27267770241379D-03)*X -
     +     2.36954961381262D-01)*X+1.54330657903756D+01)*X -
     +     5.22799159267808D+02)*X+1.05951216669313D+04)*X + (-
     +     2.51177235556236D+06/X+8.72975373557709D+05)/X -
     +     1.29194382386499D+05)*E + R14/(X-R14)
      RT2 = (((((-7.85617372254488D-02*X+6.35653573484868D+00)*X -
     +     3.38296938763990D+02)*X+1.25120495802096D+04)*X -
     +     3.16847570511637D+05)*X + ((-1.02427466127427D+09/X +
     +     3.70104713293016D+08)/X-5.87119005093822D+07)/X +
     +     5.38614211391604D+06)*E + R24/(X-R24)
      RT3 = (((((-2.37900485051067D-01*X+1.84122184400896D+01)*X -
     +     1.00200731304146D+03)*X+3.75151841595736D+04)*X -
     +     9.50626663390130D+05)*X + ((-2.88139014651985D+09/X +
     +     1.06625915044526D+09)/X-1.72465289687396D+08)/X +
     +     1.60419390230055D+07)*E + R34/(X-R34)
      RT4 = ((((((-6.00691586407385D-04*X-3.64479545338439D-01)*X +
     +     1.57496131755179D+01)*X-6.54944248734901D+02)*X +
     +     1.70830039597097D+04)*X-2.90517939780207D+05)*X + (+
     +     3.49059698304732D+07/X-1.64944522586065D+07)/X +
     +     2.96817940164703D+06)*E + R44/(X-R44)
      IF (X .LE. 25.0D+00) WW4 = ((((((( 2.33766206773151D-07*X-
     +     3.81542906607063D-05)*X +3.51416601267000D-03)*X-
     +     1.66538571864728D-01)*X +4.80006136831847D+00)*X-
     +     8.73165934223603D+01)*X +9.77683627474638D+02)*X +
     +     1.66000945117640D+04/X -6.14479071209961D+03)*E + W44*WW1
      IF (X .GT. 25.0D+00) WW4 = (((((( 5.74245945342286D-06*X-
     +     7.58735928102351D-05)*X +2.35072857922892D-04)*X-
     +     3.78812134013125D-03)*X +3.09871652785805D-01)*X-
     +     7.11108633061306D+00)*X +5.55297573149528D+01)*E + W44*WW1
      WW3 = (((((( 2.36392855180768D-04*X-9.16785337967013D-03)*X +
     +     4.62186525041313D-01)*X-1.96943786006540D+01)*X +
     +     4.99169195295559D+02)*X-6.21419845845090D+03)*X + ((+
     +     5.21445053212414D+07/X-1.34113464389309D+07)/X +
     +     1.13673298305631D+06)/X-2.81501182042707D+03)*E + W34*WW1
      WW2 = (((((( 7.29841848989391D-04*X-3.53899555749875D-02)*X +
     +     2.07797425718513D+00)*X-1.00464709786287D+02)*X +
     +     3.15206108877819D+03)*X-6.27054715090012D+04)*X + (+
     +     1.54721246264919D+07/X-5.26074391316381D+06)/X +
     +     7.67135400969617D+05)*E + W24*WW1
      WW1 = (( 1.9623264149430D-01/X-4.9695241464490D-01)/X -
     +     6.0156581186481D-05)*E + WW1-WW2-WW3-WW4
      RETURN
C
  220 IF (X .GT. 53.0D+00) GO TO 240
C     X=35.0 TO 53.0                             NROOTS = 4
      E = EXP(-X)*(X*X)**2
      RT4 = ((-2.19135070169653D-03*X-1.19108256987623D-01)*X -
     +     7.50238795695573D-01)*E + R44/(X-R44)
      RT3 = ((-9.65842534508637D-04*X-4.49822013469279D-02)*X +
     +     6.08784033347757D-01)*E + R34/(X-R34)
      RT2 = ((-3.62569791162153D-04*X-9.09231717268466D-03)*X +
     +     1.84336760556262D-01)*E + R24/(X-R24)
      RT1 = ((-4.07557525914600D-05*X-6.88846864931685D-04)*X +
     +     1.74725309199384D-02)*E + R14/(X-R14)
      WW4 = (( 5.76631982000990D-06*X-7.89187283804890D-05)*X +
     +     3.28297971853126D-04)*E + W44*WW1
      WW3 = (( 2.08294969857230D-04*X-3.77489954837361D-03)*X +
     +     2.09857151617436D-02)*E + W34*WW1
      WW2 = (( 6.16374517326469D-04*X-1.26711744680092D-02)*X +
     +     8.14504890732155D-02)*E + W24*WW1
      WW1 = WW1-WW2-WW3-WW4
      RETURN
C
C     X=47.0 TO INFINITY                         NROOTS = 4
  240 RT1 = R14/(X-R14)
      RT2 = R24/(X-R24)
      RT3 = R34/(X-R34)
      RT4 = R44/(X-R44)
      WW4 = W44*WW1
      WW3 = W34*WW1
      WW2 = W24*WW1
      WW1 = WW1-WW2-WW3-WW4
      RETURN
      END
C*MODULE RYSPOL  *DECK ROOT5
      SUBROUTINE ROOT5
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C          *****   VERSION  FEBRUARY 27,1975   *****
C
      COMMON /ROOT  / X,U(13),W(13),NROOTS
C
      EQUIVALENCE (U(1),RT1),(U(2),RT2),(U(3),RT3),(U(4),RT4),(U(5),RT5)
      EQUIVALENCE (W(1),WW1),(W(2),WW2),(W(3),WW3),(W(4),WW4),(W(5),WW5)
C
      DATA R15,PIE4/1.17581320211778D-01, 7.85398163397448D-01/
      DATA R25,W25/ 1.07456201243690D+00, 2.70967405960535D-01/
      DATA R35,W35/ 3.08593744371754D+00, 3.82231610015404D-02/
      DATA R45,W45/ 6.41472973366203D+00, 1.51614186862443D-03/
      DATA R55,W55/ 1.18071894899717D+01, 8.62130526143657D-06/
C
C     X LARGE (ASYMPTOTIC SOLUTION)              NROOTS=5
C
CA    IF(X.GT.59.0D+00) THEN
CA       CALL ROOTSA
CA       RETURN
CA    ENDIF
C
      IF (X .GT. 15.0D+00) GO TO 180
      IF (X .GT. 5.0D+00) GO TO 140
      IF (X .GT. 1.0D+00) GO TO 120
      IF (X .GT. 3.0D-07) GO TO 100
C     X IS APPROXIMATELY ZERO.                   NROOTS = 5
      RT1 = 2.26659266316985D-02 -2.15865967920897D-03 *X
      RT2 = 2.31271692140903D-01 -2.20258754389745D-02 *X
      RT3 = 8.57346024118836D-01 -8.16520023025515D-02 *X
      RT4 = 2.97353038120346D+00 -2.83193369647137D-01 *X
      RT5 = 1.84151859759051D+01 -1.75382723579439D+00 *X
      WW1 = 2.95524224714752D-01 -1.96867576909777D-02 *X
      WW2 = 2.69266719309995D-01 -5.61737590184721D-02 *X
      WW3 = 2.19086362515981D-01 -9.71152726793658D-02 *X
      WW4 = 1.49451349150580D-01 -1.02979262193565D-01 *X
      WW5 = 6.66713443086877D-02 -5.73782817488315D-02 *X
      RETURN
C
C     X=0.0 TO 1.0                               NROOTS = 5
  100 RT1 = ((((((-4.46679165328413D-11*X+1.21879111988031D-09)*X-
     +     2.62975022612104D-08 )*X+5.15106194905897D-07 )*X-
     +     9.27933625824749D-06 )*X+1.51794097682482D-04 )*X-
     +     2.15865967920301D-03 )*X+2.26659266316985D-02
      RT2 = (((((( 1.93117331714174D-10*X-4.57267589660699D-09)*X+
     +     2.48339908218932D-08 )*X+1.50716729438474D-06 )*X-
     +     6.07268757707381D-05 )*X+1.37506939145643D-03 )*X-
     +     2.20258754419939D-02 )*X+2.31271692140905D-01
      RT3 = ((((( 4.84989776180094D-09*X+1.31538893944284D-07)*X-
     +     2.766753852879D-06)*X-7.651163510626D-05)*X+
     +     4.033058545972D-03)*X-8.16520022916145D-02 )*X+
     +     8.57346024118779D-01
      RT4 = ((((-2.48581772214623D-07*X-4.34482635782585D-06)*X-
     +     7.46018257987630D-07 )*X+1.01210776517279D-02 )*X-
     +     2.83193369640005D-01 )*X+2.97353038120345D+00
      RT5 = (((((-8.92432153868554D-09*X+1.77288899268988D-08)*X+
     +     3.040754680666D-06)*X+1.058229325071D-04)*X+
     +     4.596379534985D-02)*X-1.75382723579114D+00 )*X+
     +     1.84151859759049D+01
      WW1 = ((((((-2.03822632771791D-09*X+3.89110229133810D-08)*X-
     +     5.84914787904823D-07 )*X+8.30316168666696D-06 )*X-
     +     1.13218402310546D-04 )*X+1.49128888586790D-03 )*X-
     +     1.96867576904816D-02 )*X+2.95524224714749D-01
      WW2 = ((((((( 8.62848118397570D-09*X-1.38975551148989D-07)*X+
     +     1.602894068228D-06)*X-1.646364300836D-05)*X+
     +     1.538445806778D-04)*X-1.28848868034502D-03 )*X+
     +     9.38866933338584D-03 )*X-5.61737590178812D-02 )*X+
     +     2.69266719309991D-01
      WW3 = ((((((((-9.41953204205665D-09*X+1.47452251067755D-07)*X-
     +     1.57456991199322D-06 )*X+1.45098401798393D-05 )*X-
     +     1.18858834181513D-04 )*X+8.53697675984210D-04 )*X-
     +     5.22877807397165D-03 )*X+2.60854524809786D-02 )*X-
     +     9.71152726809059D-02 )*X+2.19086362515979D-01
      WW4 = ((((((((-3.84961617022042D-08*X+5.66595396544470D-07)*X-
     +     5.52351805403748D-06 )*X+4.53160377546073D-05 )*X-
     +     3.22542784865557D-04 )*X+1.95682017370967D-03 )*X-
     +     9.77232537679229D-03 )*X+3.79455945268632D-02 )*X-
     +     1.02979262192227D-01 )*X+1.49451349150573D-01
      WW5 = ((((((((( 4.09594812521430D-09*X-6.47097874264417D-08)*X+
     +     6.743541482689D-07)*X-5.917993920224D-06)*X+
     +     4.531969237381D-05)*X-2.99102856679638D-04 )*X+
     +     1.65695765202643D-03 )*X-7.40671222520653D-03 )*X+
     +     2.50889946832192D-02 )*X-5.73782817487958D-02 )*X+
     +     6.66713443086877D-02
      RETURN
C
C     X=1.0 TO 5.0                               NROOTS = 5
  120 Y = X-3.0D+00
      RT1 = ((((((((-2.58163897135138D-14*Y+8.14127461488273D-13)*Y-
     +     2.11414838976129D-11 )*Y+5.09822003260014D-10 )*Y-
     +     1.16002134438663D-08 )*Y+2.46810694414540D-07 )*Y-
     +     4.92556826124502D-06 )*Y+9.02580687971053D-05 )*Y-
     +     1.45190025120726D-03 )*Y+1.73416786387475D-02
      RT2 = ((((((((( 1.04525287289788D-14*Y+5.44611782010773D-14)*Y-
     +     4.831059411392D-12)*Y+1.136643908832D-10)*Y-
     +     1.104373076913D-09)*Y-2.35346740649916D-08 )*Y+
     +     1.43772622028764D-06 )*Y-4.23405023015273D-05 )*Y+
     +     9.12034574793379D-04 )*Y-1.52479441718739D-02 )*Y+
     +     1.76055265928744D-01
      RT3 = (((((((((-6.89693150857911D-14*Y+5.92064260918861D-13)*Y+
     +     1.847170956043D-11)*Y-3.390752744265D-10)*Y-
     +     2.995532064116D-09)*Y+1.57456141058535D-07 )*Y-
     +     3.95859409711346D-07 )*Y-9.58924580919747D-05 )*Y+
     +     3.23551502557785D-03 )*Y-5.97587007636479D-02 )*Y+
     +     6.46432853383057D-01
      RT4 = ((((((((-3.61293809667763D-12*Y-2.70803518291085D-11)*Y+
     +     8.83758848468769D-10 )*Y+1.59166632851267D-08 )*Y-
     +     1.32581997983422D-07 )*Y-7.60223407443995D-06 )*Y-
     +     7.41019244900952D-05 )*Y+9.81432631743423D-03 )*Y-
     +     2.23055570487771D-01 )*Y+2.21460798080643D+00
      RT5 = ((((((((( 7.12332088345321D-13*Y+3.16578501501894D-12)*Y-
     +     8.776668218053D-11)*Y-2.342817613343D-09)*Y-
     +     3.496962018025D-08)*Y-3.03172870136802D-07 )*Y+
     +     1.50511293969805D-06 )*Y+1.37704919387696D-04 )*Y+
     +     4.70723869619745D-02 )*Y-1.47486623003693D+00 )*Y+
     +     1.35704792175847D+01
      WW1 = ((((((((( 1.04348658616398D-13*Y-1.94147461891055D-12)*Y+
     +     3.485512360993D-11)*Y-6.277497362235D-10)*Y+
     +     1.100758247388D-08)*Y-1.88329804969573D-07 )*Y+
     +     3.12338120839468D-06 )*Y-5.04404167403568D-05 )*Y+
     +     8.00338056610995D-04 )*Y-1.30892406559521D-02 )*Y+
     +     2.47383140241103D-01
      WW2 = ((((((((((( 3.23496149760478D-14*Y-5.24314473469311D-13)*Y+
     +     7.743219385056D-12)*Y-1.146022750992D-10)*Y+
     +     1.615238462197D-09)*Y-2.15479017572233D-08 )*Y+
     +     2.70933462557631D-07 )*Y-3.18750295288531D-06 )*Y+
     +     3.47425221210099D-05 )*Y-3.45558237388223D-04 )*Y+
     +     3.05779768191621D-03 )*Y-2.29118251223003D-02 )*Y+
     +     1.59834227924213D-01
      WW3 = ((((((((((((-3.42790561802876D-14*Y+5.26475736681542D-13)*Y-
     +     7.184330797139D-12)*Y+9.763932908544D-11)*Y-
     +     1.244014559219D-09)*Y+1.472744068942D-08)*Y-
     +     1.611749975234D-07)*Y+1.616487851917D-06)*Y-
     +     1.46852359124154D-05 )*Y+1.18900349101069D-04 )*Y-
     +     8.37562373221756D-04 )*Y+4.93752683045845D-03 )*Y-
     +     2.25514728915673D-02 )*Y+6.95211812453929D-02
      WW4 = ((((((((((((( 1.04072340345039D-14*Y-1.60808044529211D-13)*
     +     Y+2.183534866798D-12)*Y-2.939403008391D-11)*Y+
     +     3.679254029085D-10)*Y-4.23775673047899D-09 )*Y+
     +     4.46559231067006D-08 )*Y-4.26488836563267D-07 )*Y+
     +     3.64721335274973D-06 )*Y-2.74868382777722D-05 )*Y+
     +     1.78586118867488D-04 )*Y-9.68428981886534D-04 )*Y+
     +     4.16002324339929D-03 )*Y-1.28290192663141D-02 )*Y+
     +     2.22353727685016D-02
      WW5 = ((((((((((((((-8.16770412525963D-16*Y+1.31376515047977D-14)*
     +     Y-1.856950818865D-13)*Y+2.596836515749D-12)*Y-
     +     3.372639523006D-11)*Y+4.025371849467D-10)*Y-
     +     4.389453269417D-09)*Y+4.332753856271D-08)*Y-
     +     3.82673275931962D-07 )*Y+2.98006900751543D-06 )*Y-
     +     2.00718990300052D-05 )*Y+1.13876001386361D-04 )*Y-
     +     5.23627942443563D-04 )*Y+1.83524565118203D-03 )*Y-
     +     4.37785737450783D-03 )*Y+5.36963805223095D-03
      RETURN
C
  140 IF (X .GT. 10.0D+00) GO TO 160
C     X=5.0 TO 10.0                              NROOTS = 5
      Y = X-7.5D+00
      RT1 = ((((((((-1.13825201010775D-14*Y+1.89737681670375D-13)*Y-
     +     4.81561201185876D-12 )*Y+1.56666512163407D-10 )*Y-
     +     3.73782213255083D-09 )*Y+9.15858355075147D-08 )*Y-
     +     2.13775073585629D-06 )*Y+4.56547356365536D-05 )*Y-
     +     8.68003909323740D-04 )*Y+1.22703754069176D-02
      RT2 = (((((((((-3.67160504428358D-15*Y+1.27876280158297D-14)*Y-
     +     1.296476623788D-12)*Y+1.477175434354D-11)*Y+
     +     5.464102147892D-10)*Y-2.42538340602723D-08 )*Y+
     +     8.20460740637617D-07 )*Y-2.20379304598661D-05 )*Y+
     +     4.90295372978785D-04 )*Y-9.14294111576119D-03 )*Y+
     +     1.22590403403690D-01
      RT3 = ((((((((( 1.39017367502123D-14*Y-6.96391385426890D-13)*Y+
     +     1.176946020731D-12)*Y+1.725627235645D-10)*Y-
     +     3.686383856300D-09)*Y+2.87495324207095D-08 )*Y+
     +     1.71307311000282D-06 )*Y-7.94273603184629D-05 )*Y+
     +     2.00938064965897D-03 )*Y-3.63329491677178D-02 )*Y+
     +     4.34393683888443D-01
      RT4 = ((((((((((-1.27815158195209D-14*Y+1.99910415869821D-14)*Y+
     +     3.753542914426D-12)*Y-2.708018219579D-11)*Y-
     +     1.190574776587D-09)*Y+1.106696436509D-08)*Y+
     +     3.954955671326D-07)*Y-4.398596059588D-06)*Y-
     +     2.01087998907735D-04 )*Y+7.89092425542937D-03 )*Y-
     +     1.42056749162695D-01 )*Y+1.39964149420683D+00
      RT5 = ((((((((((-1.19442341030461D-13*Y-2.34074833275956D-12)*Y+
     +     6.861649627426D-12)*Y+6.082671496226D-10)*Y+
     +     5.381160105420D-09)*Y-6.253297138700D-08)*Y-
     +     2.135966835050D-06)*Y-2.373394341886D-05)*Y+
     +     2.88711171412814D-06 )*Y+4.85221195290753D-02 )*Y-
     +     1.04346091985269D+00 )*Y+7.89901551676692D+00
      WW1 = ((((((((( 7.95526040108997D-15*Y-2.48593096128045D-13)*Y+
     +     4.761246208720D-12)*Y-9.535763686605D-11)*Y+
     +     2.225273630974D-09)*Y-4.49796778054865D-08 )*Y+
     +     9.17812870287386D-07 )*Y-1.86764236490502D-05 )*Y+
     +     3.76807779068053D-04 )*Y-8.10456360143408D-03 )*Y+
     +     2.01097936411496D-01
      WW2 = ((((((((((( 1.25678686624734D-15*Y-2.34266248891173D-14)*Y+
     +     3.973252415832D-13)*Y-6.830539401049D-12)*Y+
     +     1.140771033372D-10)*Y-1.82546185762009D-09 )*Y+
     +     2.77209637550134D-08 )*Y-4.01726946190383D-07 )*Y+
     +     5.48227244014763D-06 )*Y-6.95676245982121D-05 )*Y+
     +     8.05193921815776D-04 )*Y-8.15528438784469D-03 )*Y+
     +     9.71769901268114D-02
      WW3 = ((((((((((((-8.20929494859896D-16*Y+1.37356038393016D-14)*Y-
     +     2.022863065220D-13)*Y+3.058055403795D-12)*Y-
     +     4.387890955243D-11)*Y+5.923946274445D-10)*Y-
     +     7.503659964159D-09)*Y+8.851599803902D-08)*Y-
     +     9.65561998415038D-07 )*Y+9.60884622778092D-06 )*Y-
     +     8.56551787594404D-05 )*Y+6.66057194311179D-04 )*Y-
     +     4.17753183902198D-03 )*Y+2.25443826852447D-02
      WW4 = ((((((((((((((-1.08764612488790D-17*Y+1.85299909689937D-16)*
     +     Y-2.730195628655D-15)*Y+4.127368817265D-14)*Y-
     +     5.881379088074D-13)*Y+7.805245193391D-12)*Y-
     +     9.632707991704D-11)*Y+1.099047050624D-09)*Y-
     +     1.15042731790748D-08 )*Y+1.09415155268932D-07 )*Y-
     +     9.33687124875935D-07 )*Y+7.02338477986218D-06 )*Y-
     +     4.53759748787756D-05 )*Y+2.41722511389146D-04 )*Y-
     +     9.75935943447037D-04 )*Y+2.57520532789644D-03
      WW5 = ((((((((((((((( 7.28996979748849D-19*Y-1.26518146195173D-17)
     +     *Y+1.886145834486D-16)*Y-2.876728287383D-15)*Y+
     +     4.114588668138D-14)*Y-5.44436631413933D-13 )*Y+
     +     6.64976446790959D-12 )*Y-7.44560069974940D-11 )*Y+
     +     7.57553198166848D-10 )*Y-6.92956101109829D-09 )*Y+
     +     5.62222859033624D-08 )*Y-3.97500114084351D-07 )*Y+
     +     2.39039126138140D-06 )*Y-1.18023950002105D-05 )*Y+
     +     4.52254031046244D-05 )*Y-1.21113782150370D-04 )*Y+
     +     1.75013126731224D-04
      RETURN
C
C     X=10.0 TO 15.0                             NROOTS = 5
  160 Y = X-12.5D+00
      RT1 = ((((((((((-4.16387977337393D-17*Y+7.20872997373860D-16)*Y+
     +     1.395993802064D-14)*Y+3.660484641252D-14)*Y-
     +     4.154857548139D-12)*Y+2.301379846544D-11)*Y-
     +     1.033307012866D-09)*Y+3.997777641049D-08)*Y-
     +     9.35118186333939D-07 )*Y+2.38589932752937D-05 )*Y-
     +     5.35185183652937D-04 )*Y+8.85218988709735D-03
      RT2 = ((((((((((-4.56279214732217D-16*Y+6.24941647247927D-15)*Y+
     +     1.737896339191D-13)*Y+8.964205979517D-14)*Y-
     +     3.538906780633D-11)*Y+9.561341254948D-11)*Y-
     +     9.772831891310D-09)*Y+4.240340194620D-07)*Y-
     +     1.02384302866534D-05 )*Y+2.57987709704822D-04 )*Y-
     +     5.54735977651677D-03 )*Y+8.68245143991948D-02
      RT3 = ((((((((((-2.52879337929239D-15*Y+2.13925810087833D-14)*Y+
     +     7.884307667104D-13)*Y-9.023398159510D-13)*Y-
     +     5.814101544957D-11)*Y-1.333480437968D-09)*Y-
     +     2.217064940373D-08)*Y+1.643290788086D-06)*Y-
     +     4.39602147345028D-05 )*Y+1.08648982748911D-03 )*Y-
     +     2.13014521653498D-02 )*Y+2.94150684465425D-01
      RT4 = ((((((((((-6.42391438038888D-15*Y+5.37848223438815D-15)*Y+
     +     8.960828117859D-13)*Y+5.214153461337D-11)*Y-
     +     1.106601744067D-10)*Y-2.007890743962D-08)*Y+
     +     1.543764346501D-07)*Y+4.520749076914D-06)*Y-
     +     1.88893338587047D-04 )*Y+4.73264487389288D-03 )*Y-
     +     7.91197893350253D-02 )*Y+8.60057928514554D-01
      RT5 = (((((((((((-2.24366166957225D-14*Y+4.87224967526081D-14)*Y+
     +     5.587369053655D-12)*Y-3.045253104617D-12)*Y-
     +     1.223983883080D-09)*Y-2.05603889396319D-09 )*Y+
     +     2.58604071603561D-07 )*Y+1.34240904266268D-06 )*Y-
     +     5.72877569731162D-05 )*Y-9.56275105032191D-04 )*Y+
     +     4.23367010370921D-02 )*Y-5.76800927133412D-01 )*Y+
     +     3.87328263873381D+00
      WW1 = ((((((((( 8.98007931950169D-15*Y+7.25673623859497D-14)*Y+
     +     5.851494250405D-14)*Y-4.234204823846D-11)*Y+
     +     3.911507312679D-10)*Y-9.65094802088511D-09 )*Y+
     +     3.42197444235714D-07 )*Y-7.51821178144509D-06 )*Y+
     +     1.94218051498662D-04 )*Y-5.38533819142287D-03 )*Y+
     +     1.68122596736809D-01
      WW2 = ((((((((((-1.05490525395105D-15*Y+1.96855386549388D-14)*Y-
     +     5.500330153548D-13)*Y+1.003849567976D-11)*Y-
     +     1.720997242621D-10)*Y+3.533277061402D-09)*Y-
     +     6.389171736029D-08)*Y+1.046236652393D-06)*Y-
     +     1.73148206795827D-05 )*Y+2.57820531617185D-04 )*Y-
     +     3.46188265338350D-03 )*Y+7.03302497508176D-02
      WW3 = ((((((((((( 3.60020423754545D-16*Y-6.24245825017148D-15)*Y+
     +     9.945311467434D-14)*Y-1.749051512721D-12)*Y+
     +     2.768503957853D-11)*Y-4.08688551136506D-10 )*Y+
     +     6.04189063303610D-09 )*Y-8.23540111024147D-08 )*Y+
     +     1.01503783870262D-06 )*Y-1.20490761741576D-05 )*Y+
     +     1.26928442448148D-04 )*Y-1.05539461930597D-03 )*Y+
     +     1.15543698537013D-02
      WW4 = ((((((((((((( 2.51163533058925D-18*Y-4.31723745510697D-17)*
     +     Y+6.557620865832D-16)*Y-1.016528519495D-14)*Y+
     +     1.491302084832D-13)*Y-2.06638666222265D-12 )*Y+
     +     2.67958697789258D-11 )*Y-3.23322654638336D-10 )*Y+
     +     3.63722952167779D-09 )*Y-3.75484943783021D-08 )*Y+
     +     3.49164261987184D-07 )*Y-2.92658670674908D-06 )*Y+
     +     2.12937256719543D-05 )*Y-1.19434130620929D-04 )*Y+
     +     6.45524336158384D-04
      WW5 = ((((((((((((((-1.29043630202811D-19*Y+2.16234952241296D-18)*
     +     Y-3.107631557965D-17)*Y+4.570804313173D-16)*Y-
     +     6.301348858104D-15)*Y+8.031304476153D-14)*Y-
     +     9.446196472547D-13)*Y+1.018245804339D-11)*Y-
     +     9.96995451348129D-11 )*Y+8.77489010276305D-10 )*Y-
     +     6.84655877575364D-09 )*Y+4.64460857084983D-08 )*Y-
     +     2.66924538268397D-07 )*Y+1.24621276265907D-06 )*Y-
     +     4.30868944351523D-06 )*Y+9.94307982432868D-06
      RETURN
C
  180 IF (X .GT. 25.0D+00) GO TO 220
      IF (X .GT. 20.0D+00) GO TO 200
C     X=15.0 TO 20.0                             NROOTS = 5
      Y = X-17.5D+00
      RT1 = (((((((((( 1.91875764545740D-16*Y+7.8357401095707D-16)*Y-
     +     3.260875931644D-14)*Y-1.186752035569D-13)*Y+
     +     4.275180095653D-12)*Y+3.357056136731D-11)*Y-
     +     1.123776903884D-09)*Y+1.231203269887D-08)*Y-
     +     3.99851421361031D-07 )*Y+1.45418822817771D-05 )*Y-
     +     3.49912254976317D-04 )*Y+6.67768703938812D-03
      RT2 = (((((((((( 2.02778478673555D-15*Y+1.01640716785099D-14)*Y-
     +     3.385363492036D-13)*Y-1.615655871159D-12)*Y+
     +     4.527419140333D-11)*Y+3.853670706486D-10)*Y-
     +     1.184607130107D-08)*Y+1.347873288827D-07)*Y-
     +     4.47788241748377D-06 )*Y+1.54942754358273D-04 )*Y-
     +     3.55524254280266D-03 )*Y+6.44912219301603D-02
      RT3 = (((((((((( 7.79850771456444D-15*Y+6.00464406395001D-14)*Y-
     +     1.249779730869D-12)*Y-1.020720636353D-11)*Y+
     +     1.814709816693D-10)*Y+1.766397336977D-09)*Y-
     +     4.603559449010D-08)*Y+5.863956443581D-07)*Y-
     +     2.03797212506691D-05 )*Y+6.31405161185185D-04 )*Y-
     +     1.30102750145071D-02 )*Y+2.10244289044705D-01
      RT4 = (((((((((((-2.92397030777912D-15*Y+1.94152129078465D-14)*Y+
     +     4.859447665850D-13)*Y-3.217227223463D-12)*Y-
     +     7.484522135512D-11)*Y+7.19101516047753D-10 )*Y+
     +     6.88409355245582D-09 )*Y-1.44374545515769D-07 )*Y+
     +     2.74941013315834D-06 )*Y-1.02790452049013D-04 )*Y+
     +     2.59924221372643D-03 )*Y-4.35712368303551D-02 )*Y+
     +     5.62170709585029D-01
      RT5 = ((((((((((( 1.17976126840060D-14*Y+1.24156229350669D-13)*Y-
     +     3.892741622280D-12)*Y-7.755793199043D-12)*Y+
     +     9.492190032313D-10)*Y-4.98680128123353D-09 )*Y-
     +     1.81502268782664D-07 )*Y+2.69463269394888D-06 )*Y+
     +     2.50032154421640D-05 )*Y-1.33684303917681D-03 )*Y+
     +     2.29121951862538D-02 )*Y-2.45653725061323D-01 )*Y+
     +     1.89999883453047D+00
      WW1 = (((((((((( 1.74841995087592D-15*Y-6.95671892641256D-16)*Y-
     +     3.000659497257D-13)*Y+2.021279817961D-13)*Y+
     +     3.853596935400D-11)*Y+1.461418533652D-10)*Y-
     +     1.014517563435D-08)*Y+1.132736008979D-07)*Y-
     +     2.86605475073259D-06 )*Y+1.21958354908768D-04 )*Y-
     +     3.86293751153466D-03 )*Y+1.45298342081522D-01
      WW2 = ((((((((((-1.11199320525573D-15*Y+1.85007587796671D-15)*Y+
     +     1.220613939709D-13)*Y+1.275068098526D-12)*Y-
     +     5.341838883262D-11)*Y+6.161037256669D-10)*Y-
     +     1.009147879750D-08)*Y+2.907862965346D-07)*Y-
     +     6.12300038720919D-06 )*Y+1.00104454489518D-04 )*Y-
     +     1.80677298502757D-03 )*Y+5.78009914536630D-02
      WW3 = ((((((((((-9.49816486853687D-16*Y+6.67922080354234D-15)*Y+
     +     2.606163540537D-15)*Y+1.983799950150D-12)*Y-
     +     5.400548574357D-11)*Y+6.638043374114D-10)*Y-
     +     8.799518866802D-09)*Y+1.791418482685D-07)*Y-
     +     2.96075397351101D-06 )*Y+3.38028206156144D-05 )*Y-
     +     3.58426847857878D-04 )*Y+8.39213709428516D-03
      WW4 = ((((((((((( 1.33829971060180D-17*Y-3.44841877844140D-16)*Y+
     +     4.745009557656D-15)*Y-6.033814209875D-14)*Y+
     +     1.049256040808D-12)*Y-1.70859789556117D-11 )*Y+
     +     2.15219425727959D-10 )*Y-2.52746574206884D-09 )*Y+
     +     3.27761714422960D-08 )*Y-3.90387662925193D-07 )*Y+
     +     3.46340204593870D-06 )*Y-2.43236345136782D-05 )*Y+
     +     3.54846978585226D-04
      WW5 = ((((((((((((( 2.69412277020887D-20*Y-4.24837886165685D-19)*
     +     Y+6.030500065438D-18)*Y-9.069722758289D-17)*Y+
     +     1.246599177672D-15)*Y-1.56872999797549D-14 )*Y+
     +     1.87305099552692D-13 )*Y-2.09498886675861D-12 )*Y+
     +     2.11630022068394D-11 )*Y-1.92566242323525D-10 )*Y+
     +     1.62012436344069D-09 )*Y-1.23621614171556D-08 )*Y+
     +     7.72165684563049D-08 )*Y-3.59858901591047D-07 )*Y+
     +     2.43682618601000D-06
      RETURN
C
C     X=20.0 TO 25.0                             NROOTS = 5
  200 Y = X-22.5D+00
      RT1 = (((((((((-1.13927848238726D-15*Y+7.39404133595713D-15)*Y+
     +     1.445982921243D-13)*Y-2.676703245252D-12)*Y+
     +     5.823521627177D-12)*Y+2.17264723874381D-10 )*Y+
     +     3.56242145897468D-09 )*Y-3.03763737404491D-07 )*Y+
     +     9.46859114120901D-06 )*Y-2.30896753853196D-04 )*Y+
     +     5.24663913001114D-03
      RT2 = (((((((((( 2.89872355524581D-16*Y-1.22296292045864D-14)*Y+
     +     6.184065097200D-14)*Y+1.649846591230D-12)*Y-
     +     2.729713905266D-11)*Y+3.709913790650D-11)*Y+
     +     2.216486288382D-09)*Y+4.616160236414D-08)*Y-
     +     3.32380270861364D-06 )*Y+9.84635072633776D-05 )*Y-
     +     2.30092118015697D-03 )*Y+5.00845183695073D-02
      RT3 = (((((((((( 1.97068646590923D-15*Y-4.89419270626800D-14)*Y+
     +     1.136466605916D-13)*Y+7.546203883874D-12)*Y-
     +     9.635646767455D-11)*Y-8.295965491209D-11)*Y+
     +     7.534109114453D-09)*Y+2.699970652707D-07)*Y-
     +     1.42982334217081D-05 )*Y+3.78290946669264D-04 )*Y-
     +     8.03133015084373D-03 )*Y+1.58689469640791D-01
      RT4 = (((((((((( 1.33642069941389D-14*Y-1.55850612605745D-13)*Y-
     +     7.522712577474D-13)*Y+3.209520801187D-11)*Y-
     +     2.075594313618D-10)*Y-2.070575894402D-09)*Y+
     +     7.323046997451D-09)*Y+1.851491550417D-06)*Y-
     +     6.37524802411383D-05 )*Y+1.36795464918785D-03 )*Y-
     +     2.42051126993146D-02 )*Y+3.97847167557815D-01
      RT5 = ((((((((((-6.07053986130526D-14*Y+1.04447493138843D-12)*Y-
     +     4.286617818951D-13)*Y-2.632066100073D-10)*Y+
     +     4.804518986559D-09)*Y-1.835675889421D-08)*Y-
     +     1.068175391334D-06)*Y+3.292234974141D-05)*Y-
     +     5.94805357558251D-04 )*Y+8.29382168612791D-03 )*Y-
     +     9.93122509049447D-02 )*Y+1.09857804755042D+00
      WW1 = (((((((((-9.10338640266542D-15*Y+1.00438927627833D-13)*Y+
     +     7.817349237071D-13)*Y-2.547619474232D-11)*Y+
     +     1.479321506529D-10)*Y+1.52314028857627D-09 )*Y+
     +     9.20072040917242D-09 )*Y-2.19427111221848D-06 )*Y+
     +     8.65797782880311D-05 )*Y-2.82718629312875D-03 )*Y+
     +     1.28718310443295D-01
      WW2 = ((((((((( 5.52380927618760D-15*Y-6.43424400204124D-14)*Y-
     +     2.358734508092D-13)*Y+8.261326648131D-12)*Y+
     +     9.229645304956D-11)*Y-5.68108973828949D-09 )*Y+
     +     1.22477891136278D-07 )*Y-2.11919643127927D-06 )*Y+
     +     4.23605032368922D-05 )*Y-1.14423444576221D-03 )*Y+
     +     5.06607252890186D-02
      WW3 = ((((((((( 3.99457454087556D-15*Y-5.11826702824182D-14)*Y-
     +     4.157593182747D-14)*Y+4.214670817758D-12)*Y+
     +     6.705582751532D-11)*Y-3.36086411698418D-09 )*Y+
     +     6.07453633298986D-08 )*Y-7.40736211041247D-07 )*Y+
     +     8.84176371665149D-06 )*Y-1.72559275066834D-04 )*Y+
     +     7.16639814253567D-03
      WW4 = (((((((((((-2.14649508112234D-18*Y-2.45525846412281D-18)*Y+
     +     6.126212599772D-16)*Y-8.526651626939D-15)*Y+
     +     4.826636065733D-14)*Y-3.39554163649740D-13 )*Y+
     +     1.67070784862985D-11 )*Y-4.42671979311163D-10 )*Y+
     +     6.77368055908400D-09 )*Y-7.03520999708859D-08 )*Y+
     +     6.04993294708874D-07 )*Y-7.80555094280483D-06 )*Y+
     +     2.85954806605017D-04
      WW5 = ((((((((((((-5.63938733073804D-21*Y+6.92182516324628D-20)*Y-
     +     1.586937691507D-18)*Y+3.357639744582D-17)*Y-
     +     4.810285046442D-16)*Y+5.386312669975D-15)*Y-
     +     6.117895297439D-14)*Y+8.441808227634D-13)*Y-
     +     1.18527596836592D-11 )*Y+1.36296870441445D-10 )*Y-
     +     1.17842611094141D-09 )*Y+7.80430641995926D-09 )*Y-
     +     5.97767417400540D-08 )*Y+1.65186146094969D-06
      RETURN
C
  220 WW1 = SQRT(PIE4/X)
      IF (X .GT. 40.0D+00) GO TO 240
C     X=25.0 TO 40.0                             NROOTS = 5
      E = EXP(-X)
      RT1 = ((((((((-1.73363958895356D-06*X+1.19921331441483D-04)*X -
     +     1.59437614121125D-02)*X+1.13467897349442D+00)*X -
     +     4.47216460864586D+01)*X+1.06251216612604D+03)*X -
     +     1.52073917378512D+04)*X+1.20662887111273D+05)*X -
     +     4.07186366852475D+05)*E + R15/(X-R15)
      RT2 = ((((((((-1.60102542621710D-05*X+1.10331262112395D-03)*X -
     +     1.50043662589017D-01)*X+1.05563640866077D+01)*X -
     +     4.10468817024806D+02)*X+9.62604416506819D+03)*X -
     +     1.35888069838270D+05)*X+1.06107577038340D+06)*X -
     +     3.51190792816119D+06)*E + R25/(X-R25)
      RT3 = ((((((((-4.48880032128422D-05*X+2.69025112122177D-03)*X -
     +     4.01048115525954D-01)*X+2.78360021977405D+01)*X -
     +     1.04891729356965D+03)*X+2.36985942687423D+04)*X -
     +     3.19504627257548D+05)*X+2.34879693563358D+06)*X -
     +     7.16341568174085D+06)*E + R35/(X-R35)
      RT4 = ((((((((-6.38526371092582D-05*X-2.29263585792626D-03)*X -
     +     7.65735935499627D-02)*X+9.12692349152792D+00)*X -
     +     2.32077034386717D+02)*X+2.81839578728845D+02)*X +
     +     9.59529683876419D+04)*X-1.77638956809518D+06)*X +
     +     1.02489759645410D+07)*E + R45/(X-R45)
      RT5 = ((((((((-3.59049364231569D-05*X-2.25963977930044D-02)*X +
     +     1.12594870794668D+00)*X-4.56752462103909D+01)*X +
     +     1.05804526830637D+03)*X-1.16003199605875D+04)*X -
     +     4.07297627297272D+04)*X+2.22215528319857D+06)*X -
     +     1.61196455032613D+07)*E + R55/(X-R55)
      WW5 = (((((((((-4.61100906133970D-10*X+1.43069932644286D-07)*X -
     +     1.63960915431080D-05)*X+1.15791154612838D-03)*X -
     +     5.30573476742071D-02)*X+1.61156533367153D+00)*X -
     +     3.23248143316007D+01)*X+4.12007318109157D+02)*X -
     +     3.02260070158372D+03)*X+9.71575094154768D+03)*E + W55*WW1
      WW4 = (((((((((-2.40799435809950D-08*X+8.12621667601546D-06)*X -
     +     9.04491430884113D-04)*X+6.37686375770059D-02)*X -
     +     2.96135703135647D+00)*X+9.15142356996330D+01)*X -
     +     1.86971865249111D+03)*X+2.42945528916947D+04)*X -
     +     1.81852473229081D+05)*X+5.96854758661427D+05)*E + W45*WW1
      WW3 = (((((((( 1.83574464457207D-05*X-1.54837969489927D-03)*X +
     +     1.18520453711586D-01)*X-6.69649981309161D+00)*X +
     +     2.44789386487321D+02)*X-5.68832664556359D+03)*X +
     +     8.14507604229357D+04)*X-6.55181056671474D+05)*X +
     +     2.26410896607237D+06)*E + W35*WW1
      WW2 = (((((((( 2.77778345870650D-05*X-2.22835017655890D-03)*X +
     +     1.61077633475573D-01)*X-8.96743743396132D+00)*X +
     +     3.28062687293374D+02)*X-7.65722701219557D+03)*X +
     +     1.10255055017664D+05)*X-8.92528122219324D+05)*X +
     +     3.10638627744347D+06)*E + W25*WW1
      WW1 = WW1-0.01962D+00*E-WW2-WW3-WW4-WW5
      RETURN
C
  240 IF (X .GT. 59.0D+00) GO TO 260
C     X=40.0 TO 59.0                             NROOTS = 5
      XXX = X**3
      E = XXX*EXP(-X)
      RT1 = (((-2.43758528330205D-02*X+2.07301567989771D+00)*X -
     +     6.45964225381113D+01)*X+7.14160088655470D+02)*E + R15/(X-R15)
      RT2 = (((-2.28861955413636D-01*X+1.93190784733691D+01)*X -
     +     5.99774730340912D+02)*X+6.61844165304871D+03)*E + R25/(X-R25)
      RT3 = (((-6.95053039285586D-01*X+5.76874090316016D+01)*X -
     +     1.77704143225520D+03)*X+1.95366082947811D+04)*E + R35/(X-R35)
      RT4 = (((-1.58072809087018D+00*X+1.27050801091948D+02)*X -
     +     3.86687350914280D+03)*X+4.23024828121420D+04)*E + R45/(X-R45)
      RT5 = (((-3.33963830405396D+00*X+2.51830424600204D+02)*X -
     +     7.57728527654961D+03)*X+8.21966816595690D+04)*E + R55/(X-R55)
      E = XXX*E
      WW5 = (( 1.35482430510942D-08*X-3.27722199212781D-07)*X +
     +     2.41522703684296D-06)*E + W55*WW1
      WW4 = (( 1.23464092261605D-06*X-3.55224564275590D-05)*X +
     +     3.03274662192286D-04)*E + W45*WW1
      WW3 = (( 1.34547929260279D-05*X-4.19389884772726D-04)*X +
     +     3.87706687610809D-03)*E + W35*WW1
      WW2 = (( 2.09539509123135D-05*X-6.87646614786982D-04)*X +
     +     6.68743788585688D-03)*E + W25*WW1
      WW1 = WW1-WW2-WW3-WW4-WW5
      RETURN
C
C     X=59.0 TO INFINITY                         NROOTS = 5
  260 RT1 = R15/(X-R15)
      RT2 = R25/(X-R25)
      RT3 = R35/(X-R35)
      RT4 = R45/(X-R45)
      RT5 = R55/(X-R55)
      WW2 = W25*WW1
      WW3 = W35*WW1
      WW4 = W45*WW1
      WW5 = W55*WW1
      WW1 = WW1-WW2-WW3-WW4-WW5
      RETURN
      END
C
C*MODULE RYSPOL  *DECK ROOT6
      SUBROUTINE ROOT6
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MR=13, MXAUX=55)
C
      DIMENSION RGRID(MXAUX),WGRID(MXAUX),
     *          P0(MXAUX),P1(MXAUX),P2(MXAUX),
     *          RTS(MR),WTS(MR),WRK(MR),
     *          ALPHA(0:MR-1),BETA(0:MR-1)
C
      COMMON /ROOT  / XX,UF(MR),WF(MR),NROOTS
      COMMON /RYSPAR/ XASYMP(MR),RTSASY(MR,MR),WTSASY(MR,MR),
     *                NAUXS(MR),MAPRYS(MR),RTSAUX(55,8),WTSAUX(55,8)
C
      PARAMETER (ONE=1.0D+00)
C
C           GENERAL CASE (NROOTS=1 TO 13) OF RYS ROOTS/WEIGHTS
C
C           PROGRAMMED NOVEMBER 2005 BY MIKE SCHMIDT TO REPLACE
C           THE VENERABLE HONDO CODE DUE TO MICHEL DUPUIS, USED
C           FOR THREE DECADES.
C           SEE KING/DUPUIS  J.COMPUT.PHYS. 21, 144-165(1976)
C
      EPS=1.0D-14
C
C           AYMPTOTICALLY LARGE -XX- IS SIMPLE
C           SEE ISHIDA FOR THE DETAILS: JCP 98, 2176-2181(1993)
C
      IF(XX.GE.XASYMP(NROOTS)) THEN
         CALL RYSASY(NROOTS,XX,RTS,WTS)
C
C           GAUTSCHI'S DISCRETIZED STIELTJES "BOOTSTRAPPING", WHICH
C           WORKS WITH NUMERICAL VALUES OF THE RYS POLYNOMIALS ON
C           THE GRID, YIELDING NUMERICALLY STABLE ALPHA,BETA VALUES.
C           SEE SAGAR/SMITH FOR THE DETAILS: IJQC 42, 827-836(1992)
C
      ELSE
         NAUX=NAUXS(NROOTS)
         MAP=MAPRYS(NROOTS)
         DO I=1,NAUX
            T2 = RTSAUX(I,MAP)*RTSAUX(I,MAP)
            RGRID(I) = T2
            WGRID(I) = WTSAUX(I,MAP)*EXP(-XX*T2)
         ENDDO
         CALL RYSDS(NROOTS,NAUX,RGRID,WGRID,ALPHA,BETA,IERR,P0,P1,P2)
         CALL RYSGW(NROOTS,ALPHA,BETA,EPS,RTS,WTS,IERR,WRK)
      END IF
C
      DO K=1,NROOTS
         DUM  =RTS(K)
         UF(K)=DUM/(ONE-DUM)
         WF(K)=WTS(K)
      ENDDO
      RETURN
      END
C
C*MODULE RYSPOL  *DECK RYSASY
      SUBROUTINE RYSASY(NROOTS,XX,RTS,WTS)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RTS(NROOTS),WTS(NROOTS)
      PARAMETER (MR=13)
      COMMON /RYSPAR/ XASYMP(MR),RTSASY(MR,MR),WTSASY(MR,MR),
     *                NAUXS(MR),MAPRYS(MR),RTSAUX(55,8),WTSAUX(55,8)
      PARAMETER (ONE=1.0D+00)
C
C     THE ROOTS AND WEIGHTS FOR XX = BIG ARE GIVEN BY:
C       T*T = S*S/XX    W = V/SQRT(XX)
C     WHERE S AND V ARE THE ROOTS AND WEIGHTS OF THE
C     HERMITE POLYNOMIALS, OF ORDER 2*NROOTS.
C     SEE ISHIDA-SAN'S SECOND PAPER FOR DETAILS, ALTHOUGH HE APPEARS
C     TO HAVE BEEN A BIT TOO CONSERVATIVE ABOUT WHEN TO GO ASYMPTOTIC.
C
      FACTR = ONE/XX
      FACTW = SQRT(FACTR)
      DO I=1,NROOTS
         RTS(I)= FACTR * RTSASY(I,NROOTS)
         WTS(I)= FACTW * WTSASY(I,NROOTS)
      ENDDO
      RETURN
      END
C
C*MODULE RYSPOL  *DECK RYSGW
      SUBROUTINE RYSGW(N,ALPHA,BETA,EPS,ROOTS,WEIGHT,IERR,WRK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION ALPHA(N),BETA(N),ROOTS(N),WEIGHT(N),WRK(N)
C
C     THIS ROUTINE IMPLEMENTS THE QL PROCEDURE DUE TO GOLUB/WELSCH,
C             "CALCULATION OF GAUSS QUADRATURE RULES"
C                 GENE H. GOLUB, JOHN H. WELSCH
C           MATHEMATICS OF COMPUTATION, 23, 221-230(1969)
C
C     THIS ROUTINE (ORIGINAL NAME=GAUSS) IS A SUPPLEMENT TO THE ARTICLE,
C      ALGORITHM 726: ORTHPOL - A PACKAGE OF ROUTINES FOR GENERATING
C        ORTHOGONAL POLYNOMIALS AND GAUSS-TYPE QUADRATURE RULES.
C                         WALTER GAUTSCHI,
C      ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, 20, 21-62(1994)
C
C            A FOOTNOTE IN GAUTSCHI'S ARTICLE REVEALS THAT
C            THIS PARTICULAR ROUTINE CAME FROM GENE GOLUB.
C
C GIVEN  N  AND A MEASURE  DLAMBDA, THIS ROUTINE GENERATES THE N-POINT
C GAUSSIAN QUADRATURE FORMULA
C
C     INTEGRAL OVER SUPP(DLAMBDA) OF F(X)DLAMBDA(X)
C
C        = SUM FROM K=1 TO K=N OF W(K)F(X(K)) + R(N;F).
C
C THE NODES ARE RETURNED AS  ROOS(K)=X(K) AND THE WEIGHTS AS
C WEIGHT(K)=W(K), K=1,2,...,N. THE USER HAS TO SUPPLY THE RECURSION
C COEFFICIENTS  ALPHA(K), BETA(K), K=0,1,2,...,N-1, FOR THE MEASURE
C DLAMBDA. THE ROUTINE COMPUTES THE NODES AS EIGENVALUES, AND THE
C WEIGHTS IN TERM OF THE FIRST COMPONENT OF THE RESPECTIVE NORMALIZED
C EIGENVECTORS OF THE N-TH ORDER JACOBI MATRIX ASSOCIATED WITH  DLAMBDA.
C IT USES A TRANSLATION AND ADAPTATION OF THE ALGOL PROCEDURE  IMTQL2,
C NUMER. MATH. 12, 1968, 377-383, BY MARTIN AND WILKINSON, AS MODIFIED
C BY DUBRULLE, NUMER. MATH. 15, 1970, 450. SEE ALSO HANDBOOK FOR
C AUTOM. COMPUT., VOL. 2 - LINEAR ALGEBRA, PP.241-248, AND THE EISPACK
C ROUTINE  IMTQL2.
C
C        INPUT:  N - - THE NUMBER OF POINTS IN THE GAUSSIAN QUADRATURE
C                      FORMULA; TYPE INTEGER
C                ALPHA,BETA - - ARRAYS OF DIMENSION  N  TO BE FILLED
C                      WITH THE VALUES OF  ALPHA(K-1), BETA(K-1), K=1,2,
C                      ...,N
C                EPS - THE RELATIVE ACCURACY DESIRED IN THE NODES
C                      AND WEIGHTS
C
C        OUTPUT: ROOTS- ARRAY OF DIMENSION  N  CONTAINING THE GAUSSIAN
C                      NODES (IN INCREASING ORDER)  ROOTS(K)=X(K),
C                       K=1,2,...,N
C                WEIGHT - ARRAY OF DIMENSION  N  CONTAINING THE
C                      GAUSSIAN WEIGHTS  WEIGHT(K)=W(K), K=1,2,...,N
C                IERR- AN ERROR FLAG EQUAL TO  0  ON NORMAL RETURN,
C                      EQUAL TO  I  IF THE QR ALGORITHM DOES NOT
C                      CONVERGE WITHIN 30 ITERATIONS ON EVALUATING THE
C                      I-TH EIGENVALUE, EQUAL TO  -1  IF  N  IS NOT IN
C                      RANGE, AND EQUAL TO  -2  IF ONE OF THE BETA'S IS
C                      NEGATIVE.
C
C THE ARRAY  WRK  IS NEEDED FOR WORKING SPACE.
C
      IF(N.LT.1) THEN
        IERR=-1
        RETURN
      END IF
C
      IERR=0
      ROOTS(1)=ALPHA(1)
      IF(BETA(1).LT.0.0D+00) THEN
        IERR=-2
        RETURN
      END IF
      WEIGHT(1)=BETA(1)
      IF (N.EQ.1) RETURN
C
      WEIGHT(1)=1.0D+00
      WRK(N)=0.0D+00
      DO 100 K=2,N
        ROOTS(K)=ALPHA(K)
        IF(BETA(K).LT.0.0D+00) THEN
          IERR=-2
          RETURN
        END IF
        WRK(K-1)=SQRT(BETA(K))
        WEIGHT(K)=0.0D+00
  100 CONTINUE
C
      DO 240 L=1,N
        J=0
C
C LOOK FOR A SMALL SUBDIAGONAL ELEMENT.
C
  105   DO 110 M=L,N
          IF(M.EQ.N) GO TO 120
          IF(ABS(WRK(M)).LE.EPS*(ABS(ROOTS(M))+ABS(ROOTS(M+1))))
     *      GO TO 120
  110   CONTINUE
  120   DP=ROOTS(L)
        IF(M.EQ.L) GO TO 240
        IF(J.EQ.30) GO TO 400
        J=J+1
C
C FORM SHIFT.
C
        DG=(ROOTS(L+1)-DP)/(2.0D+00*WRK(L))
        DR=SQRT(DG*DG+1.0D+00)
        DG=ROOTS(M)-DP+WRK(L)/(DG+SIGN(DR,DG))
        DS=1.0D+00
        DC=1.0D+00
        DP=0.0D+00
        MML=M-L
C
C FOR I=M-1 STEP -1 UNTIL L DO ...
C
        DO 200 II=1,MML
          I=M-II
          DF=DS*WRK(I)
          DB=DC*WRK(I)
          IF(ABS(DF).LT.ABS(DG)) GO TO 150
          DC=DG/DF
          DR=SQRT(DC*DC+1.0D+00)
          WRK(I+1)=DF*DR
          DS=1.0D+00/DR
          DC=DC*DS
          GO TO 160
  150     DS=DF/DG
          DR=SQRT(DS*DS+1.0D+00)
          WRK(I+1)=DG*DR
          DC=1.0D+00/DR
          DS=DS*DC
  160     DG=ROOTS(I+1)-DP
          DR=(ROOTS(I)-DG)*DS+2.0D+00*DC*DB
          DP=DS*DR
          ROOTS(I+1)=DG+DP
          DG=DC*DR-DB
C
C FORM FIRST COMPONENT OF VECTOR.
C
          DF=WEIGHT(I+1)
          WEIGHT(I+1)=DS*WEIGHT(I)+DC*DF
          WEIGHT(I)=DC*WEIGHT(I)-DS*DF
  200   CONTINUE
        ROOTS(L)=ROOTS(L)-DP
        WRK(L)=DG
        WRK(M)=0.0D+00
        GO TO 105
  240 CONTINUE
C
C ORDER EIGENVALUES AND EIGENVECTORS.
C
      DO 300 II=2,N
        I=II-1
        K=I
        DP=ROOTS(I)
        DO 260 J=II,N
          IF(ROOTS(J).GE.DP) GO TO 260
          K=J
          DP=ROOTS(J)
  260   CONTINUE
        IF(K.EQ.I) GO TO 300
        ROOTS(K)=ROOTS(I)
        ROOTS(I)=DP
        DP=WEIGHT(I)
        WEIGHT(I)=WEIGHT(K)
        WEIGHT(K)=DP
  300 CONTINUE
      DO 310 K=1,N
        WEIGHT(K)=BETA(1)*WEIGHT(K)*WEIGHT(K)
  310 CONTINUE
      RETURN
C
C SET ERROR - NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS.
C
  400 IERR=L
      RETURN
      END
C
C*MODULE RYSPOL  *DECK RYSDS
      SUBROUTINE RYSDS(N,NCAP,X,W,ALPHA,BETA,IERR,P0,P1,P2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(NCAP),W(NCAP),ALPHA(N),BETA(N),
     *          P0(NCAP),P1(NCAP),P2(NCAP)
C
C     THIS ROUTINE (ORIGINAL NAME=STI) IS A SUPPLEMENT TO THE ARTICLE,
C      ALGORITHM 726: ORTHPOL - A PACKAGE OF ROUTINES FOR GENERATING
C        ORTHOGONAL POLYNOMIALS AND GAUSS-TYPE QUADRATURE RULES.
C                         WALTER GAUTSCHI,
C      ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, 20, 21-62(1994)
C
C THIS ROUTINE APPLIES STIELTJES'S PROCEDURE (CF. SECTION 2.1 OF
C W. GAUTSCHI, 'ON GENERATING ORTHOGONAL POLYNOMIALS', SIAM J. SCI.
C STATIST. COMPUT. 3, 1982, 289-317) TO GENERATE THE RECURSION
C COEFFICIENTS  ALPHA(K), BETA(K) , K=0,1,...,N-1, FOR THE DISCRETE
C (MONIC) ORTHOGONAL POLYNOMIALS ASSOCIATED WITH THE INNER PRODUCT
C
C     (F,G)=SUM OVER K FROM 1 TO NCAP OF W(K)*F(X(K))*G(X(K)).
C
C THE INTEGER  N  MUST BE BETWEEN  1  AND  NCAP, INCLUSIVE; OTHERWISE,
C THERE IS AN ERROR EXIT WITH  IERR=1. THE RESULTS ARE STORED IN THE
C ARRAYS  ALPHA, BETA; THE ARRAYS  P0, P1, P2  ARE WORKING ARRAYS.
C
C IF THERE IS A THREAT OF UNDERFLOW OR OVERFLOW IN THE CALCULATION
C OF THE COEFFICIENTS  ALPHA(K)  AND  BETA(K), THE ROUTINE EXITS WITH
C THE ERROR FLAG  IERR  SET EQUAL TO  -K  (IN THE CASE OF UNDERFLOW)
C OR  +K  (IN THE CASE OF OVERFLOW), WHERE  K  IS THE RECURSION INDEX
C FOR WHICH THE PROBLEM OCCURS. THE FORMER [LATTER] CAN OFTEN BE AVOIDED
C BY MULTIPLYING ALL WEIGHTS  W(K)  BY A SUFFICIENTLY LARGE [SMALL]
C SCALING FACTOR PRIOR TO ENTERING THE ROUTINE, AND, UPON EXIT, DIVIDE
C THE COEFFICIENT  BETA(0)  BY THE SAME FACTOR.
C
C THIS ROUTINE SHOULD BE USED WITH CAUTION IF  N  IS RELATIVELY CLOSE
C TO  NCAP, SINCE THERE IS A DISTINCT POSSIBILITY OF NUMERICAL
C INSTABILITY DEVELOPING. (SEE W. GAUTSCHI,IS THE RECURRENCE RELATION
C FOR ORTHOGONAL POLYNOMIALS ALWAYS STABLE?'', BIT, 1993, TO APPEAR.)
C IN THAT CASE, THE ROUTINE  LANCZ  SHOULD BE USED.
C
      TINY = 1.0D-40
      HUGE = 1.0D+40
C
      IERR=0
      IF(N.LE.0 .OR. N.GT.NCAP) THEN
        IERR=1
        RETURN
      END IF
      NM1=N-1
C
C COMPUTE THE FIRST ALPHA- AND BETA-COEFFICIENT.
C
      SUM0=0.0D+00
      SUM1=0.0D+00
      DO 10 M=1,NCAP
        SUM0=SUM0+W(M)
        SUM1=SUM1+W(M)*X(M)
   10 CONTINUE
      ALPHA(1)=SUM1/SUM0
      BETA(1)=SUM0
      IF(N.EQ.1) RETURN
C
C COMPUTE THE REMAINING ALPHA- AND BETA-COEFFICIENTS.
C
      DO 20 M=1,NCAP
        P1(M)=0.0D+00
        P2(M)=1.0D+00
   20 CONTINUE
      DO 40 K=1,NM1
        SUM1=0.0D+00
        SUM2=0.0D+00
        DO 30 M=1,NCAP
C
C THE FOLLOWING STATEMENT IS DESIGNED TO AVOID AN OVERFLOW CONDITION
C IN THE COMPUTATION OF  P2(M)  WHEN THE WEIGHTS  W(M)  GO TO ZERO
C FASTER (AND UNDERFLOW) THAN THE  P2(M)  GROW.
C
          IF(W(M).EQ.0.0D+00) GOTO 30
          P0(M)=P1(M)
          P1(M)=P2(M)
          P2(M)=(X(M)-ALPHA(K))*P1(M)-BETA(K)*P0(M)
C
C CHECK FOR IMPENDING OVERFLOW.
C
          IF(ABS(P2(M)).GT.HUGE .OR. ABS(SUM2).GT.HUGE) THEN
            IERR=K
            RETURN
          END IF
          T=W(M)*P2(M)*P2(M)
          SUM1=SUM1+T
          SUM2=SUM2+T*X(M)
   30   CONTINUE
C
C CHECK FOR IMPENDING UNDERFLOW.
C
        IF(ABS(SUM1).LT.TINY) THEN
          IERR=-K
          RETURN
        END IF
        ALPHA(K+1)=SUM2/SUM1
        BETA(K+1)=SUM1/SUM0
        SUM0=SUM1
   40 CONTINUE
C
      RETURN
      END
C
C*MODULE RYSPOL  *DECK SETRYS
      SUBROUTINE SETRYS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MR=13, MAUX=55)
      DIMENSION RTS(MAUX),WTS(MAUX),WRK(MAUX),
     *          ALPHA(0:MAUX-1),BETA(0:MAUX-1)
      COMMON /RYSPAR/ XASYMP(MR),RTSASY(MR,MR),WTSASY(MR,MR),
     *                NAUXS(MR),MAPRYS(MR),RTSAUX(55,8),WTSAUX(55,8)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00)
C
C     ----- INITIALIZE THE RYS QUADRATURE PROCEDURE -----
C
C     NUMERICAL PARAMETERS WERE DETERMINED BY EXPERIMENTATION:
C     A) THE -XX- VALUE AT WHICH WE CAN GO TO THE ASYMPTOTIC FORMULA
C        WAS FOUND BY MATCHING ROUGHLY 13 DIGITS IN ROOTS/WEIGHTS.
C     B) DENSITY OF THE AUXILIARY GRID FOR THE DISCRETIZATION WAS
C        DETERMINED BY EYEBALL, AT THE -XX- VALUE WHERE WE CAN BEGIN
C        TO USE THE ASYMPTOTIC FORMULA.  FOR A GIVEN -NROOTS-, THE
C        REQUIRED GRID DENSITY GRADUALLY INCREASES WITH -XX-.
C
      QUART = ONE/FOUR
      PI    = 3.141592653589793238D+00
      EPS   = 1.0D-14
C
C     SET UP ASYMPTOTIC ROOT COMPUTATION, BY GENERATING
C     THE ROOTS OF THE HERMITE POLYNOMIALS OF ORDER 2N.
C
      XASYMP( 1)=29.0D+00
      XASYMP( 2)=37.0D+00
      XASYMP( 3)=43.0D+00
      XASYMP( 4)=49.0D+00
      XASYMP( 5)=55.0D+00
      XASYMP( 6)=60.0D+00
      XASYMP( 7)=65.0D+00
      XASYMP( 8)=71.0D+00
      XASYMP( 9)=76.0D+00
      XASYMP(10)=81.0D+00
      XASYMP(11)=86.0D+00
      XASYMP(12)=91.0D+00
      XASYMP(13)=96.0D+00
C
      DO I=1,MR
C            NOTE THAT MAUX MUST BE AT LEAST TWICE MR, FROM NEXT LINE
         N=2*I
         DO J=0,N-1
            ALPHA(J) = ZERO
         ENDDO
         BETA(0)=SQRT(PI)
         DO J=1,N-1
            BETA(J) = J/TWO
         END DO
         CALL RYSGW(N,ALPHA,BETA,EPS,RTS,WTS,IERR,WRK)
         DO J=1,I
            RTSASY(J,I) = RTS(I+J)*RTS(I+J)
            WTSASY(J,I) = WTS(I+J)
         ENDDO
      ENDDO
C
C        GENERATE AUXILIARY GRIDS, AT 8 DISTINCT POINT DENSITIES
C
      NAUXS( 1)=20
      NAUXS( 2)=25
      NAUXS( 3)=30
      NAUXS( 4)=30
      NAUXS( 5)=35
      NAUXS( 6)=40
      NAUXS( 7)=40
      NAUXS( 8)=40
      NAUXS( 9)=45
      NAUXS(10)=50
      NAUXS(11)=50
      NAUXS(12)=55
      NAUXS(13)=55
C
      MAPRYS( 1)=1
      MAPRYS( 2)=2
      MAPRYS( 3)=3
      MAPRYS( 4)=3
      MAPRYS( 5)=4
      MAPRYS( 6)=5
      MAPRYS( 7)=5
      MAPRYS( 8)=5
      MAPRYS( 9)=6
      MAPRYS(10)=7
      MAPRYS(11)=7
      MAPRYS(12)=8
      MAPRYS(13)=8
C
C        THE AUXILIARY QUADRATURE IS TAKEN TO BE "SHIFTED LEGENDRE"
C
      NAUXSV=0
      IGRID=0
      DO 100 M=1,MR
         NAUX=NAUXS(M)
         IF(NAUX.EQ.NAUXSV) GO TO 100
         IGRID=IGRID+1
         NAUXSV=NAUX
         DO I=0,NAUX-1
            ALPHA(I) = ONE/TWO
         ENDDO
         BETA(0)= ONE
         DO I=1,NAUX-1
            BETA(I) = QUART/(FOUR-(ONE/(I*I)))
         END DO
         CALL RYSGW(NAUX,ALPHA,BETA,EPS,RTS,WTS,IERR,WRK)
         DO I=1,NAUX
            RTSAUX(I,IGRID) = RTS(I)
            WTSAUX(I,IGRID) = WTS(I)
         ENDDO
  100 CONTINUE
      RETURN
      END
