//======================================================
// Thermodynamic Package in java
// Class Ga //Properties of perfect Gases
// Dr. Turhan Coban
// File Name   : GasA.java
// This file contains the GasA class
// this class sets basic properties of perfect GasAes
// required data is read from GasA.dat.
// =====================================================
// Description  : This file contains the GasA class
//               class GasA calculates thermophysical properties of
//               perfect GasAses
//               following properties can be calculated
//               T()     : Temperature degree K
//               h(T)    : enthalpy  KJ/kmol
//               hf      : formation enthalpy KJ/kg
//               ht(T)   : total enthalpy KJ/kg (h+hf)
//               M       : molar mass kg/kmol
//               HT(T)   : total enthalpy KJ : M*ht(T)
//               P()     : presuure bar
//               s(T,P)  : entropy KJ/kmol K
//               Cp(T)   : specific heat at constant pressure KJ/kmol K
//               Cv(T)   : specific heat at constant volume KJ/kg K
//               gamma(T): adiabatic constant Cp/Cv
//               u(T)    : Internal energy KJ/kmol
//               c(T)    : speed of sound m/s
//               vis(T)  : viscosity      Ns/m^2
//               k(T)    : thermal conductivity KJ/kg K
// DATA FILE DEFINATION
// GasA datas are written in the data file "GasA.dat"
// if GasA data is not given in the data file, it can be curve fitted
// and added to the data file. Additional curve fitting programs supplied
// in Numerical Analysis package. Each data has the following form :
//------------------
// GasAName
// n_equation M h0 hf sf
// xa[0] xb[0] xc[0] xd[0] tl[0] th[0]
// .........
// xa[n_equation-1] xb[n_equation-1]......th[n_equation-1]
// n_vis
// xvis[0]
// .........
// xvis[n_vis-1]
// n_k
// xk[0]
// .........
// xk[n_k-1]
//-------------------
// unit of the xa :  kcal/kmol
// note : if any curvefitting applied for a new GasA temperature values
// in K nad enthalpy values in the unit of Kcal.kmol should be supply
// for the Cp curve fitting
//============================================================
//  VARIABLE IDENTIFICATION
// all the variables that type is not defined is a double variable
// PROTECTED VARIABLES :
// xa,xb,xc,xd ,tl,th : double pointers. This values used to calculate
// specific heat at constant pressure  from the following equation :
// Cp(T) = xa[i]+xb[i]*1e-3*T+xc[i]*1.0e5/T^2+xd[i]*1e-6*T^3
// where tl[i] <= T <= th[i]
// n_equation : number of equations (xa,xb,xc,xd,tl,th) for a GasA
// xvis : real pointers to define viscosity according to formula :
// vis(T)=sum(xvis(i)*T^i) , for(i=0;i<n_vis;i++)
// n_vis : number of coefficients in polynomial viscosity curve fitting
// xk :real pointers to define thermal conductivity according to formula
// k(t)=sum(xk(i)*T^i) ,  for(i=0;i<n_k;i++)
// PUBLIC VARIABLES :
// GasAName : name of the GasA example : H2O   :variable class str
// (class str is defined at file cstr.h and str.cpp, written by Timotyhy A. Budd)
// M : mol number of the GasA example mol number of H2O is 18.016 kg/kmol
// h0 : value of enthalpy at 298 K in the unit of Kcal/kmol
// hf : formation enthalpy at 298 K in the unit of Kcal/kmol
// sf : value of enthalpy at 298 K and 1 bar pressure
// N  : molar weight of the GasA, kmol
// ierror : integer variable, error flag
//============================================================

import java.io.*;
import java.awt.*;
import java.applet.Applet;
import Text;
import Atom;

