#include <includes.h>
#include <utils.h>
#include <methnum.h>

#include "s_param.h"
#include "integral.h"

double Ivac     ( Param * P , double m )
{
  double m2 = m*m ; 
  double LE = sqrt ( m2 + P->L2 );
  return 0.5 * ONE_OVER_8PI_2 * ( P->L * LE * ( m2 + 2*P->L ) + m2*m2 * log ( m/(P->L + LE) ) ) ;
}

double dm_Ivac  ( Param * P , double m )
{
  double m2 = m*m ; 
  double LE = sqrt ( m2 + P->L2 );
  
  return m * ONE_OVER_4PI_2 * ( P->L * LE + m2 * log ( m/(P->L + LE) ) );
}

double dm2_Ivac ( Param * P , double m )
{
  double m2 = m*m ; 
  double LE = sqrt ( m2 + P->L2 ); 
  
  return ONE_OVER_4PI_2 * ( P->L*(3*m2 + P->L2)/LE + 3*m2 * log ( m / (P->L + LE) ) ) ;
}

double dm3_Ivac ( Param * P , double m )
{
  double m2 = m*m ; 
  double LE2 =  m2 + P->L2 ; 
  double LE  = sqrt( LE2 );  
  
  return 3 * m * ONE_OVER_2PI_2 * ( P->L*( 3*m2 + 4*P->L2 )/(3*LE2*LE) + log ( m / (P->L + LE) ) );
}

double Imed      ( Param * P , double m , double T , double mu )
{
  double m2 = m*m ;
  double b = 1./T ;
  double integ ( double p )
  {
    double p2 = p*p ;
    double E2 = p2 + m2 ; 
    double E = sqrt ( E2 );

    double x = -(E - mu) * b ;
    double y = -(E + mu) * b ;

    double a = log ( 1 + exp ( x ) ) ;
    double b = log ( 1 + exp ( y ) ) ;
    
    return p2 * ( a + b );
  }
  double I = ONE_OVER_2PI_2 * integ_dp ( integ , 0. , P->L , cutoff );
  return I ; 
}

double dm_Imed   ( Param * P , double m , double T , double mu )
{
  double m2 = m*m ;
  double b = 1./T ;
  double integ ( double p )
  {
    double p2 = p*p ;
    double E2 = p2 + m2 ; 
    double E = sqrt ( E2 );

    double x = (E - mu) * b ;
    double y = (E + mu) * b ;
    
    double ex = exp ( x ) ; 
    double ey = exp ( y );

    double f = 1. / ( 1 + ex );
    double fb = 1. / ( 1 + ey );
    
    return p2 * ( - f - fb ) / E ; 
  }

  double I = integ_dp ( integ , 0. , P->L , cutoff );
  return m * I ;
}
double dT_Imed   ( Param * P , double m , double T , double mu )
{

}
double dmu_Imed  ( Param * P , double m , double T , double mu );

/* double dm2_Imed  ( Param * P , double m , double T , double mu ); */
/* double dmT_Imed  ( Param * P , double m , double T , double mu ); */
/* double dmmu_Imed ( Param * P , double m , double T , double mu ); */

/* double dT2_Imed  ( Param * P , double m , double T , double mu ); */
/* double dTmu_Imed ( Param * P , double m , double T , double mu ); */

/* double dmu2_Imed ( Param * P , double m , double T , double mu ); */

