#include <perms.h>
int _imp_mainep(int _imp_argc, char **_imp_argv) {
  void Impdrunge(void F(double X, double Y, double Z), double H, double X0, double Y, double Res, int M, int N);
  void Impdavden(void F(double X, double Y), int N, int *Iflag, double Eps, double X, double Y);
  void Define(int *Chan, _imp_string *Char);
  void Emas3prompt(_imp_string * Text);
  void Valvecons(double *A, double *B);
  void Setup(double *Amp, double *Per, double *Sal, double *Tf, double *Mf, double *Area, double *Cd, double *Total);
  double Dens[21 /*0:20*/][21 /*0:20*/];
  double Ti;
  void Physprop(double Var) {
    double Density(double S, double T);
    double Gamma(double T);
    double Heatcap(double S, double T);
    double Latheat(double S, double T);
    double Tcond(double S, double T);
    double Volucv(double T);
    double Visc(double S, double T);
    double Satv(double S, double Vg);
    double Vg1;
    double Vg2;
    double Tsat1;
    double Tsat2;
    Vg1 = Var / Var;
    Vg2 = Var / Var;
    Tsat1 = Satv(0, Vg1);
    Tsat2 = Satv(Var, Vg2);
    Var = Latheat(0, Tsat1);
    Var = Latheat(Var, Tsat2);
    Var = Volucv(Tsat1);
    Var = Volucv(Tsat2);
    Var = Heatcap(0, Tsat1);
    Var = Heatcap(Var, Tsat2);
    Var = Density(0, Tsat1);
    Var = Density(Var, Tsat2);
    Var = Gamma(Var);
    Var = Gamma(Var);
    Var = Gamma(Var);
    Var = Gamma(Var);
    Var = Visc(0, Tsat1);
    Var = Visc(Var, Tsat2);
    Var = Tcond(0, Tsat1);
    Var = Tcond(Var, Tsat2);
    return;
  }
  void Recover(double Var) {
    void Ndnewton(void F(double X, double Y), double X, double E, int N, int *Kmax, int *Flag);
    double X[17 /*1:17*/];
    double Err[17 /*1:17*/];
    double S;
    int Kmax;
    int Flag;
    int I;
    void Bal(double X, double Y) {
      double Del[5 /*1:5*/];
      double Cp;
      double U1;
      double U2;
      double A1;
      double A2;
      double Lm1;
      double Lm2;
      double Lm3;
      double S;
      S = Var / 35;
      Del[1] = X - X;
      Del[2] = X - X;
      Del[3] = X - X;
      Del[4] = X - X;
      Del[5] = X - X;
      if (Del[1] * Del[2] <= 0)
        Lm1 = (Del[1] + Del[2]) / 2;
      else
        Lm1 = (Del[1] - Del[2]) / Log(Mod(Del[1] / Del[2]));
      if (Del[2] * Del[3] <= 0)
        Lm2 = (Del[2] + Del[3]) / 2;
      else
        Lm2 = (Del[2] - Del[3]) / Log(Mod(Del[2] / Del[3]));
      if (Del[4] * Del[5] <= 0)
        Lm3 = (Del[4] + Del[5]) / 2;
      else
        Lm3 = (Del[4] - Del[5]) / Log(Mod(Del[4] / Del[5]));
      Cp = 4.186 @3;
      U1 = 3000;
      U2 = 3000;
      A1 = 0;
      A2 = 1000;
      Y = Cp * X * (X - X) - X;
      Y = Cp * X * (X - X) - X;
      Y = U1 * A1 * Lm1 - X;
      Y = Cp * X * (X - X) - X;
      Y = 0.5 * U2 * A2 * Lm2 - X;
      Y = Cp * X * (X - X) - X;
      Y = 0.5 * U2 * A2 * Lm3 - X;
      Y = Cp * X * (X - X) - X - X;
      Y = X - Lm1;
      Y = X - Lm2;
      Y = X - Lm3;
      Y = Var - X;
      Y = (S * Var / (S - 1)) - X;
      Y = X + Var - X;
      Y = Var - X;
      Y = 20 - X;
      Y = Var - X;
      return;
    }
    S = Var / 35;
    X[1] = Var;
    X[2] = 100;
    X[3] = 24;
    X[4] = 20;
    X[5] = 95;
    X[6] = Var;
    X[7] = Var;
    X[8] = Var;
    X[9] = Var;
    X[10] = S * X[9] / (S - 1);
    X[11] = X[10] - X[9];
    X[12] = 5;
    X[13] = 5;
    X[14] = 5;
    X[15] = 150 @3;
    X[16] = 6 @6;
    X[17] = 6 @6;
    for (I = 1; I <= 14; I++) Err[I] = 0.01;
    Err[3] = 0.001;
    Err[15] = 15;
    Err[16] = 6 @2;
    Err[17] = 6 @2;
    Kmax = 200;
    Flag = 0;
    Ndnewton(Bal(), X, Err, 17, Kmax, Flag);
    if (!Flag) {
      Selectoutput(0);
      Printstring(_imp_str_literal("NDNEWTON fails...unconverged after "));
      Print(Kmax, 3, 1);
      Printstring(_imp_str_literal(" iterations."));
      Newline();
      exit(0);
    }
    if (Flag == 2) {
      Selectoutput(0);
      Printstring(_imp_str_literal("NDNEWTON fails...indeterminate problem."));
      Newline();
      exit(0);
    }
    Var = X[6];
    Var = X[3];
    Var = X[8];
    Var = X[10];
    return;
  }
  void Algebraic(double Var) {
    double Filmthick(int I, double M, double Tsat, double Tw, double S, double D);
    void Equil(double X, double Y) {
      double Satv(double S, double V);
      double Vg1;
      double Vg2;
      double Ts1;
      double Ts2;
      double M1;
      double M2;
      Vg1 = Var / Var;
      Vg2 = Var / Var;
      Ts1 = Satv(0, Vg1);
      Ts2 = Satv(Var, Vg2);
      M1 = Var + Var;
      M2 = Var + Var;
      Y = Var * Var * (Var - Ts1) + Var * Var * (Var - Ts1);
      Y = ((Y + Var * M1 * (Ts1 - X)) / Var) + Var - X;
      Y = Satv(0, Var / X) - X;
      Y = Var * Var * (Var - Ts2) + Var * Var * (Var - Ts2);
      Y = ((Y + Var * M2 * (Ts2 - X)) / Var) + Var - X;
      Y = Satv(Var, Var / X) - X;
    }
    double X[12 /*1:12*/];
    double Res[12 /*1:12*/];
    double Tw;
    double Alpha;
    int I;
    int J;
    X[1] = Var;
    X[2] = Var;
    X[3] = Var;
    X[4] = Var;
    Impdavden(Equil(), 4, I, 0.05, X, Res);
    if (I == 1) {
      Selectoutput(0);
      Printstring(_imp_str_literal("DAVDEN fails"));
      Newline();
      for (J = 1; J <= 12; J++) Print(Res[J], 1, 3);
      Newline();
      Closestream;
    }
    Var = X[2];
    Var = X[2];
    Var = X[1];
    Var = X[4];
    Var = X[4];
    Var = X[3];
    Tw = 0.5 * (Var + Var);
    Var = Filmthick(1, 0, Var, Tw, 0, 2.5) * Var * Var;
    Var = Filmthick(0, 0.5, Var, Tw, Var, 2.5) * Var * Var;
    return;
  }
  void Mathmod(double Ti, double Var, double Deriv) {
    void Volumes(double Ti, double A, double T, double Theta, double Dthetadt, double *Va, double *Vd, double *Dvadt,
                 double *Dvddt);
    double Orifice(double Ph, double Pl, double T);
    double Valvechar(double Ph, double Pl, double A, double B);
    double Satv(double S, double Vg);
    double Uval(double Hc, double He, double K, double X, double Ff1, double Ff2);
    double Htcond(int I, double Tsat, double Tw, double L, double D, double N, double H);
    double Htevap(int I, double S, double Tsat, double Tw, double L, double D, double M, double N, double H);
    double Vf[4 /*1:4*/];
    double R;
    double K;
    double Va;
    double Vb;
    double Vc;
    double Vd;
    double Dvadt;
    double Dvddt;
    double Vgb;
    double Vgc;
    double Pa;
    double Pb;
    double Pc;
    double Pd;
    double Cdao;
    double Tsatb;
    double Tsatc;
    double Twb;
    double Twc;
    double Xb;
    double Xc;
    double Xw;
    double Dnadt;
    double Dnbdt;
    double Dncdt;
    double Dnddt;
    double Phi;
    double Eps;
    double Hb;
    double Hc;
    double U;
    double Spring;
    double Inertia;
    double Kf;
    double G;
    double Ta;
    double Tb;
    double Tc;
    double Td;
    int I;
    R = 461.889;
    K = 0.392;
    Xw = 0.5 @-3;
    G = 9.812;
    Kf = 60;
    Spring = 3.352 @6;
    Inertia = 1.2906 @6;
    Ta = Var + 273.15;
    Tb = Var + 273.15;
    Tc = Var + 273.15;
    Td = Var + 273.15;
    Vb = Var;
    Vc = Var;
    Volumes(Ti, Var, Var, Var, Var, Va, Vd, Dvadt, Dvddt);
    Pa = Var * Ta * R / Va;
    Pb = Var * Tb * R / Vb;
    Pc = Var * Tc * R / Vc;
    Pd = Var * Td * R / Vd;
    Cdao = Var;
    Vf[1] = Cdao * Orifice(Pa, Pb, Var) * Valvechar(Pa / 1 @5, Pb / 1 @5, Var, Var);
    Vf[2] = Cdao * Orifice(Pd, Pb, Var) * Valvechar(Pd / 1 @5, Pb / 1 @5, Var, Var);
    Vf[3] = Cdao * Orifice(Pc, Pa, Var) * Valvechar(Pc / 1 @5, Pa / 1 @5, Var, Var);
    Vf[4] = Cdao * Orifice(Pc, Pd, Var) * Valvechar(Pc / 1 @5, Pd / 1 @5, Var, Var);
    Vgb = Var / Var;
    Vgc = Var / Var;
    Tsatb = Satv(0, Vgb);
    Tsatc = Satv(Var, Vgc);
    Twb = 0.25 * (3 * Var + Var);
    Twc = 0.25 * (Var + 3 * Var);
    Hb = Htcond(1, Tsatb, Twb, 3, 2.5, 0, 0);
    Hc = Htevap(1, Var, Tsatc, Twc, 3, 2.5, 0.5, 0, 0);
    Hb = Hb / 1000;
    Hc = Hc / 1000;
    U = Uval(Hb, Hc, K, Xw, 0, 0);
    Twb = Var - (U * (Var - Var) / Hb);
    Twc = Var + (U * (Var - Var) / Hc);
    Phi = Hb * Var * (Var - Twb) / Var;
    Eps = Hc * Var * (Twc - Var) / Var;
    Dnadt = Vf[3] - Vf[1];
    Dnddt = Vf[4] - Vf[2];
    Dnbdt = Vf[1] + Vf[2] - Phi;
    Dncdt = Eps - Vf[3] - Vf[4];
    Deriv = Var;
    Deriv = (Kf * (Pa - Pd) - Spring * Sin(Var)) / Inertia;
    Deriv = (Var - 1) * Ta * (-Dvadt / Va) + (Dnadt / Var);
    Deriv += Vf[3] * (Var - Var) / Var;
    Deriv = (Var - 1) * Tb * (Dnbdt / Var);
    Deriv = Deriv + (Vf[1] * (Var - Var) / Var) + (Vf[2] * (Var - Var) / Var);
    Deriv = (Var - 1) * Tc * (Dncdt / Var);
    Deriv = (Var - 1) * Td * ((-Dvddt / Vd) + (Dnddt / Var));
    Deriv += Vf[4] * (Var - Var) / Var;
    Deriv = Dnadt;
    Deriv = Dnbdt;
    Deriv = Dncdt;
    Deriv = Dnddt;
    Deriv = (Var * Phi - (K * Var * (Twb - Twc) / Xw)) / (Var * Var);
    Deriv = (-(Var * Eps - Var * Var * (Var - Var)));
    Deriv = (Deriv + (K * Var * (Twb - Twc) / Xw)) / (Var * Var);
    for (I = 13; I <= 44; I++) Deriv = 0;
    return;
  }
  double Htcond(int I, double Tsat, double Tw, double L, double D, double N, double H);
  double Y[50 /*1:50*/];
  double Result[44 /*1:44*/][2 /*0:1*/];
  double H;
  double T;
  double Total;
  double Tout;
  double Mass1;
  double Mass2;
  double Tw;
  double Hc;
  int M;
  int N;
  int I;
  int J;
  int Count;
  int Step1;
  int Step2;
  int Max;
  Tout = 0;
  Define(5, _imp_str_literal("emct12:general.densities"));
  Selectinput(5);
  for (I = 0; I <= 17; I++)
    for (J = 0; J <= 8; J++) Read(Dens[J][I]);
  Selectinput(0);
  Closestream;
  Setup(Y[18], Y[19], Y[20], Y[21], Y[24], Y[22], Y[23], Total);
  Define(9, _imp_str_literal("cons"));
  Selectoutput(9);
  for (I = 18; I <= 24; I++) {
    Print(Y[I], 0, 3);
    Space();
  }
  Selectoutput(0);
  Closestream;
  Valvecons(Y[41], Y[42]);
  H = 0.02;
  M = 44;
  N = 1;
  Y[1] = 0.00;
  Y[2] = 0.00;
  Y[3] = 102;
  Y[4] = 102;
  Y[5] = 102;
  Y[6] = 102;
  Y[7] = 56.1;
  Y[8] = 26.14;
  Y[9] = 26.14;
  Y[10] = 56.1;
  Y[11] = 102;
  Y[12] = 102;
  Y[13] = 35;
  Y[14] = 700;
  Y[39] = (90 - Y[22] * 0.5 @-3) / 2;
  Y[40] = Y[39];
  Y[15] = 0;
  Y[16] = 0;
  Y[17] = 0;
  Y[46] = 25;
  Y[47] = 23;
  Max = Int(Total / H);
  Ti = 0.00;
  Step1 = 0;
  Step2 = 0;
  Mass1 = 0;
  Define(10, _imp_str_literal("result"));
  Selectoutput(10);
  for (Count = 1; Count <= Max; Count++) {
    Physprop(Y);
    Algebraic(Y);
    Tw = 0.25 * (3 * Y[11] + Y[12]);
    Hc = Htcond(1, Y[4], Tw, 3, 2.5, 0, 0) / 1000;
    Mass2 = Hc * Y[22] * (Y[4] - Tw) / Y[25];
    Mass1 += Mass2;
    Y[45] = Mass1 / Count;
    if (Step2 == 10) {
      Recover(Y);
      Step2 = 0;
    }
    Tout += 4.186 * (Y[11] + 0.97 * Y[12] - 1.97 * Y[21]);
    if (Ti <= Total - 2 * Y[19]) Step1 = 0;
    if (Step1 == 5) {
      Print(Ti, 3, 3);
      Space();
      for (I = 1; I <= 14; I++) {
        Print(Y[I], 3, 3);
        Space();
      }
      Newline();
      Step1 = 0;
    }
    Impdrunge(Mathmod(), H, Ti, Y, Result, M, 1);
    for (J = 1; J <= 44; J++) Y[J] = Result[J][1];
    Step1++;
    Step2++;
    Ti += H;
  }
  Print(101010.10, 6, 2);
  Newline();
  Selectoutput(0);
  Closestream;
  Newlines(3);
  Printstring(_imp_str_literal("Finished"));
  exit(0);
  return (1);
}