class GasA extends Applet
{
double xa[],xb[],xc[],xd[],tl[],th[];
int n_equation;
int n_vis;
int n_k;
double xvis[];
double xk[];
int n_atom;
String GasAName;
Atom a[];
double M;  // molar mass of atom kg/kmol
double h0; // enthalpy at T=298 K
double hf; // enthalpy of formation
double sf; // entropy of formation    kJ/kmol K
double N;  // moles of GasA kmol
int ierror;

// definations of class functions
// constructors
//===================================================================
public GasA()
{
//empty constructor
int i;
GasAName="******************************";
n_atom=1;
a=new Atom[n_atom];
M=0;
N=1.0;
h0=0;
hf=0;
sf=0;
n_equation=6;
n_vis=10;
n_k=10;
xa=new double[6];
xb=new double[6];
xc=new double[6];
xd=new double[6];
tl=new double[6];
th=new double[6];
xvis=new double[10];
xk=new double[10];
for(i=0;i<n_equation;i++)
{xa[i]=0.0;xb[i]=0.0;xc[i]=0.0;xd[i]=0.0;tl[i]=293.0;th[i]=293.0;}
for(i=0;i<n_vis;i++)
{xvis[i]=0.0;}
for(i=0;i<n_k;i++)
{xk[i]=0.0;}
}
//===================================================================

public GasA(String gName,double Nnew) throws IOException
{
int i;
String aName;
double aN;
ierror=1;
N=Nnew;
String tempGasAName="";
//ifstream fin("GasA1.dat",ios::in);
DataInputStream fin=new DataInputStream(new FileInputStream("GasA.dat"));
try{
while(fin!=null)
{
tempGasAName=Text.readString(fin);
if (tempGasAName.equals(gName)) {ierror=0;break;}
}
  } catch(EOFException e_eof)
     {
     showStatus("error required GasA "+tempGasAName+" is not found");
     fin.close();return;
     }
GasAName=gName;
n_atom=Text.readInt(fin);
a=new Atom[n_atom];
M=0;
for(i=0;i<n_atom;i++)
 {
 aName=Text.readString(fin);
 aN=Text.readDouble(fin);
 a[i]=new Atom(aName,aN);
 M+=a[i].mass;
 }
n_equation=Text.readInt(fin);
h0=Text.readDouble(fin);
hf=Text.readDouble(fin);
sf=Text.readDouble(fin);
hf*=4.1868;
xa=new double[n_equation];
xb=new double[n_equation];
xc=new double[n_equation];
xd=new double[n_equation];
tl=new double[n_equation];
th=new double[n_equation];

for(i=0;i<n_equation;i++)
  {
    xa[i]=Text.readDouble(fin);
    xb[i]=Text.readDouble(fin);
    xc[i]=Text.readDouble(fin);
    xd[i]=Text.readDouble(fin);
    tl[i]=Text.readDouble(fin);
    th[i]=Text.readDouble(fin);
  }
n_vis=Text.readInt(fin);
xvis=new double[n_vis];
for(i=0;i<n_vis;i++)
  {
  xvis[i]=Text.readDouble(fin);
  }
n_k=Text.readInt(fin);
xk=new double[n_k];
for(i=0;i<n_k;i++)
  {
  xk[i]=Text.readDouble(fin);
  }
}

public GasA(String gName) throws IOException
{
double Nnew=1;
int i;
String aName;
double aN;
ierror=1;
N=Nnew;
String tempGasAName="";
//ifstream fin("GasA1.dat",ios::in);
DataInputStream fin=new DataInputStream(new FileInputStream("GasA.dat"));
try{
while(fin!=null)
{
tempGasAName=Text.readString(fin);
if (tempGasAName.equals(gName)) {ierror=0;break;}
}
  } catch(EOFException e_eof)
     {
     showStatus("error required GasA "+tempGasAName+" is not found");
     fin.close();return;
     }
GasAName=gName;
//fin>>n_atom;
n_atom=Text.readInt(fin);
a=new Atom[n_atom];
M=0;
for(i=0;i<n_atom;i++)
 {
 //fin>>aName>>aN;
 aName=Text.readString(fin);
 aN=Text.readDouble(fin);
 a[i]=new Atom(aName,aN);
 M+=a[i].mass;
 }
//fin>>n_equation>>h0>>hf>>sf;
n_equation=Text.readInt(fin);
h0=Text.readDouble(fin);
hf=Text.readDouble(fin);
sf=Text.readDouble(fin);
hf*=4.1868;
xa=new double[n_equation];
xb=new double[n_equation];
xc=new double[n_equation];
xd=new double[n_equation];
tl=new double[n_equation];
th=new double[n_equation];

for(i=0;i<n_equation;i++)
  {
    xa[i]=Text.readDouble(fin);
    xb[i]=Text.readDouble(fin);
    xc[i]=Text.readDouble(fin);
    xd[i]=Text.readDouble(fin);
    tl[i]=Text.readDouble(fin);
    th[i]=Text.readDouble(fin);
  }
n_vis=Text.readInt(fin);
xvis=new double[n_vis];
for(i=0;i<n_vis;i++)
  {
  xvis[i]=Text.readDouble(fin);
  }
n_k=Text.readInt(fin);
xk=new double[n_k];
for(i=0;i<n_k;i++)
  {
  xk[i]=Text.readDouble(fin);
  }
}

public GasA(GasA g) throws IOException
{
int i;
GasAName=g.GasAName;
n_atom=g.n_atom;
a=new Atom[n_atom];
for(i=0;i<n_atom;i++)
{
a[i]=g.a[i];
}
M=g.M;
N=g.N;
h0=g.h0;
hf=g.hf;
sf=g.sf;
n_equation=g.n_equation;
n_vis=g.n_vis;
n_k=g.n_k;
xa=new double[n_equation];
xb=new double[n_equation];
xc=new double[n_equation];
xd=new double[n_equation];
tl=new double[n_equation];
th=new double[n_equation];
xvis=new double[n_vis];
xk=new double[n_k];
for(i=0;i<n_equation;i++)
{xa[i]=g.xa[i];xb[i]=g.xb[i];xc[i]=g.xc[i];xd[i]=g.xd[i];tl[i]=g.tl[i];th[i]=g.th[i];}
for(i=0;i<n_vis;i++)
{xvis[i]=g.xvis[i];}
for(i=0;i<n_k;i++)
{xk[i]=g.xk[i];}

}


//===================================================================

//change/assign new molar mass
// see also * operator for the similar function
public void changeN(double Nnew)
{
N=Nnew;
}
//===================================================================
public String toString()
{
//return the chemical symbol of the GasA
String s="";
for(int i=0;i<n_atom;i++)
  s=s+a[i].toString();
return s;
}
//===================================================================
public double vis(double T)
{
// dynamic viscosity of the GasA
double visg=0;
if(n_vis!=0.0)
{
  visg=xvis[n_vis-1];
  for(int i=n_vis-2;i>=0;i--)
   { visg=visg*T+xvis[i]; }
  visg*=1.0e-7;
}
else
  visg=0;
return visg;
}
//===================================================================

public double k(double T)
{
// thermal conductivity of the GasA
double kg;
if(n_k!=0.0)
{
  int nk=n_k-1;
  kg=xk[nk];
  for(int i=n_k-2;i>=0;i--)
   { kg+=kg*T+xk[i]; }
	kg*=1.0e-3;
}
else
  kg=0;
return kg;
}
//===================================================================

