alphas.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file alphas.F

      FUNCTION alphas (Q2, NAORD)
*********************************************************************
*                                                                   *
*   THE ALPHA_S ROUTINE.                                            *
*                                                                   *
*   INPUT :  Q2    =  scale in GeV**2  (not too low, of course);    *
*            NAORD =  1 (LO),  2 (NLO).                             *
*                                                                   *
*   OUTPUT:  alphas_s/(4 pi) for use with the GRV(98) partons.      *  
*                                                                   *
*******************************************************i*************
*
      IMPLICIT DOUBLE PRECISION (a - z)
      INTEGER nf, k, i, naord
      dimension lambdal(3:6),  lambdan(3:6), q2thr(3)
*
*...HEAVY QUARK THRESHOLDS AND LAMBDA VALUES :
      DATA q2thr   /  1.960,  20.25,  30625. /
      DATA lambdal / 0.2041, 0.1750, 0.1320, 0.0665 /
      DATA lambdan / 0.2994, 0.2460, 0.1677, 0.0678 /
*
*...DETERMINATION OF THE APPROPRIATE NUMBER OF FLAVOURS :
      nf = 3
      DO 10 k = 1, 3
      IF (q2 .GT. q2thr(k)) THEN
         nf = nf + 1
      ELSE
          go to 20
       END IF
  10   CONTINUE
*
*...LO ALPHA_S AND BETA FUNCTION FOR NLO CALCULATION :
  20   b0 = 11.- 2./3.* nf
       b1 = 102.- 38./3.* nf
       b10 = b1 / (b0*b0)
       IF (naord .EQ. 1) THEN
         lam2 = lambdal(nf) * lambdal(nf)
         alp  = 1./(b0 * dlog(q2/lam2))
         go to 1
       ELSE IF (naord .EQ. 2) then
         lam2 = lambdan(nf) * lambdan(nf)
         b1 = 102.- 38./3.* nf
         b10 = b1 / (b0*b0)
       ELSE
         WRITE (6,91)
  91     FORMAT ('INVALID CHOICE FOR ORDER IN ALPHA_S')
         stop
       END IF
*
*...START VALUE FOR NLO ITERATION :
       lq2 = dlog(q2 / lam2)
       alp = 1./(b0*lq2) * (1.- b10*dlog(lq2)/lq2)
*
*...EXACT NLO VALUE, FOUND VIA NEWTON PROCEDURE :
       DO 2 i = 1, 6
       xl  = dlog(1./(b0*alp) + b10)
       xlp = dlog(1./(b0*alp*1.01) + b10)
       xlm = dlog(1./(b0*alp*0.99) + b10)
       y  = lq2 - 1./ (b0*alp) + b10 * xl
       y1 = (- 1./ (b0*alp*1.01) + b10 * xlp
     1       + 1./ (b0*alp*0.99) - b10 * xlp) / (0.02d0*alp)
       alp = alp - y/y1
  2    CONTINUE
*
*...OUTPUT :
  1    alphas = alp
       RETURN
       END