#include <stdio.h>
#include <math.h>
// Calcula isentálpicas de van der Waals
/* 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,s,h,u;
  double dv,vmax, vmin;
  double c;
  char hP;
 FILE *archivo;
 printf(" Constante c del gas:");
 scanf("%lf",&c);
 printf(" Entalpía (h/RTc):");
 scanf("%lf",&h);
 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_isoental.dat","w");
 fprintf(archivo,"Isoentálpicas 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,"h/RTcr = %lf, Energías/RTcr para c = %lf\n",h,c);
 fprintf(archivo,"  P/Pcr    v/vcr     T/Tcr     s/R       u/RTcr\n");
 fprintf(archivo,"==============================================\n");
 while(v<=vmax)
 {
  T =(h+9./4./v)/(c+3*v/(3*v-1)); //Obtenemos T a partir de h dada
  p= 8*T/(3*v-1)-3/v/v;  //Obtenemos P
  s = log(3.*v-1)+c*log(T);
  u = c*T -9./8./v;
  fprintf(archivo,"%lf %lf %lf %lf %lf\n",p,v,T,s,u);
  v+= dv;                
 }
 fclose(archivo);

}