   public double h(double T)
   {
   //  enthalpy KJ/kmol
   //integration of function dh=Cp(T)*dT
   double hh = - h0;
     for(int i=0;i<n_equation;i++)
     {
	     if(((T>th[i]) && (i== (n_equation-1) ) )
		|| ((T<tl[i]) && (i== 0 ) ) )
             {
	     hh+=   xa[i]*(T- tl[i])
	          + xb[i]*1.0e-3/2.0*(T*T-tl[i]*tl[i])
		  - xc[i]*1e5*(1/T-1/tl[i])
		  + xd[i]*1e-6*(T*T*T-tl[i]*tl[i]*tl[i])/3.0;
      }
             else if((T<= th[i]) && (T> tl[i]))
				 {
	     hh+=   xa[i]*(T- tl[i])
	          + xb[i]*1.0e-3/2.0*(T*T-tl[i]*tl[i])
		  - xc[i]*1e5*(1/T-1/tl[i])
		  + xd[i]*1e-6*(T*T*T-tl[i]*tl[i]*tl[i])/3.0;
             }
             else if(T>th[i])
				 {
	     hh+=   xa[i]*(th[i]- tl[i])
	          + xb[i]*1.0e-3/2.0*(th[i]*th[i] - tl[i]*tl[i])
		  - xc[i]*1e5*(1/th[i]-1/tl[i])
		  + xd[i]*1e-6*(th[i]*th[i]*th[i]-tl[i]*tl[i]*tl[i])/3.0;
             }
	  }
   return (hh*4.1868);
   }
//===================================================================

   public double ht( double t)
   {
   return h(t)+hf;
   }
//===================================================================

   public double H(double t)
   {
   return h(t)*N;
   }
//===================================================================

