#include <perms.h>
int _imp_mainep(int _imp_argc, char **_imp_argv) {
  void Prompt(_imp_string S);
  void Define(_imp_string S) {
    void Emas3(_imp_string * Command, _imp_string * Params, int *Flag);
    int Flag;
    Emas3(_imp_str_literal("DEFINE"), S, Flag);
  }
  int Kmax;
  int N;
  int S;
  int I;
  int K;
  int J;
  int Flag;
  int U;
  int Count;
  int B;
  int P;
  double Ph;
  double Pl;
  double Th;
  double Tl;
  double Phi;
  double Pli;
  double Thi;
  double Tli;
  double Tfhi;
  double Tfli;
  double Ktheta;
  double Ksec;
  double Dsec;
  double Dtheta;
  double Z;
  double Y;
  double Theta;
  double Vh;
  double Dvh;
  double Vl;
  double Dvl;
  double Sec;
  double Old;
  double Gamma;
  double Htc;
  float Dndt;
  float Dvdt;
  float Aa;
  float Bb;
  float Cc;
  float A1;
  float A2;
  float A3;
  float A4;
  float A5;
  float A6;
  float Vho;
  float Vlo;
  float Dvho;
  float Dvlo;
  float A;
  float Psec;
  float Dmass;
  float Mass;
  float Dpower;
  float Power;
  float Tclose;
  float Salt;
  float Area;
  float Height;
  double X[100 /*1:100*/];
  double Oldx[100 /*1:100*/];
  double Delta[100 /*1:100*/];
  double E[100 /*1:100*/];
  double Dens[21 /*0:20*/][21 /*0:20*/];
  Aa = 7.9186968;
  Bb = 1636.909;
  Cc = 224.92;
  Gamma = 1.135;
  void Adiabat(float Dvh, float Dvl, float Vh, float Vl, float Phi, float Pli, float Thi, float Tli, double *Ph,
               double *Pl, double *Th, double *Tl) {
    float Aph;
    float Apl;
    float Gamma;
    Gamma = 1.135;
    Aph = -1 * Gamma * Phi * Dvh / Vh;
    Apl = -1 * Gamma * Pli * Dvl / Vl;
    *Th = ((Gamma - 1) * Thi * Aph / (Gamma * Phi)) + Thi;
    *Tl = ((Gamma - 1) * Tli * Apl / (Gamma * Pli)) + Tli;
    *Ph = Aph + Phi;
    *Pl = Apl + Pli;
  }
  void Physprop(float *S, double *T, double *R, double *Cp, double *Bpe, double A) {
    float Dx;
    float Dy;
    float X;
    float Y;
    float P;
    float Q;
    int I;
    int J;
    int Aa;
    int Bb;
    int Cc;
    int Dd;
    *T = *T - 273.15;
    Dx = A - A;
    Dy = A - A;
    X = ((*S - A) / Dx) + 1;
    Y = ((*T - A) / Dy) + 1;
    Aa = Intpt(X);
    Bb = Intpt(X + 1);
    Cc = Intpt(Y);
    Dd = Intpt(Y + 1);
    if (Aa > 8) Aa = 8;
    if (Aa < 0) Aa = 1;
    if (Bb > 8) Bb = 8;
    if (Bb < 0) Bb = 1;
    if (Cc > 17) Cc = 17;
    if (Cc < 0) Cc = 1;
    if (Dd > 17) Dd = 17;
    if (Dd < 0) Dd = 1;
    P = (A - A) * (Y - Cc) + A;
    Q = (A - A) * (Y - Cc) + A;
    *R = (Q - P) * (X - Aa) + P;
    *T = *T + 273.15;
    float A0;
    float A1;
    float A2;
    float A3;
    float B0;
    float B1;
    float B2;
    float B3;
    float C0;
    float C1;
    float C2;
    float C3;
    float D0;
    float D1;
    float D2;
    float D3;
    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;
    A0 = A1 + A2 * *S + A3 * REXP(*S, 2);
    B0 = B1 + B2 * *S + B3 * REXP(*S, 2);
    C0 = C1 + C2 * *S + C3 * REXP(*S, 2);
    D0 = D1 + D2 * *S + D3 * REXP(*S, 2);
    *Cp = A0 + B0 * *T + C0 * REXP(*T, 2) + D0 * REXP(*T, 3);
    *Cp = *Cp * 1000;
    float Conc;
    float Bpe1;
    float Bpe2;
    float Bpe3;
    float Bpe4;
    float Bpe5;
    float Bpe6;
    Conc = *S / 1000;
    Bpe1 = Conc * REXP(*T, 2) / 13832;
    Bpe2 = 0.001373 * *T;
    Bpe3 = -0.00272 * *T * REXP(Conc, 0.5);
    Bpe4 = 17.86 * Conc;
    Bpe5 = -0.0152 * Conc * *T * ((*T - 225.9) / (*T - 236));
    Bpe6 = -(2583 * Conc * (1 - Conc)) / *T;
    *Bpe = Bpe1 * (1 + Bpe2 + Bpe3 + Bpe4 + Bpe5 + Bpe6);
  }
  void Evaluate(double F, double X, float Phi, float Pli, float Thi, float Tli, float Tfhi, float Tfli, float Tclose,
                float Salt, float A) {
    double Aa;
    double Bb;
    double Cc;
    double A1;
    double A2;
    double A3;
    double Ca;
    double Ka;
    double Kb;
    double Cpw;
    double Cpb;
    double Gamma;
    double Rhow;
    double Rhob;
    double Lamw;
    double Lamb;
    double Hc;
    double He;
    double W;
    double Xw;
    double Xb;
    double R;
    double K;
    double Mr;
    double Mf;
    double Tmix;
    double Bpe;
    double Tfeed;
    double Conc;
    Aa = 7.9186968;
    Bb = 1636.909;
    Cc = 224.92;
    Xw = 0.3 @-3;
    Xb = 0.27 @-3;
    Lamw = 2258 @3;
    Lamb = 2258 @3;
    R = 461.89;
    Gamma = 1.135;
    K = 104;
    W = 0.85 @-3;
    Conc = Salt / 34.5;
    Ca = (Gamma - 1) / Gamma;
    Cpw = 4216;
    Rhow = 962;
    Physprop(Salt, X, Rhob, Cpb, Bpe, Dens);
    Mf = 12 * Salt / (Salt - 34.5);
    Mr = 0.085 * A - Mf;
    Tfeed = ((Conc - 1) * (X - Tclose) + (X - Tclose)) / Conc;
    Tmix = (X * Mr + Tfeed * Mf) / (Mr + Mf);
    F = (Ca * Thi * (X - Phi) / Phi) + Thi - X;
    F = (Ca * Tli * (X - Pli) / Pli) + Tli - X;
    F = A * Cpw * Rhow * X * (Tfhi - X) + X * A * (X - X) * Dsec + (K * A * Dsec * (X - X) / W);
    F = Cpb * (Mf + Mr) * Dsec * (Tfli - Tmix) + Cpb * A * Xb * Rhob * (X - Tfli) + 6911 * A * (X - X) * Dsec +
        (K * A * Dsec * (X - X) / W);
    F = (-1 * Gamma * Phi * Dvh / Vh) - (X * A * (X - X) * Dsec * R * X / (Vh * Lamw)) + Phi - X;
    F = (-1 * Gamma * Pli * Dvl / Vl) + (6911 * A * (X - X) * Dsec * R * X / (Vl * Lamb)) + Pli - X;
    F = (Bb / (Aa - 0.434 * Log(X * 760 / 1.013 @5))) - Cc - X + 273.15;
    F = (Bb / (Aa - 0.434 * Log(X * 760 / 1.013 @5))) - Cc - X + 273.15 + Bpe;
    if (X < X) X = X;
    if (X < X) X = X;
    F = (REXP(X, 4)) * (X - X) - 1.26 @16;
    F = 1.2262 @-16 * X * (X - X) - (REXP(X, 3));
    return;
  }
  void Solver(double X, double E, int N, int *Kmax, int *Flag, float Phi, float Pli, float Thi, float Tli, float Tfhi,
              float Tfli) {
    double Det;
    int I;
    int J;
    int Iter;
    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] = 10 * E;
    for (Iter = 1; Iter <= *Kmax; Iter++) {
      *Flag = 1;
      Evaluate(F, X, Phi, Pli, Thi, Tli, Tfhi, Tfli, Tclose, Salt, Area);
      for (M = 1; M <= N; M++) B[M] = -F[M];
      for (J = 1; J <= N; J++) {
        X += Dx[J];
        Evaluate(F, X, Phi, Pli, Thi, Tli, Tfhi, Tfli, Tclose, Salt, Area);
        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) return;
    }
    return;
  }
  Prompt(_imp_str_literal("Initial pressure on h.p. side "));
  Read(Phi);
  Prompt(_imp_str_literal("Initial pressure on l.p. side "));
  Read(Pli);
  Prompt(_imp_str_literal("Initial temperature on h.p. side "));
  Read(Thi);
  Prompt(_imp_str_literal("Initial temperature on l.p. side "));
  Read(Tli);
  Prompt(_imp_str_literal("Initial water film temperature "));
  Read(Tfhi);
  Prompt(_imp_str_literal("Initial brine film temperature "));
  Read(Tfli);
  Prompt(_imp_str_literal("Temp. approach in H.Recov. "));
  Read(Tclose);
  Prompt(_imp_str_literal("Area for heat exchange "));
  Read(Area);
  Prompt(_imp_str_literal("Brine concentration factor "));
  Read(Salt);
  Salt = Salt * 34.5;
  Prompt(_imp_str_literal("Swept angle "));
  Read(Ktheta);
  Ktheta = Ktheta * Pi / 180;
  Prompt(_imp_str_literal("Period of oscillation "));
  Read(Ksec);
  Prompt(_imp_str_literal("Number of timesteps "));
  Read(S);
  Define(_imp_str_literal("1,theta"));
  Define(_imp_str_literal("2,temp"));
  Define(_imp_str_literal("3,delp"));
  Define(_imp_str_literal("4,time"));
  Define(_imp_str_literal("5,mass"));
  Define(_imp_str_literal("6,data"));
  Define(_imp_str_literal("10,emct12:general.densities"));
  Selectinput(10);
  for (I = 0; I <= 17; I++)
    for (J = 0; J <= 8; J++) Read(Dens[J][I]);
  Closestream;
  Dsec = Ksec / (S);
  N = 10;
  X[1] = Tli;
  X[2] = Thi;
  X[3] = Tfhi;
  X[4] = Tfli;
  X[5] = Pli;
  X[6] = Phi;
  X[9] = 10000;
  X[10] = 0.2 @-3;
  E[1] = 0.05;
  E[2] = 0.05;
  E[3] = 0.05;
  E[4] = 0.05;
  E[5] = 1000;
  E[6] = 1000;
  E[7] = 0.05;
  E[8] = 0.05;
  E[9] = 100;
  E[10] = 0.01 @-3;
  Count = 0;
  Old = -Ktheta / 2;
  Kmax = 500;
  U = 1;
  B = 1;
  Dmass = 0;
  Mass = 0;
  Dpower = 0;
  Power = 0;
  P = 0;
  K = 0;
  Psec = 0;
  A = Ktheta / Ksec;
  Y = 0;
  Vho = 20 * (4.5 * ((Pi / 2) - Old) + 0.125 * Old - 1.5);
  Vlo = 20 * (4.5 * ((Pi / 2) + Old) - 1.5 - 0.125 * Tan(Old));
  Dvho = 0;
  Dvlo = 0;
  for (I = 0; I <= S + 1; I++) {
    Sec = Dsec * I;
    if (Mod(Sec - (Ksec / 2)) < 0.01 * Dsec && B == 1) {
      Selectinput(0);
      Prompt(_imp_str_literal("New swept angle "));
      Read(Ktheta);
      Ktheta = Ktheta * Pi / 180;
      Ksec = (Ktheta / A);
      Dsec = Ksec / S;
    }
    if (I == 1 + S) {
      if (K > 4) break;
      K++;
      Thi = Tl;
      Tli = Th;
      Pli = Ph;
      Phi = Pl;
      Count = 0;
      Old = -Ktheta / 2;
      B = -B;
      U = 1;
      Vho = 20 * (4.5 * ((Pi / 2) - Old) + 0.125 * Old - 1.5);
      Vlo = 20 * (4.5 * ((Pi / 2) + Old) - 1.5 - 0.125 * Tan(Old));
      I = 0;
    }
    Sec = Dsec * I;
    Theta = -(Ktheta / 2) * Cos(3.142 * Sec / Ksec);
    Dtheta = Theta - Old;
    Old = Theta;
    X[7] = (Bb / (Aa - 0.434 * Log(X[5] * 7.502 @-3))) - Cc + 273;
    X[8] = (Bb / (Aa - 0.434 * Log(X[6] * 7.502 @-3))) - Cc + 273;
    for (J = 1; J <= 10; J++) Oldx[J] = X[J];
    if (U == 1) {
      Z = 1;
      Dvh = 0;
      Dvl = 0;
      Vh = 34;
      Vl = 34;
      A1 = X[5];
      A2 = X[6];
      A3 = X[1];
      A4 = X[2];
      A5 = X[3];
      A6 = X[4];
      Kmax = 500;
      Solver(X, E, N, Kmax, Flag, A1, A2, A3, A4, A5, A6);
      if (Flag == 2) {
        Newline();
        Printstring(_imp_str_literal("System is ill-conditioned."));
        exit(0);
      }
      if (!Flag) {
        Newline();
        Printstring(_imp_str_literal("System has failed to converge."));
        break;
      }
      if (Theta < 0) {
        Vho = 20 * (4.5 * ((Pi / 2) - Theta) + 0.125 * Theta - 1.5);
        Vlo = 20 * (4.5 * ((Pi / 2) + Theta) - 0.125 * Tan(Theta) - 1.5);
        Dvho = -87.5 * Dtheta;
        Dvlo = 20 * (4.5 - (0.125 / (REXP(Cos(Theta), 2)))) * Dtheta;
      }
      if (Theta >= 0) {
        Vho = 20 * (4.5 * ((Pi / 2) - Theta) + 0.125 * Tan(Theta) - 1.5);
        Vlo = 20 * (4.5 * ((Pi / 2) + Theta) - 0.125 * Theta - 1.5);
        Dvho = 20 * ((0.125 / (REXP(Cos(Theta), 2))) - 4.5) * Dtheta;
        Dvlo = 87.5 * Dtheta;
      }
      Adiabat(Dvho, Dvlo, Vho, Vlo, Phi, Pli, Thi, Tli, Ph, Pl, Th, Tl);
    }
    if (Count == 5) {
      Count = 0;
      Selectoutput(1);
      Print(Psec, 2, 2);
      Space();
      Print(Theta, 1, 2);
      Space();
      Print(Th, 3, 2);
      Space();
      Print(Tl, 3, 2);
      Space();
      Print(X[3], 3, 2);
      Space();
      Print(X[4], 3, 2);
      Space();
      Print(Ph / 1 @5, 1, 3);
      Space();
      Print(Pl / 1 @5, 1, 3);
      Space();
      Print(X[5] / 1 @5, 1, 3);
      Space();
      Print(X[6] / 1 @5, 1, 3);
      Newline();
      Selectoutput(2);
      Newline();
      Print(X[9], 4, 2);
      Space();
      Print(X[10] * 1000, 1, 3);
      Selectoutput(3);
      Newline();
      Print(Psec, 1, 2);
      Space();
      Print(B * (Ph - Pl) / 1 @5, 1, 3);
      Selectoutput(4);
      Newline();
    }
    if (X[5] + 0.01 @5 < Ph && X[6] > Pl + 0.01 @5) U = 0;
    for (J = 1; J <= 10; J++) {
      Delta[J] = X[J] - Oldx[J];
      X[J] = X[J] + 0.75 * Delta[J];
    }
    Phi = Ph;
    Pli = Pl;
    Thi = Th;
    Tli = Tl;
    Tfhi = X[3];
    Tfli = X[4];
    if (!U) {
      if (Theta < 0) {
        Dvh = -87.5 * Dtheta;
        Dvl = 20 * (4.5 - (0.125 / (REXP(Cos(Theta), 2)))) * Dtheta;
      }
      if (Theta >= 0) {
        Dvh = 20 * ((0.125 / (Cos(Theta) * 2)) - 4.5) * Dtheta;
        Dvl = 87.5 * Dtheta;
      }
      Kmax = 200;
      if (Z == 1) {
        Vh = Vho + 34;
        Vl = Vlo + 34;
        Kmax = 300;
        Thi = 0.5 * Thi + 0.5 * X[1];
        Tli = 0.5 * Tli + 0.5 * X[2];
        Phi = 0.5 * Phi + 0.5 * X[5];
        Pli = 0.5 * Pli + 0.5 * X[6];
        X[1] = Thi;
        X[2] = Tli;
        X[5] = Phi;
        X[6] = Pli;
        X[7] = (Bb / (Aa - 0.434 * Log(Phi * 7.502 @-3))) - Cc + 273;
        X[8] = (Bb / (Aa - 0.434 * Log(Pli * 7.502 @-3))) - Cc + 273;
        X[9] = 9000;
        X[10] = 0.2 @-3;
      }
      Z = 0;
      Solver(X, E, N, Kmax, Flag, Phi, Pli, Thi, Tli, Tfhi, Tfli);
      if (Flag == 2) {
        Newline();
        Printstring(_imp_str_literal("System is ill conditioned."));
        exit(0);
      }
      if (!Flag) {
        Newline();
        Printstring(_imp_str_literal("System has failed to converge."));
        break;
      }
      Ph = X[5];
      Pl = X[6];
      Th = X[1];
      Tl = X[2];
      Tfhi = X[3];
      Tfli = X[4];
    }
    if (Theta < 0) {
      Vh = 20 * (4.5 * ((Pi / 2) - Theta) + 0.125 * Theta - 1.5) + 34;
      Vl = 20 * (4.5 * ((Pi / 2) + Theta) - 0.125 * Tan(Theta) - 1.5) + 34;
    }
    if (Theta >= 0) {
      Vh = 20 * (4.5 * ((Pi / 2) - Theta) + 0.125 * Tan(Theta) - 1.5) + 34;
      Vl = 20 * (4.5 * ((Pi / 2) + Theta) - 0.125 * Theta - 1.5) + 34;
    }
    Count++;
    Psec += Dsec;
    P++;
    Dmass = X[9] * Area * (X[7] - X[3]) / 2258 @3;
    Dmass = REXP(Dmass, 2);
    Mass += Dmass;
    Dpower = B * (Ph - Pl) * 90 * Dtheta / Dsec;
    Dpower = REXP(Dpower, 2);
    Power += Dpower;
  }
  Mass = Mass / P;
  Mass = REXP(Mass, 0.5);
  Power = Power / P;
  Power = REXP(Power, 0.5);
  Selectoutput(6);
  Printstring(_imp_str_literal("Area for heat exchange ="));
  Print(Area, 4, 1);
  Newline();
  Printstring(_imp_str_literal("Temp. approach ="));
  Print(Tclose, 3, 2);
  Newline();
  Printstring(_imp_str_literal("RMS Fresh Water Flowrate ="));
  Print(Mass, 3, 3);
  Printstring(_imp_str_literal("kg/s"));
  Newline();
  Printstring(_imp_str_literal("RMS Power Consumption ="));
  Print(Power / 1000, 3, 3);
  Printstring(_imp_str_literal("kW"));
  Newline();
  Printstring(_imp_str_literal("RMS Specific Energy Usage ="));
  Print(Power / (Mass * 1000), 3, 3);
  Printstring(_imp_str_literal("kJ/kg"));
  Newline();
  Printstring(_imp_str_literal("Mean Wave Climate ="));
  Print(Power / 20 @3, 2, 2);
  Printstring(_imp_str_literal("kW/m"));
  Newline();
  Height = REXP((Power / (22 @3 * Ksec)), 0.5);
  Printstring(_imp_str_literal("RMS Wave Height ="));
  Print(Height, 2, 2);
  Printstring(_imp_str_literal("m"));
  Newline();
  Closestream;
  exit(0);
  return (1);
}