#include <stdio.h>
#include <math.h>
// Calcula isobaras de van der Waals monoatómico
/* Magnitudes adimensionales:
 v =  vol(m3/mol)/3b
 P = P(Pa)/Pcr, Pcr = a/27 b^2
 T = T(K) Tcr, Tcr = (8a /27bR)
 ecuación adimensional de vdW: P = 8T(3v-1)-3/v^2
*/

main()
{double p,v,T,Cp,s,u;
  double dv,vmax, vmin, T1,s1,T2,s2;
  double c;
  char i=0;
 FILE *archivo;
 printf(" Cosntante c del gas:");
 scanf("%lf",&c);
 printf(" Presion reducida:");
 scanf("%lf",&p);
 printf("\nVolumen reducido inicial y final");
 printf("\n(calcula 1000 puntos):  ");
 scanf("%lf%lf",&vmin,&vmax);
 dv=(vmax-vmin)/1000.;
 v=vmin;
 while(v <= 1./3.) {v+= dv;}
 archivo=fopen("vdw_isobar.dat","w");
 fprintf(archivo,"Isobaras de van der Waals\n");
 fprintf(archivo,"=========================\n");
 fprintf(archivo,"Unidades: v = vol(m^3)/3b, P =P(Pa)/Pcr, Pcr = a/27b^2 \n");
 fprintf(archivo,"          T = T(K)/Tcr, Tcr = 8a/27bR \n\n");     
 fprintf(archivo,"p = %lf, Energías/RTcr para c = %lf\n",p,c);
 fprintf(archivo,"  T/Tcr    v/vcr     s/R     Cp/R       u/RTcr\n");
 fprintf(archivo,"==============================================\n");
 while(v<=vmax)
 {
  T1=T;  //T anterior
  s1=s;  // s anterior
  T =(p+3/v/v)*(3*v-1.)/8.;
  s = log(3.*v-1)+c*log(T);
  u = c*T -9./8./v;
  T2=(p+3/(v+dv)/(v+dv))*(3*(v+dv)-1.)/8.; //T siguiente
  s2= log(3.*(v+dv)-1)+c*log(T2);          //s siguiente
  if(i!=0)
  {Cp=T*(s2-s1)/(T2-T1); // Deriva numericamente
  fprintf(archivo,"%lf %lf %lf %lf %lf\n",T,v,s,Cp,u);
  }
  else 
  {Cp=0.;i=1;} // i=0 la primera vez que se ejecuta el bucle sobre v
  v+= dv;                
 }
 fclose(archivo);

}