   public double HT(double t)
   {
   return ht(t)*N;
   }
//===================================================================

   public double u(double T)
   {
   // internal energy KJ/kmol
   // Integration of function du = Cv(T)*dT
   return (h(T) - 8.3145*T);
   }
//===================================================================

   public double v(double T,double P)
   {
   // specific volume of the GasA m^3/kmol
   return 8314.5*T/(P*1e5);
   }

   public double v(double T)
   {
   double P=1.0;
   // specific volume of the GasA m^3/kmol
   return 8314.5*T/(P*1e5);
   }

//===================================================================
   public double c(double t)
   {
   // speed of sound m/s
   return Math.sqrt(8314.5/M*t*gamma(t));
   }
//===================================================================
   public double s(double T, double P)
   {
   //entropy KJ/kmol K
   //integration of function
   // ds = Cp(T) * dt/T - R dP/P
   double ss=sf;
   for(int i=0;i<n_equation;i++)
		{
           if( ( T > th[i]  && i==n_equation - 1 )
	     ||( T < tl[i]  && i==0              ))
		       {
	       ss+=xa[i]*Math.log(T/tl[i])
	       + xb[i]*1.0E-3*(T-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i]))
	       + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]);
               }
	   else if((T <= th[i]) && (T > tl[i]))
	       {
	       ss+=xa[i]*Math.log(T/tl[i])
	       + xb[i]*1.0E-3*(T-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i]))
               + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]);
					}
	   else if( T > th[i] )
	       {
	       ss+=xa[i]*Math.log(th[i]/tl[i])
	       + xb[i]*1.0E-3*(th[i]-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(th[i]*th[i]) - 1.0/(tl[i]*tl[i]))
	       + xd[i]*1e-6/2.0*(th[i]*th[i] - tl[i]*tl[i]);
					}
	   }
   ss*=4.1868;
   return (ss-8.3145*Math.log(P));
   }

   public double s(double T)
   {
   //entropy KJ/kmol K
   //integration of function
   // ds = Cp(T) * dt/T - R dP/P
   double ss=sf;
   for(int i=0;i<n_equation;i++)
		{
           if( ( T > th[i]  && i==n_equation - 1 )
	     ||( T < tl[i]  && i==0              ))
		       {
	       ss+=xa[i]*Math.log(T/tl[i])
	       + xb[i]*1.0E-3*(T-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i]))
	       + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]);
               }
	   else if((T <= th[i]) && (T > tl[i]))
	       {
	       ss+=xa[i]*Math.log(T/tl[i])
	       + xb[i]*1.0E-3*(T-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(T*T) - 1.0/(tl[i]*tl[i]))
               + xd[i]*1e-6/2.0*(T*T-tl[i]*tl[i]);
					}
	   else if( T > th[i] )
	       {
	       ss+=xa[i]*Math.log(th[i]/tl[i])
	       + xb[i]*1.0E-3*(th[i]-tl[i])
	       - xc[i]*1e5/2.0*(1.0/(th[i]*th[i]) - 1.0/(tl[i]*tl[i]))
	       + xd[i]*1e-6/2.0*(th[i]*th[i] - tl[i]*tl[i]);
					}
	   }
   ss*=4.1868;
   return ss;
   }
//===================================================================
public double s0(double T)
{
  return s(T);
}
//===================================================================
public double g(double T,double P)
{
  return h(T)-T*s(T,P);
}

public double g(double T)
{
  double P=1.0;
  return h(T)-T*s(T,P);
}

