#include <stdio.h>
#include <math.h>

/* Este programa calcula la presi¢n de vapor de un fluido de van der Waals
   los vol£mens espec¡ficos del l¡quido y gas y los l¡mites de metaestabilidad
   de ambos a varias temperaturas dadas

   Y tambien el calor latente

   Guarda los datos en archivos de texto llamados "vdw_coex.dat" y "clatente.dat"
   
   Las tempearturas, presiones y vol£menes espec¡ficos van en unidades
   de los valores cr¡ticos Tcr= 8*a/(27*b*R),  Pcr=a/(27*b^2) y vcr=3b

   Por tanto sirven para cualquier gas
   basta multiplicar los datos de salida por los valores cr¡ticos
   medidos experimentalmente (o deducidos a partir de a y b)
   para pasar a K, Pa y m^3/mol y obtener valores te¢ricos
   para un gas real dados.
*/


//Variables globales:

double t=0.999; // T/Tcr
double vgmin,vlmax,pgmax,plmin,vg,vl;
double err =1.e-8; // error maximo en las soluciones numericas

// funciones:
  /* fun2 obtiene los l¡mites de estabilidad del l¡quido y gas,
     cuando (dp/dv)t =0, m ximo y m¡nimo de la funci¢n de van der Waals p(v)
  */
  double fun2(double v)  { return 4*t*v*v*v-9*v*v+6*v-1;}


  /* vdw sirve para obtener el volumen a p y t dadas,
  de la ecuaci¢n de estado de van der waals:
  la funci¢n es :  p - 8t/(3*v-1)+3/v^2, que debe anularse */
  double vdw( double p, double v)  { return p-8*t/(3*v-1)+3/v/v;}

  /*   fun(p) es la funci¢n de la p gina 241 de Callen (corregida),
   que a temperatura dada
  DEBE ANULARSEe cuando el l¡quido y el gas tienen el mismo potencial qu¡mico
     y por tanto est n en equilibrio. Se obtienen lso vol£menes espec¡ficos del
     l¡quido y del gas para presi¢n dada y se resuelve la ecuaci¢n por el m‚todo
     de bisecci¢n. Es m s r pido el de Newton, pero este es m s sencillo */
  double fun(double p)
{  double v1,v2;
 // obtenemos el volumen de l¡quido a esa p
    v1=0.333334; v2=vlmax;
//    printf("\n vmax liquido (en fun): %lf",vlmax);
    do
     {
      vl=(v1+v2)/2.;
      if(vdw(p,v1)*vdw(p,vl)>0) v1=vl;
      else v2=vl;
     }
     while(fabs(v2-v1)>err);
 // obtenemos el volumen de gas a esa p     v1=vgmin; v2=10;
//     printf("\n vmin gas  (en fun): %lf",vgmin);
     v1=vgmin;v2=300;
     do
     {
      vg=(v1+v2)/2.;
      if(vdw(p,v1)*vdw(p,vg)>0) v1=vg;
      else v2=vg;
     }
     while(fabs(v2-v1)>err);
 //retorna el valor de la funci¢n de Callen p.241
 // return log((3*vg-1)/(3*vl-1))+9./4./t*(1./vg-1./vl)-1./(3*vg-1)+1./(3*vl-1);
   return log((3.*vg-1.)/(3.*vl-1.))+9./4./t*(1./vg-1./vl)-3.*vg/(3.*vg-1.)+3.*vl/(3.*vl-1.);
}

