#include <stdio.h>
#include <math.h>
double alfa; // alfa = 4*epsilon/(k_B*T)
double Fun(double x) // x= sigma/r
{// Lennard-Jones
double x6,x2;
x2=x*x;
x6=x2*x2*x2;
return (1.-exp(-alfa*(x6*(x6-1))))/x2/x2;
}       
main()
{
      double epsilon, sigma, pi=3.14159265358969;
      double T,dT,T1,T2,B;
      double x1=10.,x2=0.1,x,dx=0.001,x23,x29;
       FILE *archivo;
printf("\nPotencial de Lennard-Jones: U(r) = 4e{(s/r)^12-(s/r)^6}");
printf("\nDefinimos x= s/r)");
printf("\n el programa calcula B(T) resolviendo numéricamente");
printf("\n B(T)= 2*pi*NA* integral[1-exp{-U(r)/kBT}]r^2dr");
printf("\nde cero a infinito");
printf("\n\n Resultados en archivo de texto: 'Lennard.dat'");
printf("\n\nparámetros epsilon/k_B (K) y sigma(A): ");
scanf("%lf %lf",&epsilon,&sigma);
//epsilon=5.; //Ne: epsilon/k_B
//sigma=2.74; //Angstrom
T1=epsilon;
T2=epsilon*5.;
dT= (T2-T1)/100.;
archivo=fopen("Lennard.dat","w");
fprintf(archivo,"#epsilon/KB: %lf (K), sigma: %lf (A)",epsilon,sigma);
fprintf(archivo,"\n# T(K)  1/T(K-1)   B(cm3/mol)");
fprintf(archivo,"\n#==============================");
for (T=T1;T<=T2;T+=dT)
{
alfa = 4*epsilon/T;
//x1=x2=1.;
B=0.;
if(fabs(x1-x2)/x1>1.e-5)
{

B=0.;
for (x=x2;x<=x1;x+=dx) B+= Fun(x);
B=dx*(B-Fun(x1)/2.-Fun(x2)/2);
}
x23=x2*x2*x2;
x29=x23*x23*x23;
B=B+1./3./x1/x1/x1+alfa*(x29/9.-x23/3.);
B=B*sigma*sigma*sigma*2*pi*0.6022; //cm^3/mol
fprintf(archivo,"\n%.2lf %lf %lf",T,1./T,B);
//getchar();
}
 fclose(archivo);

}