#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include "globales1.h"

using namespace std; 

int j = 1;

double f1(double u, void * params){
   double *par =  (double *) params;
   double f= 1/par[9]*(1/(2*par[8]))*pow(((1/u)-par[7])*((1/u)-par[7])-par[7]*par[7],1/2)*
  ( par[6]*pow(log(1/u),par[4])/((1-par[7]*u)*(1/par[9]-u)));
   return f;
}

double integral3log (double s, const double *parm){
double par[10];
   par[0]=parm[0];
   par[1]=parm[1];
   par[2]=parm[2];
   par[3]=parm[3];
   par[4]=parm[4];
   par[5]=parm[5];
   par[6]=parm[6];
   par[7]=2*mp*mp;
   par[8]=mp;

   gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);

   int niter = 100000;
   double liminf = 0.0;
   double limsup =1/s-1e-9;
   double errabs = 1e-11;
   double errrel = 1e-11;
   double result, error;

   gsl_function F1;
   F1.function = &f1;
   F1.params = par;
   par[9]=s;

   gsl_integration_qags (&F1, liminf, limsup, errabs, errrel , niter, w, &result, &error);

   return result;

}


double f2 (double u, void * params){
double *par =  (double *) params;
   double f = - 1/par[9]*(1/(2*par[8]))*pow(((1/u)-par[7])*((1/u)-par[7])-par[7]*par[7],1/2)*
      par[6]*pow(log(1/u),par[4])/(1-par[7]*u);
   return f;

}

double integral3log1 (double s,const double *parm){
   double par[10];
   par[0]=parm[0];
   par[1]=parm[1];
   par[2]=parm[2];
   par[3]=parm[3];
   par[4]=parm[4];
   par[5]=parm[5];
   par[6]=parm[6];
   par[7]=2*mp*mp;
   par[8]=mp;
   gsl_integration_workspace * w = gsl_integration_workspace_alloc (100000);

   int niter = 100000;
   double liminf = 1/s-1e-9;
   double limsup = 0.25/(mp*mp);
   double errabs = 1e-11;
   double errrel = 1e-11;
   double result1, error1;
//double  par [10] = {41.27, 23.12, -0.4, -0.54, 2.169, 29.29, 0.13, 2*mp*mp, mp,s};

   gsl_function F2;
   F2.function = &f2;
   F2.params = par;
   par[9]=s;
 

 gsl_integration_qawc (&F2, liminf, limsup, 1/s, errabs, errrel, niter, w, &result1, &error1);

 gsl_integration_workspace_free(w);

  
 return result1;

}

double rhopp(double s,const double *parm){

   double p=(1/(2*mp)*sqrt((s-2*ma*mp)*(s-2*ma*ma)-4*mp*mp*ma*ma));
   double suma = -integral3log(s,parm)-integral3log1(s,parm);
         suma=suma*(s-(ma*ma)-(mp*mp))/(M_PI*p);
   	 return suma;  
}


int main(){
double par[8]={40.27,23.12,0.4,0.54,2.16976,29.29,0.13};
double s[1000000];  

  for (int j = 1; j <= 4000; j++ ){
    s[j] = 4000+j*2; 


    cout<<rhopp(s[j],par)<<endl;
    
 
 }

 
  
  return 0;
}