#include <perms.h>
void Emas3(_imp_string *Command, _imp_string *Params, _imp_string *Integer, _imp_name Flag);
void Emas3itos(int *I, _imp_string *S);
double Dens[21 /*0:20*/][21 /*0:20*/];
void Define(int *Chan, _imp_string *File) {
  _imp_string Number;
  int Flag;
  Emas3itos(*Chan, Number);
  Emas3(_imp_str_literal("define"), _imp_join(_imp_str_literal("sto"), Number), _imp_join(_imp_str_literal(","), *File),
        Flag);
  return;
}
double Density(double S, double T) {
  void Emas3(_imp_string * Comm, _imp_string * Parms, int Flag);
  int Aa;
  int B;
  int C;
  int D;
  int I;
  int J;
  int Flag;
  double X;
  double Y;
  double P;
  double Q;
  double Dx;
  double Dy;
  double R;
  if (T > 180) return (1000);
  if (T < 20) return (1000);
  if (S <= 0) {
    R = 1.00386 @3 - (2.23115 @-1 * T) - (2.40169 @-3 * REXP(T, 2));
    return (R);
  } else {
    Dx = Dens[2][0] - Dens[1][0];
    Dy = Dens[0][2] - Dens[0][1];
    X = ((S - Dens[1][0]) / Dx) + 1;
    Y = ((T - Dens[0][1]) / Dy) + 1;
    Aa = Intpt(X);
    B = Intpt(X + 1);
    C = Intpt(Y);
    D = Intpt(Y + 1);
    P = (Dens[Aa][D] - Dens[Aa][C]) * (Y - C) + Dens[Aa][C];
    Q = (Dens[B][D] - Dens[B][C]) * (Y - C) + Dens[B][C];
    R = (Q - P) * (X - Aa) + P;
  }
  return (R);
}
double Visc(double S, double T) {
  double Aa;
  double A1;
  double A2;
  double Bb;
  double B1;
  double B2;
  double C1;
  double C2;
  double N20;
  double N;
  Aa = 1.37220;
  A1 = -0.001015;
  A2 = 0.000005;
  Bb = 0.000813;
  B1 = 0.006102;
  B2 = -0.000040;
  C1 = 0.001550;
  C2 = 0.0000093;
  N20 = 1.002 + C1 * S + C2 * (REXP(S, 2));
  N = ((T - 20) / (T + 109)) *
      (Aa * (1 + A1 * S + A2 * (REXP(S, 2))) + Bb * (T - 20) * (1 + B1 * S + B2 * (REXP(S, 2))));
  N = N20 / REXP(10, (N));
  N = N / 1000;
  return (N);
}
double Heatcap(double S, double T) {
  double Aa;
  double A1;
  double A2;
  double A3;
  double Bb;
  double B1;
  double B2;
  double B3;
  double Cc;
  double C1;
  double C2;
  double C3;
  double Dd;
  double D1;
  double D2;
  double D3;
  double Ta;
  double Cp;
  A1 = 5.328;
  A2 = -9.76 @-2;
  A3 = 4.04 @-4;
  B1 = -6.913 @-3;
  B2 = 7.351 @-4;
  B3 = -3.15 @-6;
  C1 = 9.6 @-6;
  C2 = -1.927 @-6;
  C3 = 8.23 @-9;
  D1 = 2.5 @-9;
  D2 = 1.666 @-9;
  D3 = -7.125 @-12;
  Ta = T + 273.15;
  Aa = A1 + A2 * S + A3 * (REXP(S, 2));
  Bb = B1 + B2 * S + B3 * (REXP(S, 2));
  Cc = C1 + C2 * S + C3 * (REXP(S, 2));
  Dd = D1 + D2 * S + D3 * (REXP(S, 2));
  Cp = Aa + Bb * Ta + Cc * (REXP(Ta, 2)) + Dd * (REXP(Ta, 3));
  return (Cp);
}
double Tcond(double S, double T) {
  double Lamc;
  double Tc;
  double X;
  double Y;
  double G;
  double P;
  double Q;
  double Ta;
  double K;
  Ta = T + 273.15;
  Tc = 647;
  Lamc = 240;
  X = S * 0.0002;
  Y = S * 0.03;
  G = 343.5 + (S * 0.37);
  P = Log(Lamc + X);
  Q = (2.3 - (G / Ta)) * (REXP((1 - (Ta / (Tc + Y))), (1 / 3)));
  K = Exp;
  K = K / 1000;
  return (K);
}
double Bpe(double S, double T) {
  double A;
  double Bb;
  double C;
  double D;
  double E;
  double F;
  double X;
  double Ta;
  double B;
  X = S / 1000;
  Ta = T + 273.15;
  A = X * (REXP(Ta, 2)) / 13832;
  Bb = 0.001373 * Ta;
  C = -0.00272 * Ta * (REXP(X, 0.5));
  D = 17.86 * X;
  E = -0.0152 * X * Ta * ((Ta - 225.9) / (Ta - 236));
  F = -(2583 * X * (1 - X)) / Ta;
  B = A * (1 + Bb + C + D + E + F);
  return (B);
}
double Recov(double Rec, double T) {
  double A[21 /*0:20*/][21 /*0:20*/];
  int Aa;
  int B;
  int C;
  int D;
  int I;
  int J;
  double X;
  double Y;
  double P;
  double Q;
  double Dx;
  double Dy;
  double Delg;
  Define(7, _imp_str_literal("general.delg"));
  Selectinput(7);
  for (I = 0; I <= 8; I++)
    for (J = 0; J <= 4; J++) Read(A[I][J]);
  Closestream;
  Dx = A[2][0] - A[1][0];
  Dy = A[0][2] - A[0][1];
  X = ((T - A[1][0]) / Dx) + 1;
  Y = ((Rec - A[0][1]) / Dy) + 1;
  Aa = Intpt(X);
  B = Intpt(X + 1);
  C = Intpt(Y);
  D = Intpt(Y + 1);
  P = (A[Aa][D] - A[Aa][C]) * (Y - C) + A[Aa][C];
  Q = (A[B][D] - A[B][C]) * (Y - C) + A[B][C];
  Delg = (Q - P) * (X - Aa) + P;
  return (Delg);
}
double Barcp(double T) {
  double A;
  double B;
  double C;
  double D;
  A = -1.28115;
  B = 9.76656 @-2;
  C = -9.90726 @-4;
  D = 3.44861 @-6;
  double Cp;
  Cp = A + B * T + C * REXP(T, 2) + D * REXP(T, 3);
  return (Cp);
}
double Volucv(double T) {
  double Barcp(double T);
  double R;
  R = 0.462;
  double Cv;
  double Cp;
  Cp = Barcp(T);
  Cv = Cp - R;
  return (Cv);
}
double Gamma(double T) {
  double Cp;
  double Cv;
  double Gamma;
  Cp = Barcp(T);
  Cv = Volucv(T);
  Gamma = Cp / Cv;
  return (Gamma);
}
void Sateq(double X, double Y) {
  double Bpe(double S, double T);
  double A;
  double B;
  double C;
  A = 7.9186968;
  B = 1636.909;
  C = 224.92;
  Y = (B / (A - 0.434294 * Log(X * 750.06))) - C + X - X;
  Y = Bpe(X, X) - X;
  return;
}
double Sattemp(double S, double P) {
  double Bpe(double S, double T);
  void Impdavden(void Ff(double X, double Y), int N, int *Ifault, double Eps, double X, double Y);
  void Sateq(double X, double Y);
  double X[4 /*1:4*/];
  double Y[4 /*1:4*/];
  double A;
  double B;
  double C;
  double Tol;
  int Flag;
  A = 7.9186968;
  B = 1636.909;
  C = 224.92;
  Tol = 0.000001;
  if (P <= 0) return (60);
  X[1] = (B / (A - 0.434294 * Log(750.06 * P))) - C;
  if (S <= 0) return (X[1]);
  X[2] = Bpe(S, X[1]);
  X[3] = P;
  X[4] = S;
  Impdavden(Sateq(), 2, Flag, Tol, X, Y);
  if (Flag == 1) Printstring(_imp_str_literal("Something wrong in DAVDEN"));
  return (X[1]);
}
double Vappress(double S, double T) {
  double Bpe(double S, double T);
  double A;
  double B;
  double C;
  double P;
  A = 7.9186968;
  B = 1636.909;
  C = 224.92;
  if (S > 0) C -= Bpe(S, T);
  P = REXP(10, (A - (B / (C + T))));
  P = P * 1.01325 / 760.15;
  return (P);
}
double Latheat(double S, double T) {
  double Vappress(double S, double T);
  double Bpe(double S, double T);
  double Density(double S, double T);
  double A;
  double B;
  double C;
  double R;
  double P;
  double Vg;
  double Vl;
  double L;
  double Ta;
  A = 18.233473;
  B = 3769.122;
  C = 224.902;
  R = 8.314;
  Ta = T + 273.15;
  if (S > 0) C -= Bpe(S, T);
  P = Vappress(S, T);
  Vg = R * Ta / (P * 1.8 @3);
  Vl = 1 / Density(S, T);
  L = -B * Ta * (Vl - Vg) * P * 1 @2 / (REXP((T + C), 2));
  return (L);
}
double Htcond(int I, double Tsat, double Tw, double L, double D, double Nf, double Hf) {
  double Latheat(double S, double T);
  double Visc(double S, double T);
  double Density(double S, double T);
  double Tcond(double S, double T);
  double H;
  double K;
  double Rho;
  double G;
  double Lam;
  double Tref;
  double Delt;
  double U;
  double Re;
  G = 9.812;
  Delt = Tsat - Tw;
  Tref = Tsat - (3 * Delt / 4);
  K = Tcond(0, Tref);
  Rho = Density(0, Tref);
  Lam = Latheat(0, Tref) * 1000;
  U = Visc(0, Tref);
  if (Delt <= 0) return (1);
  if (I == 1) {
    H = (REXP(K, 3)) * (REXP(Rho, 2)) * G * Lam / (Delt * L * U);
    H = REXP(H, 0.25);
    H = H * 0.943;
    Re = 4 * Delt * H * L / (Lam * U);
    if (Re > 2100) {
      H = REXP(Re, 0.4);
      H = H * 0.0076;
      H = H * REXP(((REXP(K, 3)) * (REXP(Rho, 2)) * G / (REXP(U, 2))), 0.333);
    }
    return (H);
  }
}
double Htevap(int I, double S, double Tsat, double Tw, double L, double D, double M, double Nf, double Hf) {
  double Heatcap(double S, double T);
  double Tcond(double S, double T);
  double Visc(double S, double T);
  double Density(double S, double T);
  double Nu;
  double Pr;
  double Re;
  double H;
  double K;
  double U;
  double G;
  double Rho;
  double Gamma;
  double Tref;
  double Delt;
  double Cp;
  if (I == 1) Gamma = M / D;
  if (I == 2) Gamma = (M / 3.142) * D;
  Delt = Tsat - Tw;
  Tref = Tsat - (3 * Delt / 4);
  G = 9.812;
  U = Visc(S, Tref);
  K = Tcond(S, Tref);
  Rho = Density(S, Tref);
  Cp = Heatcap(S, Tref);
  Re = 4 * Gamma / U;
  Pr = U * Cp / K;
  if (I <= 2) {
    Nu = 0.565 / (Pr + 5.47);
    Nu -= 0.110;
    Nu = Nu * (REXP(Re, 0.231));
    H = Nu * K * (REXP((G * (REXP(Rho, 2)) / (REXP(U, 2))), 0.333));
    return (H);
  }
}
double Orifice(double Ph, double Pl, double T) {
  double R;
  double Ta;
  double F;
  R = 461.89;
  Ta = T + 273.15;
  if (Ph >= Pl) {
    F = REXP((Mod(REXP(Ph, 2) - REXP(Pl, 2)) / (R * Ta)), 0.5);
    return (F);
  }
  if (Ph < Pl) {
    F = REXP(-(Mod(REXP(Pl, 2) - REXP(Ph, 2)) / (R * Ta)), 0.5);
    return (F);
  }
}
double Uval(double Hc, double He, double K, double X, double Ff1, double Ff2) {
  double U;
  U = (1 / Hc) + (1 / He) + (X / K) + Ff1 + Ff2;
  U = 1 / U;
  return (U);
}
void Volumes(double T, double A, double Ti, double Theta, double Dthetadt, double *Va, double *Vd, double *Dvadt,
             double *Dvddt) {
  double Arg;
  double Omega;
  double Epsi;
  double Depsdt;
  double Dwdt;
  double Hold1;
  double Hold2;
  int Aa;
  int B;
  Arg = 2 * Pi / Ti;
  Omega = A * Sin(Arg * T);
  Dwdt = A * Arg * Cos(Arg * T);
  Epsi = Omega - Theta;
  Depsdt = Dwdt - Dthetadt;
  if (Epsi >= 0) {
    Aa = 1;
    B = 0;
  }
  if (Epsi < 0) {
    Epsi = Mod(Epsi);
    Aa = 0;
    B = 1;
  }
  Hold1 = 20 * (4.5 * ((Pi / 2) - Epsi) + 0.1406 * (Tan(Epsi) - Epsi) - 2.25);
  Hold2 = 20 * (4.5 * ((Pi / 2) + Epsi) - 0.1406 * (Tan(Epsi) + Epsi) - 2.25);
  *Va = Aa * Hold1 + B * Hold2;
  *Vd = B * Hold1 + Aa * Hold2;
  Hold1 = 20 * Depsdt * (0.1406 * ((1 / (REXP(Cos(Epsi), 2))) - 1) - 4.5);
  Hold2 = 20 * Depsdt * (4.5 - 0.1406 * (1 + (1 / (REXP(Cos(Epsi), 2)))));
  *Dvadt = -Depsdt * 20 * 4.5;
  *Dvddt = -*Dvadt;
  return;
}
double Filmthick(int Flag, double M, double T1, double T2, double S, double D) {
  double Visc(double S, double T);
  double Density(double S, double T);
  double Htcond(int I, double Th, double Tl, double L, double D, double Nf, double Hf);
  double Latheat(double S, double T);
  double U;
  double Rho;
  double G;
  double Delta;
  double Gamma;
  double Re;
  double H;
  double Lat;
  G = 9.812;
  if (!Flag) {
    U = Visc(S, T2);
    Rho = Density(S, T2);
    Gamma = M / D;
    Re = 4 * Gamma / U;
    if (Re <= 2000) {
      Delta = REXP(((3 * Gamma * U) / (G * (REXP(Rho, 2)))), (1 / 3));
      return (Delta);
    } else {
      Delta = 0.304 * REXP(((REXP(Gamma, 1.75)) * (REXP(U, 0.25)) / (G * (REXP(Rho, 2)))), (1 / 3));
      return (Delta);
    }
  }
  Lat = Latheat(0, T2);
  Rho = Density(0, T2);
  U = Visc(0, T2);
  if (T2 >= T1) return (0.0005 @-3);
  H = Htcond(1, T1, T2, 2, D, 0, 0) / 1000;
  Gamma = H * 2 * Mod((T1 - T2)) / Lat;
  Delta = REXP(((3 * Gamma * U) / (G * (REXP(Rho, 2)))), (1 / 3));
  return (Delta);
}
double Valvechar(double Ph, double Pi, double A, double B) {
  double Open;
  double Delp;
  Delp = Ph - Pi;
  Open = 0.5 * (1 + (2 / Pi) * (Arctan(1, (A * Delp + B))));
  return (Open);
}
void Valvecons(double *A, double *B) {
  void Emas3prompt(_imp_string * Text);
  double Xf;
  double Yf;
  double Hold1;
  double Hold2;
  Newline();
  Printstring(_imp_str_literal("Parameters determining the N/R valve behaviour–"));
  Newline();
  Emas3prompt(_imp_str_literal("First give the pressure drop across the valve........."));
  Read(Xf);
  Emas3prompt(_imp_str_literal("And now the % valve opening at that dP................."));
  Read(Yf);
  Hold1 = (Pi / 2) * (0.02 * Yf - 1);
  Hold2 = (Pi / 2) * (2.222 @-4 * Yf - 1);
  *A = 0.5 * (Tan(Hold1) - Tan(Hold2));
  *A = *A / Xf;
  *B = 0.5 * (Tan(Hold1) + Tan(Hold2));
  return;
}
double Satv(double S, double Vg) {
  double Bpe(double S, double T);
  double A;
  double B;
  double Tsat;
  A = Log(Vg);
  Tsat = 114.74 - 26.991 * A;
  if (!S) return (Tsat);
  B = Bpe(S, Tsat);
  Tsat += B;
  return (Tsat);
}
void Ndnewton(void Evaluate(double X, double Y), double X, double E, int N, int *Kmax, int *Flag) {
  double Det;
  int I;
  int J;
  int K;
  int M;
  int H;
  double B[N];
  double Dx[N];
  double F[N];
  double A[N][N];
  void Solvelneq(double A, double F, int N, double *Det);
  for (I = 1; I <= N; I++) Dx[I] = 5 * E;
  for (K = 1; K <= *Kmax; K++) {
    *Flag = 1;
    Evaluate;
    for (M = 1; M <= N; M++) B[M] = -F[M];
    for (J = 1; J <= N; J++) {
      X += Dx[J];
      Evaluate;
      for (I = 1; I <= N; I++) A[I][J] = (F[I] + B[I]) / Dx[J];
      X -= Dx[J];
    }
    Solvelneq(A, B, N, Det);
    if (Mod(Det) < 10 @-8) {
      *Flag = 2;
      return;
    }
    for (H = 1; H <= N; H++) {
      X += B[H];
      if (Mod(B[H]) > E) *Flag = 0;
    }
    if (*Flag == 1) {
      *Kmax = K;
      return;
    }
  }
  return;
}