//===================================================================
// File Name    : pipe.java
// Author       : Turhan Coban
// Description  :This file contains the pipe class
//               class pipe calculates pressure drops
//               and internal heat transfer of pipes
// DATA FILE DEFINATION
// pipe datas are written in the data file "pipe.dat"
// if gas data is not given in the data file, it can be curve fitted
// and added to the data file. Each data has the following form :

//------------------
// pipeName
// shape
// gasName
// d1 d2 (or only d1 depends on shape)
// mass_flow_rate
//-------------------
// a constructor is also available to define this data as passing values
//============================================================
//  VARIABLE IDENTIFICATION
// all the variables that type is not defined is a double variable
// pipeName : str class variable for name
// shape    : str class variable for shape of the cross-section
//            options :
//                    circular             
//                    rectangular         
//                    concentric_annulus 
//                    circular_segment
//                    circular_sector
//                    right_triangle
// gasName : str class variable for name of the gas mixtures
// g      : gmix class variable
// mass_flow_rate : mass flow rate of gas mixture flowing
//                  through pipe, kg/s
// eps     : pipe equivalent roughness m
//          -------------------------------------------
//          value of equivalent roughness for some pipes:
//          Pipe             equivalent rougness m
//          -------------    ---------------------
//          Riveted steel    0.9e-3 - 9e-3
//          Concrete         0.3e-3 - 3e-3
//          Wood stave       0.18e-3 - 0.9 e-3
//          Cast iron        0.26e-3
//          Galvanised iron  0.15e-3
//          Commercial steel 0.045e-3
//          wrought iron     0.045e-3
//          drawn tubing     0.0015e-3
//          Plastic,glass    0.0 (smooth)
//          -------------------------------------------
// P       : Pressure bar
// l       : length m
// dh      : hydraulic diameter m
// A       : pipe area m^2
// KL      : local pressure drop coefficient
//
//============================================================

import java.io.*;
import Gmix;
 

class pipe{
protected double d1,d2;  // m pipe dimensions
	                 // actual meaning of d1 and d2
	                 // can be change according to shape
public String pipeName;
public String shape;    //pipe shape "circular","rectangular"...
str gasName;  //name of a gas mixture
gmix *g;
double mass_flow_rate; // kg/s
double eps ;            // m
double P  ;             // Pressure bar
double l ;              // length
double lentry(double);         // hydraulic pipe entry length
double dh;              // hydraulic diameter
double A;               // pipe area m^2
double KL;               //additional local pressure drop
// pipe functions
// constructors
pipe (const str &,double length);
pipe();
pipe(const str & pn,
     const str & sh,
     const str & gn,
     double dim1,double dim2,
     double m,
     double e,
     double pressure,
     double kl,
     double length);
// pn : pipe name (search pipename through gmix.cpp file and assign gas mixture
// sh : pipe shape  "circular","rectangular"...
// gn : gas name  (search gasname through gmix.cpp file and assign gas mixture
// dim1 dim2 pipe dimensions
// e  : surface
//length : length
pipe(const str & pn,
     const str & sh,
	  gmix & gmix1,
     double dim1,double dim2,
     double m,
     double e,
     double pressure,
     double kl,
     double length);
// gmix1 : class gmix variable this variable should be previously defined
void operator = (const pipe & );
// destructor
~pipe();
//change properties
void change_length(double l1);
void change_mass_flow_rate(double);
void change_gas(const str& );
void change_gas(gmix &);
void change_P(double );
void change_KL(double);
void change_e(double);
void change_pipe (const str &,double);
void change_pipe(const str & pn,
     const str & sh,
     const str & gn,
     double dim1,double dim2,
     double m,double e,
     double pp,double k);
void change_d1(double dim1);
void change_d2(double dim2);
//properties
double Re(double);        // Reynold's number
double f(double);         // friction coefficient
double dP(double,double);        // pressure drop Bar
double dPl_e(double,double);      //exact laminar pressure drop
double velocity(double);  // m/s
double hc(double);        // convective heat transfer coefficient W/m^2K
double Nu(double);        //Nusselt's number

//access functions
double get_d1(); //getting the first dimension
double get_d2(); //getting the second dimension
		  };

 //=============================================
