#include "bwrm.h"
#include <momoki/r22.h>
#include <momoki/r114.h>
#include <momoki/bwr.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define	OMEGA_R114	0.234	/* 偏心係数 R114 */
#define	OMEGA_R22	0.205	/* 偏心係数 R22  */
#define	M_22_114	0.976	/* 異分子間(R22-R114)相互作用パラメータ */

double mij_f(double y, double t){return(M_22_114);}
double h0R22mole(double t)      {return(h0R22(t)*M_R22);}
double h0R114mole(double t)     {return(h0R114(t)*M_R114);}

double DummyCpid(double t)      {return(0.0);}
double DummySid(double t)       {return(0.0);}


int
main( argc , argv )
int	 argc ;
char *argv[] ;
{
  double t, p, h, hv, hl, tdew,tbub, pv,pl, y22, tcal, yb, x;
  BWRM    r22, r114;
  BWRM2   m;

  InitBWRM(&r22, TcR22,  RHOcR22/M_R22, M_R22,
		   bwrmOMEGA_CONST, OMEGA_R22, 300.0, DummyCpid, h0R22mole, DummySid);
  InitBWRM(&r114, TcR114, RHOcR114/M_R114, M_R114,
		   bwrmOMEGA_CONST, OMEGA_R114, 300.0, DummyCpid, h0R114mole, DummySid);
  yb = 0.5;
  InitBWRM2(&m, &r22, &r114, yb, mij_f);
  CalBWRM2(&m, 'v',  yb);  CalBWRM2(&m, 'l',  yb);

  t =35 + 273.15;
  for ( x = 0.01; x < 1.0; x += 0.1){
	printf("%f %f %f %f\n", x, t - 273.15,
         bwrmTs(&r22, PsR22(t)*1E+6) - 273.15,
		 bwrmTs(&r114, PsR114(t)*1E+6) - 273.15);
  }
  x = 0.5;
  p = PcR22*x + PcR114*(1-x);
  t = TcR22*x + TcR114*(1-x);
  InitBWRM2(&m, &r22, &r114, x, mij_f);
  CalBWRM2(&m, 'v',  x);  CalBWRM2(&m, 'l',  x);

  p *= 8e+5;
  printf ("%f %f %f \n", PcR22, p*1e-6, PcR114);
  tcal = bwrmTyv_Pyl(&m, p, x);
  printf (" %f : %f", tcal, p*1e-6);
  printf ("\n");
  tcal = bwrmTyl_Pyv(&m, p, x);
  printf (" %f : %f", tcal, p*1e-6);
  printf ("\n");
  
  t = 300.0;
  p = bwrmPyl_Tyv(&m, t, x);
  tcal = bwrmTyv_Pyl(&m, p, m.l.y[1]);
  printf (" %f : %f %f %f (%d)", t, tcal, m.l.y[1], p*1e-6, m.b.tmp.n_iteration);
  printf ("\n");

  p = bwrmPyv_Tyl(&m, t, x);
  tcal = bwrmTyl_Pyv(&m, p, m.v.y[1]);
  printf (" %f : %f %f %f (%d)", t, tcal, m.v.y[1], p*1e-6, m.b.tmp.n_iteration);
  printf ("\n");

  for(t = 360; t < 380; t+=1.0){
	p = bwrmPyl_Tyv(&m, t, x);
	printf (" %f %f %f (%d)", t,  m.v.y[1], p*1e-6, m.b.tmp.n_iteration);
	printf ("\n");
  }
	
}