#include <math.h>
#include <default.h>
#include <complex.h>
#include <vector.h>
#include <four_vector.h>
#include <matrix.h>
#include <gamma.h>
#include <angular.h>
#include <int.h>
#include <int_complex.h>
#include <t_gaN_piN_regge.h>

main()
{

  int choice_alpha;
double tg;
dcomplex pion_regge_trajectory(int choice_alpha, double t);
dcomplex rho_regge_trajectory(int deg_rho, double t);
dcomplex api;

 for (choice_alpha=1; choice_alpha<=3; choice_alpha++){
   printf(" choice_alpha=%d\n",choice_alpha);

   for (tg=-5.0; tg<=5.0; tg=tg+0.1) {
     api=pion_regge_trajectory(choice_alpha, tg);
     /*  printf(" %f  %f   %f\n",tg,api.re,api.im); */
   }
 }
}

dcomplex pion_regge_trajectory(int choice_alpha, double t)
{
     double alpha_lin, alpha, alpha_sat;
     dcomplex traj;
     int deg_pi;

     switch(choice_alpha)
     {
     case 1: /* linear trajectory */
           {
           alpha = ALPHAp_PI * (t - pow(PIM_M, 2.));

/*           alpha = .7 * (t - pow(PIM_M, 2.));   */
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 2: /* trajectory with saturation at SAT_PI */
           {
           alpha_lin = ALPHAp_PI * (t - pow(PIM_M, 2.));
           alpha_sat = SAT_PI;

           if ( alpha_lin > alpha_sat)
              alpha = alpha_lin;
           else
              alpha = alpha_sat;  

	   /*           alpha = alpha_sat;*/

	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 3: /* interpolating pion trajectory of Sergeenko,
                Z.Phys.C64 (1994) 315 */
           {
           double mass_q, V_o, kappa, alpha_s, np,
                  lam, c1, c2, c3, pt, qt, Qt, beta;

           mass_q = 0.231;
           V_o = - 0.5 / (4. * mass_q);
           kappa = 0.17;
           alpha_s = 4./3. * 0.892;

           np = 0.;

           lam = (4. * mass_q * (mass_q + V_o) - t) / (8. * kappa);
           c1 = 4. * np + 3.5 - alpha_s + lam;
           c2 = pow(np + 1., 2.)
              + 2. * (np + 1.) * (2. * np + 1.5 - alpha_s + lam);
           c3 = pow(np + 1., 2.) * (2. * np + 1.5 - alpha_s + lam)
              - pow(alpha_s * mass_q, 2.) / (8. * kappa);
           pt = - 1./3. * pow(c1, 2.) + c2;
           qt = 2./27. * pow(c1, 3.) - 1./3. * c1 * c2 + c3;
           Qt = 1./27. * pow(pt,3.) + 1./4. * pow(qt,2.);

           beta = acos(- qt / (2. * sqrt(- pow(pt,3.)/27.)) );

           if (Qt <= 0.)
              alpha = 2. * sqrt(- pt/3.) * cos(beta/3.) - c1/3.;
           if (Qt > 0.)
              alpha = pow(- 0.5 * qt + sqrt(Qt), 1./3.)
                    + pow(- 0.5 * qt - sqrt(Qt), 1./3.) - c1/3.;
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     }
}

dcomplex rho_regge_trajectory(int deg_rho, double t)
{
  int choice_alpha;
     double alpha_lin, alpha_sat, alpha;
     dcomplex traj;


     /* --gh--
     choice_alpha = 2;
*/
     choice_alpha = 4; 

     switch(choice_alpha)
     {
     case 1: /* linear trajectory which passes exactly through rho pole */
           {
           alpha = 1. + ALPHAp_RHO * (t - pow(RHO_M, 2.));
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 2: /* linear trajectory */
           {
           alpha = ALPHAo_RHO + ALPHAp_RHO * t;
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 3: /* experimentally determined rho traj :
                Barnes et al., PRL37 (1976) 76 */
           {
           alpha = 0.481 + 0.928 * t + 0.205 * pow(t, 2.);
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 4: /* interpolating rho trajectory of Sergeenko,
                Z.Phys.C64 (1994) 315 */
           {
           double mass_q, V_o, kappa, alpha_s, np,
                  lam, c1, c2, c3, pt, qt, Qt, beta;

           mass_q = 0.33;
           V_o = - 1.35 / (4. * mass_q);
           kappa = 0.145;
           alpha_s = 4./3. * 0.91;

           np = 0.;

           lam = (4. * mass_q * (mass_q + V_o) - t) / (8. * kappa);
           c1 = 4. * np + 3.5 - alpha_s + lam;
           c2 = pow(np + 1., 2.)
              + 2. * (np + 1.) * (2. * np + 1.5 - alpha_s + lam);
           c3 = pow(np + 1., 2.) * (2. * np + 1.5 - alpha_s + lam)
              - pow(alpha_s * mass_q, 2.) / (8. * kappa);
           pt = - 1./3. * pow(c1, 2.) + c2;
           qt = 2./27. * pow(c1, 3.) - 1./3. * c1 * c2 + c3;
           Qt = 1./27. * pow(pt,3.) + 1./4. * pow(qt,2.);

           beta = acos(- qt / (2. * sqrt(- pow(pt,3.)/27.)) );

           if (Qt <= 0.)
              alpha = 2. * sqrt(- pt/3.) * cos(beta/3.) - c1/3.;
           if (Qt > 0.)
              alpha = pow(- 0.5 * qt + sqrt(Qt), 1./3.)
                    + pow(- 0.5 * qt - sqrt(Qt), 1./3.) - c1/3.;
	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     case 5: /* trajectory with saturation at SAT_RHO */
           {
	     /*
           alpha_lin = 1. + ALPHAp_RHO * (t - pow(RHO_M, 2.));
	   /*
           alpha_lin = ALPHAo_RHO + ALPHAp_RHO * t;
           alpha_sat = SAT_RHO;

           if ( alpha_lin > alpha_sat)
              alpha = alpha_lin;
           else
              alpha = alpha_sat; 

	   /*           alpha = alpha_sat;*/

	   printf(" t=%f alpha=%f\n",t,alpha);
           break;
           }
     }
}