LHAPDF/lhapdf5.2.2/wrapcteq6.f

   1       subroutine CTEQ6evolve(x,Q,pdf)
   2       implicit real*8(a-h,o-z)
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6       common/NAME/name,nmem,ndef,mmem
   7       real*8 pdf(-6:6)
   8       integer nset,k
   9       Character Line*80
  10       PARAMETER (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
  11       PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX)
  12       Common
  13      > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
  14      > / CtqPar2 / Nx, Nt, NfMx
  15      > / XQrange / Qini, Qmax, Xmin
  16      > / QCDtable /  Alambda, Nfl, Iorder
  17      > / Masstbl / Amass(6)
  18       common/masses_LHA/cMass,bMass,tMass
  19       data pi / 3.141592653589793d0 /
  20       save
  21 c
  22       call getnset(iset)
  23       call getnmem(iset,imem)
  24 c
  25       U =         X * CtLhCtq6Pdf(imem,1,X,Q)
  26       D =         X * CtLhCtq6Pdf(imem,2,X,Q)
  27       USEA =      X * CtLhCtq6Pdf(imem,-1,X,Q)
  28       DSEA =      X * CtLhCtq6Pdf(imem,-2,X,Q)
  29       STR =       X * CtLhCtq6Pdf(imem,3,X,Q)
  30       CHM =       X * CtLhCtq6Pdf(imem,4,X,Q)
  31       BOT =       X * CtLhCtq6Pdf(imem,5,X,Q)
  32       GLU  =      X * CtLhCtq6Pdf(imem,0,X,Q)
  33       UPV=U-USEA
  34       DNV=D-DSEA
  35 c
  36       pdf(0)  = glu
  37       pdf(1)  = dnv+dsea
  38       pdf(-1) = dsea
  39       pdf(2)  = upv+usea
  40       pdf(-2) = usea
  41       pdf(3)  = str
  42       pdf(-3) = str
  43       pdf(4)  = chm
  44       pdf(-4) = chm
  45       pdf(5)  = bot
  46       pdf(-5) = bot
  47       pdf(6)  = 0.0d0
  48       pdf(-6) = 0.0d0
  49
  50       return
  51 *
  52       entry CTEQ6read(nset)
  53
  54       call CtLhbldat1           !this line was missing in previous releases (jcp)
  55       call CtLhbldat2           !this line was missing in previous releases (jcp)
  56
  57       read(1,*)nmem(nset),ndef(nset)    !*** nmem+1=number of members; ndef is not used for anything ***
  58       if(nmem(nset) .gt. nhess) then
  59          print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
  60          stop
  61       endif
  62       Read  (1, '(A)') Line
  63       Read  (1, '(A)') Line
  64       Read  (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
  65       Iorder = Nint(Dr)
  66       Nfl = Nint(Fl)
  67       Alambda = Al
  68
  69       cMass = Amass(4)                  !missing in previous releases(jcp)
  70       bMass = Amass(5)                  !missing in previous releases(jcp)
  71       tMass = Amass(6)                  !missing in previous releases(jcp)
  72
  73       Read  (1, '(A)') Line
  74       Read  (1, *) NX,  NT, NfMx
  75
  76       Read  (1, '(A)') Line
  77       Read  (1, *) QINI, QMAX, (TV(I), I =0, NT)
  78
  79       Read  (1, '(A)') Line
  80       Read  (1, *) XMIN, (XV(I), I =0, NX)
  81
  82       Do 11 Iq = 0, NT
  83          TV(Iq) = Log(Log (TV(Iq) /Al))
  84    11 Continue
  85 C
  86 C                  Since quark = anti-quark for nfl>2 at this stage,
  87 C                  we Read  out only the non-redundent data points
  88 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
  89
  90       Nblk = (NX+1) * (NT+1)
  91       Npts =  Nblk  * (NfMx+3)
  92
  93
  94 c **  do ihess = 0,nhess        !*** old version ***
  95       do ihess = 0,nmem(nset)           !*** new version: allows nmem < nhess ***
  96
  97         Read  (1, '(A)') Line
  98         Read  (1, '(A)') Line
  99         Read  (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts)
 100
 101       enddo
 102       return
 103 *
 104       entry CTEQ6alfa(alfas,Qalfa)
 105       alfas = pi*CtLhALPI(Qalfa)
 106       return
 107 *
 108       entry CTEQ6init(Eorder,Q2fit)
 109
 110       return
 111 *
 112       entry CTEQ6pdf(mem)
 113 c        imem = mem
 114         call getnset(iset)
 115         call setnmem(iset,mem)
 116       return
 117 *
 118       end
 119
 120       SUBROUTINE CtLhPOLINT (XA,YA,N,X,Y,DY)
 121
 122       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 123 C                                        Adapted from "Numerical Recipes"
 124       PARAMETER (NMAX=10)
 125       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 126       NS=1
 127       DIF=ABS(X-XA(1))
 128       DO 11 I=1,N
 129         DIFT=ABS(X-XA(I))
 130         IF (DIFT.LT.DIF) THEN
 131           NS=I
 132           DIF=DIFT
 133         ENDIF
 134         C(I)=YA(I)
 135         D(I)=YA(I)
 136 11    CONTINUE
 137       Y=YA(NS)
 138       NS=NS-1
 139       DO 13 M=1,N-1
 140         DO 12 I=1,N-M
 141           HO=XA(I)-X
 142           HP=XA(I+M)-X
 143           W=C(I+1)-D(I)
 144           DEN=W/(HO-HP)
 145           D(I)=HP*DEN
 146           C(I)=HO*DEN
 147 12      CONTINUE
 148         IF (2*NS.LT.N-M)THEN
 149           DY=C(NS+1)
 150         ELSE
 151           DY=D(NS)
 152           NS=NS-1
 153         ENDIF
 154         Y=Y+DY
 155 13    CONTINUE
 156       RETURN
 157       END
 158
 159 c===========================================================================
 160       Function CtLhPartonX6 (iset,IPRTN, XX, QQ)
 161 c  Given the parton distribution function in the array U in
 162 c  COMMON / PEVLDT / , this routine interpolates to find
 163 c  the parton distribution at an arbitray point in x and q.
 164 c
 165       Implicit Double Precision (A-H,O-Z)
 166
 167       Parameter (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
 168       Parameter (MXQX= MXQ * MXX,   MXPQX = MXQX * (MXF+3))
 169
 170       Common
 171      > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
 172      > / CtqPar2 / Nx, Nt, NfMx
 173      > / XQrange / Qini, Qmax, Xmin
 174
 175         parameter(nqvec = 4)
 176
 177       Dimension fvec(4), fij(4)
 178       Dimension xvpow(0:mxx)
 179       Data OneP / 1.00001 /
 180       Data xpow / 0.3d0 /       !**** choice of interpolation variable
 181       Data nxsave / -999 /
 182       Save nxsave, xvpow
 183
 184         save xlast, qlast
 185         save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
 186         save const1 , const2, const3, const4, const5, const6
 187         save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
 188
 189 c store the powers used for interpolation on first call...
 190       if(nx .ne. nxsave) then
 191          nxsave = nx
 192          xvpow(0) = 0.D0
 193          do i = 1, nx
 194             xvpow(i) = xv(i)**xpow
 195          enddo
 196       endif
 197
 198       X = XX
 199       Q = QQ
 200
 201 c skip the initialization in x if same as in the previous call.
 202         if(x .eq. xlast) goto 100
 203         xlast = x
 204
 205 c      -------------    find lower end of interval containing x, i.e.,
 206 c                       get jx such that xv(jx) .le. x .le. xv(jx+1)...
 207       JLx = -1
 208       JU = Nx+1
 209  11   If (JU-JLx .GT. 1) Then
 210          JM = (JU+JLx) / 2
 211          If (X .Ge. XV(JM)) Then
 212             JLx = JM
 213          Else
 214             JU = JM
 215          Endif
 216          Goto 11
 217       Endif
 218 C                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
 219 C                           |---|---|---|...|---|-x-|---|...|---|---|
 220 C                     x     0  Xmin               x                 1
 221 C
 222       If     (JLx .LE. -1) Then
 223         Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX6 x=',x
 224         Stop
 225       ElseIf (JLx .Eq. 0) Then
 226          Jx = 0
 227       Elseif (JLx .LE. Nx-2) Then
 228
 229 C                For interior points, keep x in the middle, as shown above
 230          Jx = JLx - 1
 231       Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
 232
 233 C                  We tolerate a slight over-shoot of one (OneP=1.00001),
 234 C              perhaps due to roundoff or whatever, but not more than that.
 235 C                                      Keep at least 4 points >= Jx
 236          Jx = JLx - 2
 237       Else
 238         Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX6 x=',x
 239         Stop
 240       Endif
 241 C          ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
 242
 243 C                       This is the variable to be interpolated in
 244       ss = x**xpow
 245
 246       If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
 247
 248 c     initiation work for "interior bins": store the lattice points in s...
 249       svec1 = xvpow(jx)
 250       svec2 = xvpow(jx+1)
 251       svec3 = xvpow(jx+2)
 252       svec4 = xvpow(jx+3)
 253
 254       s12 = svec1 - svec2
 255       s13 = svec1 - svec3
 256       s23 = svec2 - svec3
 257       s24 = svec2 - svec4
 258       s34 = svec3 - svec4
 259
 260       sy2 = ss - svec2
 261       sy3 = ss - svec3
 262
 263 c constants needed for interpolating in s at fixed t lattice points...
 264       const1 = s13/s23
 265       const2 = s12/s23
 266       const3 = s34/s23
 267       const4 = s24/s23
 268       s1213 = s12 + s13
 269       s2434 = s24 + s34
 270       sdet = s12*s34 - s1213*s2434
 271       tmp = sy2*sy3/sdet
 272       const5 = (s34*sy2-s2434*sy3)*tmp/s12
 273       const6 = (s1213*sy2-s12*sy3)*tmp/s34
 274
 275       EndIf
 276
 277 100     continue
 278
 279 c skip the initialization in q if same as in the previous call.
 280         if(q .eq. qlast) goto 110
 281         qlast = q
 282
 283       tt = log(log(Q/Al))
 284
 285 c         --------------Now find lower end of interval containing Q, i.e.,
 286 c                          get jq such that qv(jq) .le. q .le. qv(jq+1)...
 287       JLq = -1
 288       JU = NT+1
 289  12   If (JU-JLq .GT. 1) Then
 290          JM = (JU+JLq) / 2
 291          If (tt .GE. TV(JM)) Then
 292             JLq = JM
 293          Else
 294             JU = JM
 295          Endif
 296          Goto 12
 297        Endif
 298
 299       If     (JLq .LE. 0) Then
 300          Jq = 0
 301       Elseif (JLq .LE. Nt-2) Then
 302 C                                  keep q in the middle, as shown above
 303          Jq = JLq - 1
 304       Else
 305 C                         JLq .GE. Nt-1 case:  Keep at least 4 points >= Jq.
 306         Jq = Nt - 3
 307
 308       Endif
 309 C                                   This is the interpolation variable in Q
 310
 311       If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
 312 c                                        store the lattice points in t...
 313       tvec1 = Tv(jq)
 314       tvec2 = Tv(jq+1)
 315       tvec3 = Tv(jq+2)
 316       tvec4 = Tv(jq+3)
 317
 318       t12 = tvec1 - tvec2
 319       t13 = tvec1 - tvec3
 320       t23 = tvec2 - tvec3
 321       t24 = tvec2 - tvec4
 322       t34 = tvec3 - tvec4
 323
 324       ty2 = tt - tvec2
 325       ty3 = tt - tvec3
 326
 327       tmp1 = t12 + t13
 328       tmp2 = t24 + t34
 329
 330       tdet = t12*t34 - tmp1*tmp2
 331
 332       EndIf
 333
 334 110     continue
 335
 336 c get the pdf function values at the lattice points...
 337
 338       If (Iprtn .GE. 3) Then
 339          Ip = - Iprtn
 340       Else
 341          Ip = Iprtn
 342       EndIf
 343       jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
 344
 345       Do it = 1, nqvec
 346
 347          J1  = jtmp + it*(NX+1)
 348
 349        If (Jx .Eq. 0) Then
 350 C                          For the first 4 x points, interpolate x^2*f(x,Q)
 351 C                           This applies to the two lowest bins JLx = 0, 1
 352 C            We cannot put the JLx.eq.1 bin into the "interior" section
 353 C                           (as we do for q), since Upd(J1) is undefined.
 354          fij(1) = 0
 355          fij(2) = Upd(iset,J1+1) * XV(1)**2
 356          fij(3) = Upd(iset,J1+2) * XV(2)**2
 357          fij(4) = Upd(iset,J1+3) * XV(3)**2
 358 C
 359 C                 Use CtLhPolint which allows x to be anywhere w.r.t. the grid
 360
 361          Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
 362
 363          If (x .GT. 0D0)  Fvec(it) =  Fx / x**2
 364 C                                              Pdf is undefined for x.eq.0
 365        ElseIf  (JLx .Eq. Nx-1) Then
 366 C                                                This is the highest x bin:
 367
 368 c** fix allow 4 consecutive elements with iset... mrw 19.9.2005
 369         fij(1) = Upd(iset,j1)
 370         fij(2) = Upd(iset,j1+1)
 371         fij(3) = Upd(iset,j1+2)
 372         fij(4) = Upd(iset,j1+3)
 373         Call CtLhPolint4 (XVpow(Nx-3), Fij(1), 4, ss, Fx, Dfx)
 374
 375         Fvec(it) = Fx
 376
 377        Else
 378
 379 C for all interior points, use Jon's in-line function
 380 C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
 381 c (This is cubic spline interpolation, as used by cteq; it was
 382 c changed to polint in previous Durham releases (jcp).)
 383          sf2 = Upd(iset,J1+1)
 384          sf3 = Upd(iset,J1+2)
 385
 386          Fvec(it) = (const5*(Upd(iset,J1)
 387      &                       - sf2*const1 + sf3*const2)
 388      &             + const6*(Upd(iset,J1+3)
 389      &                       + sf2*const3 - sf3*const4)
 390      &             + sf2*sy3 - sf3*sy2) / s23
 391
 392        Endif
 393
 394       enddo
 395 C                                   We now have the four values Fvec(1:4)
 396 c     interpolate in t...
 397
 398       If (JLq .LE. 0) Then
 399 C                         1st Q-bin, as well as extrapolation to lower Q
 400         Call CtLhPolint4(TV(0), Fvec(1), 4, tt, ff, Dfq)
 401
 402       ElseIf (JLq .GE. Nt-1) Then
 403 C                         Last Q-bin, as well as extrapolation to higher Q
 404         Call CtLhPolint4(TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
 405       Else
 406 C                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
 407 C       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
 408 C                         the full range QV(0:Nt)  (in contrast to XV)
 409         tf2 = fvec(2)
 410         tf3 = fvec(3)
 411
 412         g1 = ( tf2*t13 - tf3*t12) / t23
 413         g4 = (-tf2*t34 + tf3*t24) / t23
 414
 415         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
 416      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
 417
 418         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
 419       EndIf
 420
 421       CtLhPartonX6 = ff
 422
 423       Return
 424 C                                       ********************
 425       End
 426 c===========================================================================
 427       Function CtLhCtq6Pdf (iset,Iparton, X, Q)
 428       Implicit Double Precision (A-H,O-Z)
 429       Logical Warn
 430       Common
 431      > / CtqPar2 / Nx, Nt, NfMx
 432      > / QCDtable /  Alambda, Nfl, Iorder
 433
 434       Data Warn /.true./
 435       save Warn
 436
 437       If (X .lt. 0D0 .or. X .gt. 1D0) Then
 438         Print *, 'X out of range in CtLhCtq6Pdf: ', X
 439         Stop
 440       Endif
 441       If (Q .lt. Alambda) Then
 442         Print *, 'Q out of range in CtLhCtq6Pdf: ', Q
 443         Stop
 444       Endif
 445
 446 c added to force pdf = 0.0 at x=1.0 exactly - mrw
 447       if(x .eq. 1.0d0) then
 448           CtLhCtq6Pdf = 0.0d0
 449           return
 450       endif
 451 c
 452       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
 453          If (Warn) Then
 454 C        put a warning for calling extra flavor.
 455              Warn = .false.
 456              Print *, 'Warning: Iparton out of range in CtLhCtq6Pdf: '
 457      >              , Iparton
 458          Endif
 459          CtLhCtq6Pdf = 0D0
 460          Return
 461       Endif
 462
 463       CtLhCtq6Pdf = CtLhPartonX6 (iset,Iparton, X, Q)
 464       if(CtLhCtq6Pdf.lt.0.D0)  CtLhCtq6Pdf = 0.D0
 465
 466       Return
 467
 468 C                             ********************
 469       End
 470       SUBROUTINE CtLhPOLINT4 (XA,YA,N,X,Y,DY)
 471 c fast version of polint, valid only for N=4
 472 c Have explicitly unrolled the loops.
 473
 474       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 475
 476       PARAMETER (NMAX=4)
 477       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 478
 479         if(n .ne. 4) then
 480            print *,'fatal CtLhPolint4 call',n
 481            stop
 482         endif
 483
 484       NS=1
 485       DIF=ABS(X-XA(1))
 486
 487         DIFT=ABS(X-XA(1))
 488         IF (DIFT.LT.DIF) THEN
 489           NS=1
 490           DIF=DIFT
 491         ENDIF
 492         C(1)=YA(1)
 493         D(1)=YA(1)
 494
 495         DIFT=ABS(X-XA(2))
 496         IF (DIFT.LT.DIF) THEN
 497           NS=2
 498           DIF=DIFT
 499         ENDIF
 500         C(2)=YA(2)
 501         D(2)=YA(2)
 502
 503         DIFT=ABS(X-XA(3))
 504         IF (DIFT.LT.DIF) THEN
 505           NS=3
 506           DIF=DIFT
 507         ENDIF
 508         C(3)=YA(3)
 509         D(3)=YA(3)
 510
 511         DIFT=ABS(X-XA(4))
 512         IF (DIFT.LT.DIF) THEN
 513           NS=4
 514           DIF=DIFT
 515         ENDIF
 516         C(4)=YA(4)
 517         D(4)=YA(4)
 518
 519
 520       Y=YA(NS)
 521       NS=NS-1
 522
 523
 524           HO=XA(1)-X
 525           HP=XA(2)-X
 526           W=C(2)-D(1)
 527           DEN=W/(HO-HP)
 528           D(1)=HP*DEN
 529           C(1)=HO*DEN
 530
 531
 532           HO=XA(2)-X
 533           HP=XA(3)-X
 534           W=C(3)-D(2)
 535           DEN=W/(HO-HP)
 536           D(2)=HP*DEN
 537           C(2)=HO*DEN
 538
 539
 540           HO=XA(3)-X
 541           HP=XA(4)-X
 542           W=C(4)-D(3)
 543           DEN=W/(HO-HP)
 544           D(3)=HP*DEN
 545           C(3)=HO*DEN
 546
 547         IF (2*NS.LT.3)THEN
 548           DY=C(NS+1)
 549         ELSE
 550           DY=D(NS)
 551           NS=NS-1
 552         ENDIF
 553         Y=Y+DY
 554
 555           HO=XA(1)-X
 556           HP=XA(3)-X
 557           W=C(2)-D(1)
 558           DEN=W/(HO-HP)
 559           D(1)=HP*DEN
 560           C(1)=HO*DEN
 561
 562           HO=XA(2)-X
 563           HP=XA(4)-X
 564           W=C(3)-D(2)
 565           DEN=W/(HO-HP)
 566           D(2)=HP*DEN
 567           C(2)=HO*DEN
 568
 569         IF (2*NS.LT.2)THEN
 570           DY=C(NS+1)
 571         ELSE
 572           DY=D(NS)
 573           NS=NS-1
 574         ENDIF
 575         Y=Y+DY
 576
 577           HO=XA(1)-X
 578           HP=XA(4)-X
 579           W=C(2)-D(1)
 580           DEN=W/(HO-HP)
 581           D(1)=HP*DEN
 582           C(1)=HO*DEN
 583
 584         IF (2*NS.LT.4-3)THEN
 585           DY=C(NS+1)
 586         ELSE
 587           DY=D(NS)
 588           NS=NS-1
 589         ENDIF
 590         Y=Y+DY
 591
 592       RETURN
 593       END
 594       SUBROUTINE CTLHPOLINT3 (XA,YA,N,X,Y,DY)
 595 c fast version of polint, valid only for N=3
 596 c Have explicitly unrolled the loops.
 597       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 598       PARAMETER (NMAX=3)
 599       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 600         if(n .ne. 3) then
 601            print *,'fatal CtLhPolint3 call',n
 602            stop
 603         endif
 604       NS=1
 605       DIF=ABS(X-XA(1))
 606         DIFT=ABS(X-XA(1))
 607         IF (DIFT.LT.DIF) THEN
 608           NS=1
 609           DIF=DIFT
 610         ENDIF
 611         C(1)=YA(1)
 612         D(1)=YA(1)
 613         DIFT=ABS(X-XA(2))
 614         IF (DIFT.LT.DIF) THEN
 615           NS=2
 616           DIF=DIFT
 617         ENDIF
 618         C(2)=YA(2)
 619         D(2)=YA(2)
 620         DIFT=ABS(X-XA(3))
 621         IF (DIFT.LT.DIF) THEN
 622           NS=3
 623           DIF=DIFT
 624         ENDIF
 625         C(3)=YA(3)
 626         D(3)=YA(3)
 627       Y=YA(NS)
 628       NS=NS-1
 629           HO=XA(1)-X
 630           HP=XA(2)-X
 631           W=C(2)-D(1)
 632           DEN=W/(HO-HP)
 633           D(1)=HP*DEN
 634           C(1)=HO*DEN
 635           HO=XA(2)-X
 636           HP=XA(3)-X
 637           W=C(3)-D(2)
 638           DEN=W/(HO-HP)
 639           D(2)=HP*DEN
 640           C(2)=HO*DEN
 641         IF (2*NS.LT.2)THEN
 642           DY=C(NS+1)
 643         ELSE
 644           DY=D(NS)
 645           NS=NS-1
 646         ENDIF
 647         Y=Y+DY
 648           HO=XA(1)-X
 649           HP=XA(3)-X
 650           W=C(2)-D(1)
 651           DEN=W/(HO-HP)
 652           D(1)=HP*DEN
 653           C(1)=HO*DEN
 654         IF (2*NS.LT.1)THEN
 655           DY=C(NS+1)
 656         ELSE
 657           DY=D(NS)
 658           NS=NS-1
 659         ENDIF
 660         Y=Y+DY
 661       RETURN
 662       END