//===================================================================
public double g0(double T)
{
  return h(T)-T*s0(T);
}
//===================================================================
public double Cp(double T)
{
   //specific heat at constant pressure KJ/kmol K
   double cp=0.0;
   for (int i=0;i<n_equation;i++)
     {
        if( ( T > th[i]  && i==n_equation - 1 )
		  ||( T < tl[i]  && i==0              ))
		{
        cp=xa[i]+xb[i]*1.0e-3*T+xc[i]*1.0e5/T/T+xd[i]*1.0e-6*T*T;
	break;
	}
	else if((T <= th[i]) && (T > tl[i]) )
	{
	cp=xa[i]+xb[i]*1.0e-3*T+xc[i]*1.0e5/T/T+xd[i]*1.0e-6*T*T;
		  break;
	}
}
return (cp*4.1868);
}
//===================================================================
public double Cv(double T)
{
   //specific heat at constant volume KJ/kmol K
double cv;
cv=Cp(T) - 8.3145;
return cv;
}
//===================================================================
public double gamma(double T)
{
//adiabatic constant
return Cp(T)/Cv(T);
}
//===================================================================
public double T( char name,double y0,double p)
{
// name can have values h : for enthalpy
//                      u : for internal energy
//                      s : for entropy
//                      v : specific volume
// yo : the value of the variable given by variable name
double t=300;
if(name=='v') {t= p*1e5*y0/8.3145e3;}
else
  {
  double dt=0;
  int nmax=400;
  double tolerance=1.0e-8;
  for(int i=0;i<nmax;i++)
    {
// apply newtons method for finding roots of equation
    if     (name=='h') dt=-( h(t)   - y0 ) /Cp(t);
    else if(name=='u') dt=-( u(t)   - y0 ) /Cv(t);
    else if(name=='s') dt=-( s(t,p) - y0 ) /(Cp(t)/t);
    else { showStatus("wrong name defined please try h,u,s ot v");}
	 t+=dt;
//   if error range is less than tolerance, exit
    if(Math.abs(dt)<tolerance) break;
    }
  }
return t;
}

public double T( char name,double y0)
{
// name can have values h : for enthalpy
//                      u : for internal energy
//                      s : for entropy
//                      v : specific volume
// yo : the value of the variable given by variable name
double t=300;
double p=1.0;
if(name=='v') {t= p*1e5*y0/8.3145e3;}
else
  {
  double dt=0;
  int nmax=400;
  double tolerance=1.0e-8;
  for(int i=0;i<nmax;i++)
    {
// apply newtons method for finding roots of equation
    if     (name=='h') dt=-( h(t)   - y0 ) /Cp(t);
    else if(name=='u') dt=-( u(t)   - y0 ) /Cv(t);
    else if(name=='s') dt=-( s(t,p) - y0 ) /(Cp(t)/t);
    else { showStatus("wrong name defined please try h,u,s ot v");}
	 t+=dt;
//   if error range is less than tolerance, exit
    if(Math.abs(dt)<tolerance) break;
    }
  }
return t;
}
//===================================================================
public double P( char name,double y0,double t1)
{
// name can have values v : for specific volume
//                      s : for entropy
// note : for a perfect GasA enthalpy and internal energy
// is not function of pressure
// yo : the value of the variable given by variable name
if(name=='v')       return 8.3145e3*t1/y0*1e-5;
else if (name=='s') return Math.exp((s(t1,1.0)-y0)/8.3145);
else { showStatus("wrong name defined please try s or v"); return 1.0;}
}
//===================================================================
public double Pr(double t)
{
// Prandtl number
return Cp(t)*vis(t)/k(t)/M*1e3;
}
//===================================================================
public void assign(GasA g1) throws IOException
{
// assign operator (assigning a new GasA to the GasA variable)
int i;
ierror=1;
GasAName=g1.GasAName;
N=g1.N;
n_equation=g1.n_equation;
n_k=g1.n_k;
n_vis=g1.n_vis;
M=g1.M;
h0=g1.h0;
hf=g1.hf;
sf=g1.sf;
n_atom=g1.n_atom;
a=new Atom[n_atom];
M=0;
for(i=0;i<n_atom;i++)
 {
 a[i]=new Atom(g1.a[i].name,g1.a[i].N);
 }
xa=new double[n_equation];
xb=new double[n_equation];
xc=new double[n_equation];
xd=new double[n_equation];
tl=new double[n_equation];
th=new double[n_equation];
xvis=new double[n_vis];
xk=new double[n_k];
for(i=0;i<n_equation;i++)
  { xa[i]=g1.xa[i];
    xb[i]=g1.xb[i];
    xc[i]=g1.xc[i];
    xd[i]=g1.xd[i];
    tl[i]=g1.tl[i];
    th[i]=g1.th[i];
  }
for(i=0;i<n_vis;i++)
  { xvis[i]=g1.xvis[i];}
for(i=0;i<n_k;i++)
  {xk[i]=g1.xk[i];}
}

//===================================================================
public GasA multiply(double Nnew, GasA g1) throws IOException
{
GasA g2=new GasA(g1);
g2.N*=Nnew;
return g2;
}

}

//=================================================================

