/* This file calculates the strong coupling constant of QCD. */

/**************************************************************************/
/*                                                                        */
/* IMPORTANT : all dimensional quantities in this program                 */
/*             are in units GeV                                           */  
/*                                                                        */
/**************************************************************************/


#include <strong_coupling.h>


/* the function alpha_strong() returns the running strong coupling constant 
   as function of the negative of the momentum transfer squared (in GeV^2). */

double alpha_strong(int choice, double Q_sqr)
{
     double beta, coupling, running;
     
     beta = 11. - 2./3. * N_FLAVOURS;
     running = 4. * PI / (beta * log(Q_sqr / pow(LAMBDA1_QCD, 2.)));
     
     switch(choice)
     {
     case 1 : /* constant coupling */
            {
             coupling = ALPHA_S;
             break;
            }
     case 2 : /* running coupling */
            {
             coupling = running;
             break;
            }
     case 3 : /* running coupling, 
                 which is frozen at 0.5 at low momentum transfer */
            {
             if(running > 0.5) 
               coupling = 0.5;
             else
               coupling = running;           
             break;
            } 
     case 4 : /* analytical model for the strong coupling which is IR finite, 
                 D.V.Shirkov and I.L.Solovtsov, Phys.Rev.Lett.79 (1997)1209. */
            {
             double lambda_anal;

	     lambda_anal = .280001;

             if (Q_sqr == .0784)
               coupling = 2. * PI / beta; 
             else
             {
              if(Q_sqr < pow(10.,-6.))
	      { 
               coupling = 4. * PI / beta; 
	      }
              else
	      {
               coupling = 4. * PI / beta 
                        * ( 1./log(Q_sqr/pow(lambda_anal,2.)) 
                          + pow(lambda_anal,2.)/(pow(lambda_anal,2.) - Q_sqr)
                          );
              } 
             }
             break;
            }
     }
            
     return coupling;
}