double f(double Re,double eod)
{
// friction factor for turbulent flow 5000<Re<10^8
double f1=log10(eod/3.7+5.74/pow(Re,0.9));
f1=1.325/(f1*f1);
return f1;
}

double fx(double X,double Re,double eod)
{
//colebrook equation to solve
// friction factor for turbulent flow 2000 < Re
double xx=2.0*log10(eod/3.7+2.51/Re*X)+X;
return xx;
}

double dfx(double X,double Re,double eod)
{
//derivative of colebrook equation
double xx;
xx = 1+2.0/(eod/3.7+2.51/Re*X)/log(10.0)*2.51/Re;
return  xx;
}

double fcol(double Re,double eod)
{
//solution of the colebrook equation
// by using newton method
double fi=f(Re,eod);
double x=1.0/pow(fi,0.5);
int nmax=50;
double tolerance=1.0e-10;
for(int i=0;i<nmax;i++)
{
double fx1=fx(x,Re,eod);
x-=fx1/dfx(x,Re,eod);
if(fabs(fx1)<tolerance) return 1.0/(x*x);
}
return 1.0/(x*x);
}

//==============================================
// implementation of constructors functions
// for class pipe
pipe::pipe(const str & name,double le=1.0)
{
//this function reads pipe configuration
//from file pipe.dat
l=le;
ifstream cfin("pipe.dat",ios::in);
int ierror=1;
while(cfin>>pipeName)
{
if(pipeName==name) {ierror=0;break;}
}
if(ierror==1)
{cout<<"error required pipe name cannot be found in the data file\n";return;}
cfin>>shape;
cfin>>gasName;
//gmix temp(gasName);
if(!(g=new gmix(gasName)))
{
cout<<"Insufficient memory for gmix\n";
exit(1);
}
double pi=4.0*atan(1.0);
if(shape=="circular")
{ cfin>>d1;d2=0.0;A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
{
cfin>>d1>>d2;A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
{cfin>>d1>>d2;dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
{cfin>>d1>>d2;dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
{cfin>>d1>>d2;dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
{cfin>>d1>>d2;dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
}
eps=0.03e-3;
cfin>>mass_flow_rate;
P=1.0;;
KL=0.0;
cfin.close();
}

void pipe::change_d1(double dim1)
{
d1=dim1;
double pi=4.0*atan(1.0);
if(shape=="circular")
{ A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
{dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
{dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
{
dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}

}

void pipe::change_d2(double dim2)
{
d2=dim2;
double pi=4.0*atan(1.0);
if(shape=="circular")
{ d2=0;A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
{
 A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
{dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
{dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
{dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
{
dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
}

void pipe::change_pipe(const str & name,double le=1.0)
{
//this function reads pipe configuration
//from file pipe.dat
l=le;
ifstream cfin("pipe.dat",ios::in);
int ierror=1;
while(cfin>>pipeName)
{
if(pipeName==name) {ierror=0;break;}
}
if(ierror==1)
{cout<<"error required pipe name cannot be found in the data file\n";return;}
cfin>>shape;
cfin>>gasName;
if(!(g=new gmix(gasName)))
{
cout<<"Insufficient memory for gmix\n";
exit(1);
}
double    pi= 4.0*atan(1.0);
if(shape=="circular")
//d1 is diameter of the pipe
{ cfin>>d1;d2=0;A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
// d1 and d2 are two sides of the rectangle
{cfin>>d1>>d2;A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{cfin>>d1>>d2;dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
// d1 is diameter and d2 is the angle of the circular segment
{cfin>>d1>>d2;dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
//d1 is the diameter, d2 is the height of the circular section
{cfin>>d1>>d2;dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
// d1 and d2 two sides of right triangle
{cfin>>d1>>d2;dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=0.03e-3;
cfin>>mass_flow_rate;
P=1.0;
KL=1.0;
cfin.close();
}

pipe::pipe()
{
//this function reads pipe configuration
//from file pipe.dat for pipe_1
ifstream cfin("pipe.dat",ios::in);
int ierror=1;
l=0.25;
str name="pipe_6";
while(cfin>>pipeName)
{
if(pipeName==name) {ierror=0;break;}
}
if(ierror==1)
{cout<<"error required pipe name cannot be found in the data file\n";return;}
cfin>>shape;
cfin>>gasName;
if(!(g=new gmix(gasName)))
{
cout<<"Insufficient memory for gmix\n";
exit(1);
}
double    pi= 4.0*atan(1.0);
if(shape=="circular")
//d1 is diameter of the pipe
{ cfin>>d1;d2=0.0;A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
// d1 and d2 are two sides of the rectangle
{cfin>>d1>>d2;A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{cfin>>d1>>d2;dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
// d1 is diameter and d2 is the angle of the circular segment
{cfin>>d1>>d2;dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
//d1 is the diameter, d2 is the height of the circular section
{cfin>>d1>>d2;dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
// d1 and d2 two sides of right triangle
{cfin>>d1>>d2;dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=0.03e-3;
cfin>>mass_flow_rate;
P=1.0;
KL=0.0;
cfin.close();
}

pipe::pipe(const str & pn,
	  const str & sh,
	  const str & gn,
	  double dim1,double dim2,
	  double m,double e=0.03e-3,
	  double pp=1.0,double k=0.0,
	  double length=1.0)
{
// direct assignment constructor
// this constructor can assign all
// the required valued directly without
// going to the data file
pipeName=pn;
l=length;
shape=sh;
if(!(g=new gmix(gn)))
{
cout<<"Insufficient memory for gmix\n";
exit(1);
}
double    pi= 4.0*atan(1.0);
d1=dim1;
d2=dim2;
if(shape=="circular")
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;d2=0;}
else if(shape=="rectangular")
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
// d1 and d2 two sides of right triangle
{
dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}

pipe::pipe(const str & pn,
	  const str & sh,
	  gmix & gmix1,
	  double dim1,double dim2,
	  double m,double e=0.03e-3,
	  double pp=1.0,double k=0.0,
	  double length=1.0)
{
// direct assignment constructor
// this constructor can assign all
// the required valued directly without
// going to the data file
pipeName=pn;
shape=sh;
//setup gas name
g = new gmix;
g = & gmix1;
double    pi= 4.0*atan(1.0);
d1=dim1;
d2=dim2;
l=length;
if(shape=="circular")
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;d2=0;}
else if(shape=="rectangular")
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
// d1 and d2 two sides of right triangle
{ dh=2.0*d1*sin(d2)/(1.0+sin(d2)+cos(d2)); A=pi*dh*dh/4.0;}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}

void pipe::change_gas(const str & str1)
{
delete &g;
if(!(g=new gmix(str1)))
{
cout<<"Insufficient memory for gmix\n";
}
}

void pipe::change_gas(gmix & gmix1)
{
delete &g;
g = & gmix1;
}


void pipe::change_pipe(const str & pn,
	  const str & sh,
	  const str & gn,
	  double dim1,double dim2,
	  double e,double m,
	  double pp=1,double k=0)
{
// change assignment
// this constructor can assign all
// the required valued directly without
// going to the data file
pipeName=pn;
shape=sh;
if(!(g=new gmix(gn)))
{
cout<<"Insufficient memory for gmix\n";
exit(1);
}
double    pi= 4.0*atan(1.0);
d1=dim1;
d2=dim2;
if(shape=="circular")
//d1 is diameter of the pipe
{ A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
// d1 and d2 are two sides of the rectangle
{ A=d1*d2;dh=2.0*d1*d2/(d1+d2);
}
else if(shape=="concentric_annulus")
// d1 and d2 are inrernal and external diameters
//of the concentric annulus
{ dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
// d1 is diameter and d2 is the angle of the circular segment
{ dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
//d1 is the diameter, d2 is the height of the circular section
{ dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
// d1 and d2 two sides of right triangle
{ dh=4.0*d1*d2/(d1+d2+pow((d1*d1+d2*d2),0.5));
A=d1*d2/2.0;
}
eps=e;
mass_flow_rate=m;
P=pp;
KL=k;
}

pipe::~pipe()
{
  delete g;
  delete & pipeName;
}

double pipe::get_d1()
//getting the first dimension
{
return d1;
}

double pipe::get_d2()
//getting the first dimension
{
return d2;
}

double pipe::velocity(double t)
{
// velocity of the pipe m/x
double d =g->M/g->v(t,P);
return mass_flow_rate/d/A;
}

double pipe::Re(double t)
{
// reynold's number
double d=g->M/g->v(t,1.0);
double vis=g->vis(t);
double vv=velocity(t);
return dh*vv*d/vis;
}

double pipe::f(double t)
{
// this function calculated friction coefficient 
// for both laminar and turnulent flow for various 
// cross sections
double Rey=Re(t); //Reynold's number
if(Rey<2000)
  {
  str c         = "circular";
  str r         = "rectangular";
  str co_an     = "concentric_annulus";
  str cir_seg   = "circular_segment";
  str cir_sec   = "circular_sector";
  str rig_tri   = "right_triangle";
  double C;
  if(shape==c)
     //circular cross section, laminar flow
     { C=64.0; }
  else if(shape==r)
     {
     // Laminar rectangular channel friction factor
       double x=d1/d2;
       if(x>1) x=1/x; 
       double a1[6];
       a1[0]=    .959923253176985E+02;
       a1[1]=   -.130364643234979E+03;
       a1[2]=    .188053177394400E+03;
       a1[3]=  -.166159826967448E+03 ;
       a1[4]=    .952739160606451E+02;
       a1[5]=  -.258949934514987E+02 ;
       C=a1[5];
       for(int i=4;i>=0;i--) { C=C*x+a1[i];}
     }

  else if(shape==co_an)
     {
      //concentric annulus laminar friction factor
      double x=d2/d2;     
      if(x>1) x=1/x;
      double a1[5];
      a1[0]=     .717068909524062E+02;
      a1[1]=     .932037998655195E+03;
      a1[2]=    -.947586946520906E+04;
      a1[3]=     .204346570871727E+05;
      a1[4]=    -.118665325115712E+05;
      C=a1[4];
      for(int i=4;i>=0;i--) {C=C*x+a1[i];}
     }
  else if(shape==cir_seg)
     {
     //circular segment laminar friction factor
      double a1[5];
      a1[0]=     .639999999999997E+02;
      a1[1]=    -.162338014624008E+01;
      a1[2]=     .162620455430327E+01;
      a1[3]=   -.827251633703104E+00 ;
      a1[4]=     .138590724447624E+00;
      C=a1[4];
      for(int i=3;i>=0;i--) {C=C*d2+a1[i];}
     }
  else if(shape==cir_sec)
     {
     //circular section laminar friction factor
     double a1[4];
      a1[0]=     .480000000000000E+02;
      a1[1]=     .227909935649730E+02;
      a1[2]=    -.134959915493607E+02;
      a1[3]=     .325095868081599E+01;
      C=a1[3];
      for(int i=2;i>=0;i--) {C=C*d2+a1[i];}
     }
  else if(shape==rig_tri)
     {
     //right triangle laminar friction factor
     double a1[5];           
      a1[0]=     .480001792918036E+02;
      a1[1]=     .127207559700323E+02;
      a1[2]=    -.108590518340935E+02;
      a1[4]=     .158423382029386E+01;
      a1[5]=     .115425033696128E+01;
      C=a1[5];
      for(int i=4;i>=0;i--) {C=C*d2+a1[i];}
     }
  return C/Rey;
  }
else
  {
  double eod=eps/dh;
  return fcol(Rey,eod);
  }
}

double pipe::lentry(double t)
{
//hydraulic entry length
double Rey=Re(t);
double le;
if(Rey<2000) le=0.06*dh*Rey;
else le=4.4*dh*pow(Rey,(1.0/6.0));
return le;
}

double pipe::dP(double t,double P=1)
{
// pipe pressure drop for the given length plus
// local KL factors
  double V=velocity(t);
  return (KL+f(t)*l/dh)*(g->M/g->v(t,P))*V*V/2.0*1e-5;
}

double pipe::dPl_e(double t,double P=1)
{
long double Q=mass_flow_rate*g->v(t,P)/g->M;
long double a=d1/2.0;
long double b=d2/2.0;
long double x=0;
long double pi=4.0*atan(1.0);
for(int i=1;i<8000;i+=2)
  {
  x+=tanhl(pi*b*i/d1)/(long double)(i*i*i*i*i);
  }
x=1.0-192.0*a/(pi*pi*pi*pi*pi*b)*x;
double xx=b*a*a*a/(6.0*g->vis(t))*x;
xx=Q/xx*l*1e-5;
return xx;
}

void pipe::change_length(double l1)
{ l=l1;       }

void pipe::change_mass_flow_rate(double mfr)
{ mass_flow_rate=mfr; }


void pipe::change_P(double Pin)
//change inlet pressure
{ P=Pin; }

void pipe::change_KL(double klin)
//change local pressure drop coefficient
{KL=klin; }

void pipe::change_e(double e)
{ eps=e; }

double pipe::hc(double t)
{
//convective heat transfer coefficient
double Rey=Re(t);
double Nu;
if(Rey<2000)
  {
  if(shape=="circular")
     { Nu=3.66; }
  else if(shape=="rectangular")
     {
     // Laminar rectangular channel friction factor
       double x=d1/d2;
       if(x>1) x=1/x; 
       double a1[7];
       a1[0]=    .754000000000580E+01;
       a1[1]=   -.193789080028494E+02;
       a1[2]=    .333861212302836E+02;
       a1[3]=  -.170496495070800E+02 ;
       a1[4]=   -.304226226912837E+02;
       a1[5]=    .495297442580787E+02;
       a1[6]=   -.206246852871550E+02;
       Nu=a1[6];
      for(int i=5;i>=0;i--) {Nu=Nu*x+a1[i];}
     }
  else if(shape=="concentric_annulus")
     {
      double x=d2/d2;     
      if(x>1) x=1/x;
      double a1[4];
      a1[0]=     .401698044199492E+01;
      a1[1]=     .921803331830132E+00;
      a1[2]=    -.307845192109755E+00;
      a1[3]=     .229092480046921E+00;
       Nu=a1[3];
      for(int i=2;i>=0;i--) {Nu=Nu*x+a1[i];}
     }
  else if(shape=="circular_segment")
     {
     Nu=3.66;
     }
  else if(shape=="circular_sector")
     {
     Nu=3.66;
     }
  else if(shape=="right_triangle")
     {
     Nu=2.47;
     }
  }
else
  {
  double fr  = f(t);
  double Rey = Re(t);
  double Prt = g->Pr(t);
  // Gnielinski equation
  // Gnielinski, V, Int. Chem Eng., 16, 359,1976
  // Valid for 2000 < Re <5e5, 0.5 < Pr < 2000
  Nu=fr/8.0*(Rey-1000.0)*Prt/
  (1+12.7*sqrt(fr/8.0)*(pow(Prt,(2.0/3.0))-1.0));
  }
  double kt  =  g->k(t);
//  double Prt = g.Pr(t);
  double result;
  result=Nu*kt/dh;
  return result;
}

double pipe::Nu(double t)
//Nusselt's Number
{
return hc(t)*dh/g->k(t);
}

void pipe::operator = (const pipe & p1)
{
// assign new values to the class pipe variable
l=p1.l;
pipeName=p1.pipeName;
shape=p1.shape;
if(!(g=new gmix(gasName)))
{
cout<<"Insufficient memory for class gmix\n";
exit(1);
}
double pi=4.0*atan(1.0);
if(shape=="circular")
{ d1=p1.d1;A=pi*d1*d1/4;dh=d1;}
else if(shape=="rectangular")
{
d1=p1.d1;d2=p1.d2;A=d1*d2;dh=4.0*d1*d2/(2.0*d1+2.0*d2);
}
else if(shape=="concentric_annulus")
{d1=p1.d1;d2=p1.d2;dh=abs(d2-d1);A=pi*dh*dh/4.0;}
else if(shape=="circular_segment")
{d1=p1.d1;d2=p1.d2;dh=d1*(1+sin(2.0*d2)/2.0/(pi-d2));A=pi*dh*dh/4.0;}
else if(shape=="circular_sector")
{d1=p1.d1;d2=p1.d2;dh=d2/(1+d2)*d1; A=pi*dh*dh/4.0;   }
else if(shape=="right_triangle")
{d1=p1.d1;d2=p1.d2;dh=2.0*d1*sin(d2)/(1.0+sin(d2)+cos(d2)); A=pi*dh*dh/4.0;}
eps=p1.eps;
mass_flow_rate=p1.mass_flow_rate;;
P=p1.P;
KL=p1.KL;

}