/*
Programa principal. Para cada temperatura,  t < 1 (t=1: ta. cr¡tica)
se hace:
1) Se determinan los l¡mites de estabilidad del l¡quido y del gas, es decir
   p liquido minima, v liquido m ximo, p gas m ximo y vgas m¡nimo.
2) los valores de p liquido minimo y p gas m ximo se toman como
   cota inferior y superior  para b£squeda de la soluci¢n de la ecuaci¢n
   de Callen , p.241 (corregido el error) por el m‚todo de bisecci¢n.
3) Se obtiene el valor de p que hace anularse fun por el m‚todo de bisecci¢n.
   Es decir, se toma el valor medio de pmax y pmin. Si la funci¢n en pmin
   tiene el mismo signo que en el valor medio, se sustituye pmin por el
   valor medio. Si no se sustituye pmax.

   En un paso intermedio la funci¢n fun obtiene los vol£menes de l¡quido  y gas
   para la presi¢n p probada.
4) Listado de resultados.
5) Se toma otra temperatura: bajando en incrementos de 0.01, hasta t=0.5)
*/
main()
{
FILE *archivo,*arch1;
int i;
double p,v,vmin,dt=0.001;   // p = P/Pcr, v = V/Vcr, t = T/Tcr
double v1,v2,p1,p2;
double c=3./2.,ul,sl,hl,gl,ug,sg,hg,gg; //c = cv/R del gas,por defecto monoatomico.
printf("\n Este programa calcula vl,vg sl,sg, hl,hg");
printf("y calor latente en un fluido ideal de van der Waals");
printf("\n Teclear la cte c: ");scanf("%lf",&c);
archivo=fopen("vdw_coex.dat","w");
arch1=fopen("clatente.dat","w");
//Escribe  encabezamiento (puede ser preferible no ponerlo)
fprintf(archivo,"Programa de transicion de fase en fluido de van der Waals");
fprintf(archivo,"\n");
fprintf(archivo,"\n\n Capac. calorifica cte Cv = cR,  c = %lf\n",c);
fprintf(archivo,"\n\nSignificado: Columna 1: T/Tcr, Columna 2: P/Pcr (eq. liqu-gas)");
fprintf(archivo,"\n Col 3: v liquido/vcr, Col 4: vgas/vcr (en equilibrio)");
fprintf(archivo,"\n Cols 5 y 6 v liquido maximo  y presion (limite estabilidad liquido)");
fprintf(archivo,"\n Cols 7 y 8 v gas  y presion (limite estabilidad gas)");
fprintf(archivo,"\n Col 9: Delta h *27b/8a = clatente, (multiplicar por 8a/27b para obtenerla en J/mol");

fprintf(archivo,"\n\n Temp\t p_vapor\t vliqu\t vgas\t vlmax\t  plmin\t vgmin\t pgmax\t hg-hl ");

fprintf(arch1,"Programa de transicion de fase en fluido de van der Waals");
fprintf(arch1,"\n\nSignificado: Columna 1: T/Tcr, Columna 2: P/Pcr");
fprintf(arch1,"\n Col 3: v liquido/vcr, Col 4: vgas/vcr");
fprintf(arch1,"\n Col 5 entropia liquido/R, Col 6:  entropia gas/R, Col 7:  entalp liq*27b/8a, Col8 : entalp gas*27b/8a");
fprintf(arch1,"\n Col 9: Delta h *27b/8a = clatente, Col 10: T * Delta s =Delta_g-Delta_h");
fprintf(arch1,"\n\nPara cada t hay tres lineas: 1a. Para presi¢n de equilibrio liquido-gas (Dg =0)");
fprintf(arch1,"\n 2a. Para presi¢n del limite de estabilidad del liquido");
fprintf(arch1,"\n 3a. Para presi¢n del limite de estabilidad del gas");

fprintf(arch1,"\n\n Capac. calorifica cte Cv = cR,  c = %lf\n",c);
fprintf(arch1,"\n T \t\t p  \t vliq \t vgas \t sl \t sg \t\t hl\t hg\t gl \t gg \t ul\t ug");
//fprintf(arch1,"\n T \t\t p  \t vliq \t vgas \t sl \t sg \t\t hl\t hg\t Dh \t TDs");
fprintf(arch1,"\n");
fprintf(archivo,"\n");
for(i=1;i<=90;i++) {fprintf(archivo,"=");fprintf(arch1,"=");}
fprintf(archivo,"\n");
fprintf(arch1,"\n");

while (t>0.5)
  {
   // determina los l¡mites de estabilidad en presi¢n y volumen
   // resolviendo la ec 4*t*v^3-9*v^2+6*v-1=0
    vmin=(3.+sqrt(9.-8.*t))/4.;
   // printf("\n vmin: %lf",vmin);
    v1=vmin; v2=5.;
    do
    {
     v= (v1+v2)/2.;
     if(fun2(v)*fun2(v1)>0) v1=v;
     else v2=v;
     }
   while((v2-v1)>err);
    vgmin=v;
   pgmax=8*t/(3*v-1)-3/v/v; // presi¢n maxima del gas a esa t
   v1=0.3334;v2=vmin;
if(t==0.95) printf("\n proceso para t = 0.95\n");
   do   
   {  v=(v1+v2)/2.;
      if(fun2(v)*fun2(v1)>0) v1=v;
      else v2=v;
     if(t==0.95) printf("\n v: %lf",v);
   }
   while((v2-v1)>err);
   vlmax=v;
   plmin=8*t/(3*v-1)-3/v/v; // presi¢n m¡nima del liquido
   p1=plmin;
   p2=pgmax;
    do  //buscamos la presi¢n ("de vapor") de equilibrio liquido-gas a esa t
   {
      p=(p1+p2)/2.;
      if(fun(p)*fun(p1)>0) p1=p;
      else p2=p;

   }
   while((p2-p1)>err);
   // Calcula energia, entropia, entalpia y e.l de gibbs del liquido y gas en equlibrio
   sl = log(3.*vl-1.)+c*log(t); // Se ha tomado s_cr=S(t=1,v=1)= ln2
   sg = log(3.*vg-1.)+c*log(t);
   hl = t*(c+3.*vl/(3.*vl-1.)) -9./4./vl;  //h = entalpia * 27b/8a
   hg = t*(c+3.*vg/(3.*vg-1.)) -9./4./vg;
   gl = hl-t*sl;
   gg = hg-t*sg;
   ul = c*t-9./8./vl;
   ug = c*t-9./8./vg;
  
   //Escribe resultados
  
   fprintf(archivo,"\n %lf %lf %lf %lf  %lf %lf %lf %lf %lf",t,p,vl,vg,vlmax,plmin,vgmin,pgmax,hg-hl);
   fprintf(arch1,"\n %lf eq %lf %lf %lf %lf %lf %lf %lf %lf %f %lf %lf",t,p,vl,vg,sl,sg,hl,hg,gl,gg,ul,ug);
      // calcula en el limite de estabilidad del líquido
   v1=vgmin;
   v2=300;
   vl=vlmax;
   p=plmin;
   while(v2-v1>err)
   { vg=(v1+v2)/2.;
    if (vdw(p,vg)*vdw(p,v1)>0) v1=vg;
    else v2=vg;
   }
   sl = log(3.*vl-1.)+c*log(t);
   sg = log(3.*vg-1.)+c*log(t);
   hl = t*(c+3.*vl/(3.*vl-1.)) -9./4./vl;  //h = entalpia * 27b/8a
   hg = t*(c+3.*vg/(3.*vg-1.)) -9./4./vg;
   gl = hl-t*sl;
   gg = hg-t*sg;
   ul = c*t-9./8./vl;
   ug = c*t-9./8./vg;
  fprintf(arch1,"\n %lf limliq %lf %lf %lf %lf %lf %lf %lf %lf %f %lf %lf",t,p,vl,vg,sl,sg,hl,hg,gl,gg,ul,ug);
 // fprintf(arch1,"\n %lf limliq %lf %lf %lf %lf %lf %lf %lf %lf %f ",t,p,vl,vg,sl,sg,hl,hg,hg-hl,t*(sg-sl));
      // calcula en el limite de estabilidad del gas
   v1=0.33335;
   v2=vlmax;
   p=pgmax;
   vg=vgmin;
   while(v2-v1>err)
   { vl=(v1+v2)/2.;
    if (vdw(p,vl)*vdw(p,v1)>0) v1=vl;
    else v2=vl;
   }
   sl = log(3.*vl-1.)+c*log(t);
   sg = log(3.*vg-1.)+c*log(t);
   hl = t*(c+3.*vl/(3.*vl-1.)) -9./4./vl;  //h = entalpia * 27b/8a
   hg = t*(c+3.*vg/(3.*vg-1.)) -9./4./vg;
   gl = hl-t*sl;
   gg = hg-t*sg;
   ul = c*t-9./8./vl;
   ug = c*t-9./8./vg;
 fprintf(arch1,"\n %lf limgas %lf %lf %lf %lf %lf %lf %lf %lf %f %lf %lf",t,p,vl,vg,sl,sg,hl,hg,gl,gg,ul,ug);
//   fprintf(arch1,"\n %lf limgas %lf %lf %lf %lf %lf %lf %lf %lf %f ",t,p,vl,vg,sl,sg,hl,hg,hg-hl,t*(sg-sl));
   fprintf(arch1,"\n");
   t=t-dt;
  }
  fclose(archivo);
  fclose(arch1);
}