/*********************************************************************
  nernst.cpp
  Copyright (c) 2002 Haruhiko Okumura <okumura@matsusaka-u.ac.jp>
  Compile with gcc-3.x: g++ -O2 nernst.cpp -o nernst
**********************************************************************

以下はまだ作業中です。コンパイルするとエラーになります。

一つ前のバージョンは nernst-000814.cpp という名前に変えました。

h:
    stepsize [meter], i.e., distance between nearest-neighbor grid
    points.

i, j:
    Discretized (integer) x- and y-coordinates, respectively.
    Distance between points (i,j) and (i+1,j) is h.
    So is distance between (i,j) and (i,j+1).
    i runs from 0 to imax.
    j runs from 0 to jmax.

T[i][j]:
    Temperature [Kelvin] at point (i,j)

V[i][j]:
    Electrochemical potential [Volt] at point (i,j)

Jx[i][j]:
    x-component of electric current density [Ampere/m^2]
    on line segment (i,j)-(i+1,j)

Jy[i][j]:
    y-component of electric current density [Ampere/m^2]
    on line segment (i,j)-(i,j+1)

pos[i][j]:
    Flags for point (i,j) to indicate which neighbors exist
    and how V[i][j] and T[i][j] are constrained.

          1: (i-1,j) exists
          2: (i+1,j) exists
          4: (i,j-1) exists
          8: (i,j+1) exists
         16: V fixed
         32: T fixed
         64: belongs to shortcircuited island A
        128: belongs to shortcircuited island B

    pos[i][j] = 0 represents vacancy.

    For example, for an H-shaped sample with fixed temperatures on
    the left and right sides, the grid points have the following flags:

        j
        |  6+32  7  7  5  0  0  6  7  7  5+32
        | 14+32 15 15 13  0  0 14 15 15 13+32
        | 14+32 15 15 13  0  0 14 15 15 13+32
        | 14+32 15 15 15  7  7 15 15 15 13+32
        | 14+32 15 15 15 15 15 15 15 15 13+32
        | 14+32 15 15 15 15 15 15 15 15 13+32
        | 14+32 15 15 15 11 11 15 15 15 13+32
        | 14+32 15 15 13  0  0 14 15 15 13+32
        | 14+32 15 15 13  0  0 14 15 15 13+32
        | 10+32 11 11  9  0  0 10 11 11  9+32
        +------------------------------------i
          T=T0                           T=T1

*********************************************************************/

#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <ios>
#include <iomanip>
#include <cstdlib>
#include <cmath>
#include <csignal>

using namespace std;

const int NX = 17 * 32;         // up to 17mm in x-direction for h = (1/32)mm
const int NY = 12 * 32;         // up to 12mm in y-direction for h = (1/32)mm

const double BOUNDFAC = 0.5;    // 50% of a boundary point belongs
                                // to the sample

const int SHOWSTEP = 100;       // show results every 100 iterations

int imax;                       // i = 0..imax; imax <= NX
int jmax;                       // j = 0..jmax; jmax <= NY
int resolution;                 // steps / mm

int pos[NX+1][NY+1];

double T[NX+1][NY+1];
double V[NX+1][NY+1];
double Jx[NX+1][NY+1];
double Jy[NX+1][NY+1];
// double T0, T1;                  // temperatures at i=0 and i=imax [K]
// double V0, V1;                  // voltages at i=0 and i=imax [V]
// bool V_not_fixed;               // V is chosen so that net current = 0
double length;                  // length between i=0 and i=imax [m]
double width;                   // width begween j=0 and j=jmax [m]
double h;                       // stepsize [m]
double sor;                     // SOR factor (for T updates only)
double B;                       // magnetic field [Tesla]
double fac, weight1, weight2;   // see below
volatile bool quitting = false; // turns true if SIGINT raised

double rho0, drho0;             // rho(T=0), d(rho)/dT
double seebeck0, dseebeck0;     // seebeck(T=0), d(seebeck)/dT
double kappa0, dkappa0;         // kappa(T=0), d(kappa)/dT
double hall0, dhall0;           // hall(T=0), d(hall)/dT
double nernst0, dnernst0;       // nernst(T=0), d(nernst)/dT
double righi0, drighi0;         // righi(T=0), d(righi)/dT

inline double rho(double T) { return rho0 + drho0 * T; }
inline double seebeck(double T) { return seebeck0 + dseebeck0 * T; }
inline double kappa(double T) { return kappa0 + dkappa0 * T; }
inline double hall(double T) { return hall0 + dhall0 * T; }
inline double nernst(double T) { return nernst0 + dnernst0 * T; }
inline double righi(double T) { return righi0 + drighi0 * T; }
inline double drho(double T) { return drho0; }
inline double dseebeck(double T) { return dseebeck0; }
inline double dkappa(double T) { return dkappa0; }
inline double dhall(double T) { return dhall0; }
inline double dnernst(double T) { return dnernst0; }

class savefmt {
    ios_base& f;
    ios_base::fmtflags ff;
    streamsize prec;            // or 'int'
    int c;
public:
    savefmt(ios_base& _f) : f(_f) {
        ff = f.flags();
        prec = f.precision();
    }
    ~savefmt() {
        f.flags(ff);
        f.precision(prec);
    }
};

void error(int line)
{
    cerr << __FILE__ << ": " << line << endl;
    exit(1);
}

void myhandler(int signum)      // Called when Ctrl-C is hit
{
    signal(SIGINT, SIG_DFL);    // quit on another Ctrl-C
    cerr << "Signal " << signum << " caught." << endl;
    quitting = true;            // quit on next iteration
}

double update_T_sub(int i, int j) // update temperature at (i,j)
{
    if ((pos[i][j] & 32) != 0) return 0;

    double t = T[i][j];
    double beta = nernst(t) * B;
    double k = kappa(t);
    double tt, cx, cy, tx, ty;
    switch (pos[i][j] & 15) {
    case 15:
        cx = (Jx[i-1][j] + Jx[i][j]) / 2;
        cy = (Jy[i][j-1] + Jy[i][j]) / 2;
        tx = (T[i+1][j] - T[i-1][j]) / (2 * h);
        ty = (T[i][j+1] - T[i][j-1]) / (2 * h);
        tt = (T[i-1][j] + T[i][j-1] + T[i][j+1] + T[i+1][j]) / 4;
        break;
    case 11:
        cx = (Jx[i-1][j] + Jx[i][j]) / 2;
        tx = (T[i+1][j] - T[i-1][j]) / (2 * h);
        if ((pos[i][j] & 16) != 0) { // V fixed
            cy = Jy[i][j];
            ty = (T[i][j+1] - T[i][j]) / h;
        } else {
            cy = 0;
            ty = (beta * t * cx) / k + righi(t) * B * tx;
        }
        tt = (T[i-1][j] + T[i+1][j] + 2 * (T[i][j+1] - h * ty)) / 4;
        break;
    case 7:
        cx = (Jx[i-1][j] + Jx[i][j]) / 2;
        tx = (T[i+1][j] - T[i-1][j]) / (2 * h);
        if ((pos[i][j] & 16) != 0) { // V fixed
            cy = Jy[i][j-1];
            ty = (T[i][j] - T[i][j-1]) / h;
        } else {
            cy = 0;
            ty = (beta * t * cx) / k + righi(t) * B * tx;
        }
        tt = (T[i-1][j] + T[i+1][j] + 2 * (T[i][j-1] + h * ty)) / 4;
        break;
    case 14:
        cy = (Jy[i][j-1] + Jy[i][j]) / 2;
        ty = (T[i][j+1] - T[i][j-1]) / (2 * h);
        if ((pos[i][j] & 16) != 0) { // V fixed
            cx = Jx[i][j];
            tx = (T[i+1][j] - T[i][j]) / h;
        } else {
            cx = 0;
            tx = -((beta * t * cy) / k + righi(t) * B * ty);
        }
        tt = (T[i][j-1] + T[i][j+1] + 2 * (T[i+1][j] - h * tx)) / 4;
        break;
    case 13:
        cy = (Jy[i][j-1] + Jy[i][j]) / 2;
        ty = (T[i][j+1] - T[i][j-1]) / (2 * h);
        if ((pos[i][j] & 16) != 0) { // V fixed
            cx = Jx[i-1][j];
            tx = (T[i][j] - T[i-1][j]) / h;
        } else {
            cx = 0;
            tx = -((beta * t * cy) / k + righi(t) * B * ty);
        }
        tt = (T[i][j-1] + T[i][j+1] + 2 * (T[i-1][j] + h * tx)) / 4;
        break;
    case 5:
        if ((pos[i][j] & 16) != 0) { // V fixed
            if ((pos[i][j-1] & 16) != 0) cx = Jx[i-1][j];  else cx = 0;
            if ((pos[i-1][j] & 16) != 0) cy = Jy[i][j-1];  else cy = 0;
        } else {
            cx = cy = 0;
        }
        tx = ty = 0;
        tt = (T[i-1][j] + T[i][j-1]) / 2;
        break;
    case 6:
        if ((pos[i][j] & 16) != 0) {
            if ((pos[i][j-1] & 16) != 0) cx = Jx[i][j];    else cx = 0;
            if ((pos[i+1][j] & 16) != 0) cy = Jy[i][j-1];  else cy = 0;
        } else {
            cx = cy = 0;
        }
        tx = ty = 0;
        tt = (T[i+1][j] + T[i][j-1]) / 2;
        break;
    case 9:
        if ((pos[i][j] & 16) != 0) {
            if ((pos[i][j+1] & 16) != 0) cx = Jx[i-1][j];  else cx = 0;
            if ((pos[i-1][j] & 16) != 0) cy = Jy[i][j];    else cy = 0;
        } else {
            cx = cy = 0;
        }
        tx = ty = 0;
        tt = (T[i-1][j] + T[i][j+1]) / 2;
        break;
    case 10:
        if ((pos[i][j] & 16) != 0) {
            if ((pos[i][j+1] & 16) != 0) cx = Jx[i][j];  else cx = 0;
            if ((pos[i+1][j] & 16) != 0) cy = Jy[i][j];  else cy = 0;
        } else {
            cx = cy = 0;
        }
        tx = ty = 0;
        tt = (T[i+1][j] + T[i][j+1]) / 2;
        break;
    default:
        tt = cx = cy = tx = ty = 0;
        error(__LINE__);
    }
    double c2 = cx * cx + cy * cy;
    double r = rho(t);
    double u = - r * c2;
    u += (t * dseebeck(t) + beta * t * B * dhall(t) / r)
        * (tx * cx + ty * cy);
    u += (beta * (2 - t * drho(t) / r) + t * dnernst(t) * B)
        * (ty * cx - tx * cy);
    u += (beta * t * B * dnernst(t) / r - dkappa(t)) * (tx * tx + ty * ty);
    u /= k - beta * beta * t / r; // Laplacian T
    tt = tt - 0.25 * h * h * u;
    T[i][j] = tt + sor * (tt - t);
    return fabs(tt - t);
}

double update_T()               // update all temperatures
{
    double emax = 0;
    for (int i = 0; i <= imax; i++) {
        for (int j = 0; j <= jmax; j++) {
            if (pos[i][j] == 0) continue;
            double e = update_T_sub(i, j);
            if (e > emax) emax = e;
        }
    }
    return emax;
}

bool diverging()
{
    for (int i = 0; i <= imax; i++) {
        for (int j = 0; j <= jmax; j++) {
            if (pos[i][j] == 0) continue;
            if (T[i][j] < 0 || T[i][j] > 1000 || T[i][j] != T[i][j])
                return true;
        }
    }
    return false;
}

// update x-component of electric current density
void calculate_Jx_sub(int i, int j)
{
    int p = pos[i][j];
    int q = pos[i+1][j];
    double t = (T[i][j] + T[i+1][j]) / 2;
    double alpha = seebeck(t);
    double beta = nernst(t) * B;
    double r = rho(t);
    double s = hall(t) * B;

    if ((p & q & (4|8)) == (4|8)) {
        // [i][j-1], [i][j+1], [i+1][j-1], [i+1][j+1] all exist
        double tx = T[i+1][j] - T[i][j];
        double ty = (T[i][j+1] + T[i+1][j+1] - T[i][j-1] - T[i+1][j-1]) / 4;
        double vx = V[i+1][j] - V[i][j];
        double vy = (V[i][j+1] + V[i+1][j+1] - V[i][j-1] - V[i+1][j-1]) / 4;
        double ax = - vx - alpha * tx + beta * ty;
        double ay = - vy - alpha * ty - beta * tx;
        double denom = h * (r * r + s * s);
        Jx[i][j] = (r * ax + s * ay) / denom;
        weight1 = r / denom;
        weight2 = - weight1;
        fac = 1;
    } else {
        // Jy = 0 (and in this case Qy = 0 also)
        // Note we don't consider cases with electrodes or heat sinks on
        // sides other than i=0 and i=imax.
        double k = kappa(t);
        double a = alpha - beta * righi(t) * B;
        double z = r - beta * beta * t / k;
        Jx[i][j] = (V[i][j] - V[i+1][j] + a * (T[i][j] - T[i+1][j])) / (h * z);
        weight1 = BOUNDFAC / (h * z);
        weight2 = - weight1;
        fac = BOUNDFAC;
    }
}

// update y-component of electric current density
void calculate_Jy_sub(int i, int j)
{
    double t = (T[i][j] + T[i][j+1]) / 2;
    double alpha = seebeck(t);
    double beta = nernst(t) * B;
    double r = rho(t);
    double s = hall(t) * B;
    int p = pos[i][j];
    int q = pos[i][j+1];

    if ((p & q & 3) == 3) {
        // [i+1][j], [i-1][j], [i+1][j+1], [i-1][j+1] all exist
        double tx = (T[i+1][j] + T[i+1][j+1] - T[i-1][j] - T[i-1][j+1]) / 4;
        double ty = T[i][j+1] - T[i][j];
        double vx = (V[i+1][j] + V[i+1][j+1] - V[i-1][j] - V[i-1][j+1]) / 4;
        double vy = V[i][j+1] - V[i][j];
        double ax = - vx - alpha * tx + beta * ty;
        double ay = - vy - alpha * ty - beta * tx;
        double denom = h * (r * r + s * s);
        Jy[i][j] = (r * ay - s * ax) / denom;
        weight1 = r / denom;
        weight2 = - weight1;
        fac = 1;
    } else if ((p & (16|32)) == 0 || (q & (16|32)) == 0) {
        // neither V nor T fixed at one or both of the ends
        // so that Jx = Qx = 0 on the edge
        double k = kappa(t);
        double a = alpha - beta * righi(t) * B;
        double z = r - beta * beta * t / k;
        Jy[i][j] =
            (V[i][j] - V[i][j+1] + a * (T[i][j] - T[i][j+1])) / (h * z);
        weight1 = BOUNDFAC / (h * z);
        weight2 = - weight1;
        fac = BOUNDFAC;
    } else if ((p & q & 16) != 0) {
        // V fixed at both ends, so that Jx != 0
        double tx, vx;
        if ((p & 3) == 2 && (q & 3) == 2) {
            tx = (T[i+1][j] + T[i+1][j+1] - T[i][j] - T[i][j+1]) / 2;
            vx = (V[i+1][j] + V[i+1][j+1] - V[i][j] - V[i][j+1]) / 2;
        } else if ((p & 3) == 1 && (q & 3) == 1) {
            tx = (T[i][j] + T[i][j+1] - T[i-1][j] - T[i-1][j+1]) / 2;
            vx = (V[i][j] + V[i][j+1] - V[i-1][j] - V[i-1][j+1]) / 2;
        } else {
            tx = vx = 0;
            error(__LINE__);
        }
        double ty = T[i][j+1] - T[i][j];
        double vy = V[i][j+1] - V[i][j];
        double ax = - vx - alpha * tx + beta * ty;
        double ay = - vy - alpha * ty - beta * tx;
        double denom = h * (r * r + s * s);
        Jy[i][j] = (r * ay - s * ax) / denom;
        weight1 = BOUNDFAC * r / denom;
        weight2 = - weight1;
        fac = BOUNDFAC;
    } else if ((p & 16) == 0 || (q & 16) == 0) {
        // V not fixed at one or both of the ends, so that Jx = 0
        double tx;
        if ((p & q & 2) != 0) {
            tx = (T[i+1][j] + T[i+1][j+1] - T[i][j] - T[i][j+1]) / 2;
        } else if ((p & q & 1) != 0) {
            tx = (T[i][j] + T[i][j+1] - T[i-1][j] - T[i-1][j+1]) / 2;
        } else {
            tx = 0;
            error(__LINE__);
        }
        double ty = T[i][j+1] - T[i][j];
        double vy = V[i][j+1] - V[i][j];
        double ay = - vy - alpha * ty - beta * tx;
        Jy[i][j] = ay / (h * r);
        // = (V[i][j] - V[i][j+1]) / (h * r) + ...
        weight1 = BOUNDFAC / (h * r);
        weight2 = - weight1;
        fac = BOUNDFAC;
    } else {
        error(__LINE__);
    }
}

void shortcircuit()
{
    for (int k = 2; k <= 3; k++) {
        int i1 = steps(circ[k].x - circ[k].r) - 1;
        if (i1 < 0) i1 = 0;
        int i2 = steps(circ[k].x + circ[k].r) + 1;
        if (i2 > imax) i2 = imax;
        int j1 = steps(circ[k].y - circ[k].r) - 1;
        if (j1 < 0) j1 = 0;
        int j2 = steps(circ[k].y + circ[k].r) + 1;
        if (j2 > jmax) j2 = jmax;
        for (int i = i1; i <= i2; i++) {
            for (int j = j1; j <= j2; j++) {
                if (circ[k].contains(i * h, j * h))
                    V[i][j] = circ[k].v;
            }
        }
    }
}

double update_J_V()           // update potentials and currents
{
    double emax = 0;
    double w2 = 0;
    double w3 = 0;
    double sum2 = 0;
    double sum3 = 0;
    double abssum2 = 0;
    double abssum3 = 0;
    for (int i = 0; i <= imax; i++) {
        for (int j = 0; j <= jmax; j++) {
            int p = pos[i][j];
            if (p == 0) continue;
            if ((p & 64) != 0) {
                // shortcircuited; need not be updated individually
                if ((p & 1) != 0 && (pos[i-1][j] & 64) == 0) {
                    calculate_Jx_sub(i-1, j);
                    w2 += weight2;
                    double t = fac * Jx[i-1][j];
                    sum2 += t;
                    abssum2 += fabs(t);
                }
                if ((p & 2) != 0 && (pos[i+1][j] & 64) == 0) {
                    calculate_Jx_sub(i, j);
                    w2 -= weight1;
                    double t = fac * Jx[i][j];
                    sum2 -= t;
                    abssum2 += fabs(t);
                }
                if ((p & 4) != 0 && (pos[i][j-1] & 64) == 0) {
                    calculate_Jy_sub(i, j-1);
                    w2 += weight2;
                    double t = fac * Jy[i][j-1];
                    sum2 += t;
                    abssum2 += fabs(t);
                }
                if ((p & 8) != 0 && (pos[i][j+1] & 64) == 0) {
                    calculate_Jy_sub(i, j);
                    w2 -= weight1;
                    double t = fac * Jy[i][j];
                    sum2 -= t;
                    abssum2 += fabs(t);
                }
            } else if ((p & 128) != 0) {
                // shortcircuited; need not be updated individually
                if ((p & 1) != 0 && (pos[i-1][j] & 128) == 0) {
                    calculate_Jx_sub(i-1, j);
                    w3 += weight2;
                    double t = fac * Jx[i-1][j];
                    sum3 += t;
                    abssum3 += fabs(t);
                }
                if ((p & 2) != 0 && (pos[i+1][j] & 128) == 0) {
                    calculate_Jx_sub(i, j);
                    w3 -= weight1;
                    double t = fac * Jx[i][j];
                    sum3 -= t;
                    abssum3 += fabs(t);
                }
                if ((p & 4) != 0 && (pos[i][j-1] & 128) == 0) {
                    calculate_Jy_sub(i, j-1);
                    w3 += weight2;
                    double t = fac * Jy[i][j-1];
                    sum3 += t;
                    abssum3 += fabs(t);
                }
                if ((p & 8) != 0 && (pos[i][j+1] & 128) == 0) {
                    calculate_Jy_sub(i, j);
                    w3 -= weight1;
                    double t = fac * Jy[i][j];
                    sum3 -= t;
                    abssum3 += fabs(t);
                }
            } else {
                double w = 0;
                double sum = 0;
                double abssum = 0;
                if ((p & 1) != 0) {
                    calculate_Jx_sub(i-1, j);
                    w += weight2;
                    double t = fac * Jx[i-1][j];
                    sum += t;
                    abssum += fabs(t);
                }
                if ((p & 2) != 0) {
                    calculate_Jx_sub(i, j);
                    w -= weight1;
                    double t = fac * Jx[i][j];
                    sum -= t;
                    abssum += fabs(t);
                }
                if ((p & 4) != 0) {
                    calculate_Jy_sub(i, j-1);
                    w += weight2;
                    double t = fac * Jy[i][j-1];
                    sum += t;
                    abssum += fabs(t);
                }
                if ((p & 8) != 0) {
                    calculate_Jy_sub(i, j);
                    w -= weight1;
                    double t = fac * Jy[i][j];
                    sum -= t;
                    abssum += fabs(t);
                }
                if ((p & 16) == 0) {
                    V[i][j] -= sum / w;
                    double e = fabs(sum) / (abssum + 1);
                    if (e > emax) emax = e;
                }
            }
        }
    }
    circ[2].v -= sum2 / w2;
    double e = fabs(sum2) / (abssum2 + 1);
    if (e > emax) emax = e;
    circ[3].v -= sum3 / w3;
    e = fabs(sum3) / (abssum3 + 1);
    if (e > emax) emax = e;
    shortcircuit();
    return emax;
}

istream& nextline(istream& in)
{
    char c;
    while (in.get(c)) {
        if (c == '#') {
            while (in.get(c) && c != '\n')
                ;
        } else if (!std::isspace(c)) {
            in.putback(c);
            break;
        }
    }
    return in;
}

void readinitfile(char* initfile)
{
    ifstream in(initfile);
    if (!in) error(__LINE__);
    while (nextline(in)) {
        int i, j, p;
        double t, v, cx, cy, qx, qy;
        if (!(in >> i >> j >> p >> t >> v >> cx >> cy >> qx >> qy)) error(__LINE__);
        if (i < 0 || i > NX || j < 0 || j > NY) {
            cerr << "readinitfile: i=" << i << ", j=" << j << endl;
        } else {
            T[i][j] = t;
            V[i][j] = v;
            Jx[i][j] = 0;
            Jy[i][j] = 0;
        }
    }
}

void writeoutfile(char* outfile)
{
    ofstream out(outfile);
    if (!out) error(__LINE__);
    for (int i = 0; i <= imax; i++) {
        for (int j = 0; j <= jmax; j++) {

            double cx;
            if (pos[i][j] == 0)
                cx = 0;
            else if ((pos[i][j] & (1|2)) == (1|2))
                cx = (Jx[i-1][j] + Jx[i][j]) / 2;
            else if ((pos[i][j] & 16) == 0)
                cx = 0;
            else if ((pos[i][j] & 1) != 0)
                cx = Jx[i-1][j];
            else
                cx = Jx[i][j];

            double cy;
            if (pos[i][j] == 0)
                cy = 0;
            else if ((pos[i][j] & (4|8)) == (4|8))
                cy = (Jy[i][j-1] + Jy[i][j]) / 2;
            else
                cy = 0;

            double tx;
            if (pos[i][j] == 0) {
                tx = 0;
            } else if ((pos[i][j] & (1|2)) == (1|2)) {
                tx = (T[i+1][j] - T[i-1][j]) / (2 * h);
            } else if ((pos[i][j] & (1|2)) == 1) {
                tx = (T[i][j] - T[i-1][j]) / h;
            } else {
                tx = (T[i+1][j] - T[i][j]) / h;
            }

            double ty;
            if (pos[i][j] == 0) {
                ty = 0;
            } else if ((pos[i][j] & (4|8)) == (4|8)) {
                ty = (T[i][j+1] - T[i][j-1]) / (2 * h);
            } else if ((pos[i][j] & (4|8)) == 4) {
                ty = (T[i][j] - T[i][j-1]) / h;
            } else {
                ty = (T[i][j+1] - T[i][j]) / h;
            }

            double t = T[i][j];
            double k = kappa(t);
            double at = seebeck(t) * t;
            double bt = nernst(t) * B * t;
            double g = k * righi(t) * B;
            double qx = - k * tx + at * cx - bt * cy - g * ty;
            double qy = - k * ty + at * cy + bt * cx + g * tx;
            if (pos[i][j] == 0) qx = qy = 0;

            out << i << ' ' << j << ' ' << pos[i][j];
            out << ' ' << T[i][j] << ' ' << V[i][j];
            out << ' ' << cx << ' ' <<  cy;
            out << ' ' << qx << ' ' <<  qy;
            out << '\n';
        }
        out << '\n';
    }
}

double Ix(int i)                // current in x-direction at i
{
    int j0 = 0;
    while (j0 < jmax && pos[i][j0] == 0) j0++;
    int j1 = jmax;
    while (j1 > 0 && pos[i][j1] == 0) j1--;
    double c = BOUNDFAC * (Jx[i][j0] + Jx[i][j1]);
    for (int j = j0 + 1; j < j1; j++)
        c += Jx[i][j];
    return c * h;
}

double Javg(int i)              // average current density in x-direction
{
    int j0 = jmin[i];
    int j1 = jmax[i];
    double c = BOUNDFAC * (Jx[i][j0] + Jx[i][j1]);
    for (int j = j0 + 1; j < j1; j++)
        c += Jx[i][j];
    int denom = j1 - j0;
    if (denom <= 0) denom = 1;
    return c / denom;
}

double heatflow(int i)          // average of -kappa nabla T in x-direction
{
    int j0 = jmin[i];
    int j1 = jmax[i];
    double q = 0;
    for (int j = j0; j <= j1; j++) {
        double dt;
        int p = pos[i][j] & 3;
        if (p == 3) {
            dt = (T[i+1][j] - T[i-1][j]) / (2 * h);
        } else if (p == 1) {
            dt = (T[i][j] - T[i-1][j]) / h;
        } else if (p == 2) {
            dt = (T[i+1][j] - T[i][j]) / h;
        } else {
            dt = 0;
        }
        if (j == j0 || j == j1) {
            q -= BOUNDFAC * kappa(T[i][j]) * dt;
        } else {
            q -= kappa(T[i][j]) * dt;
        }
    }
    int denom = j1 - j0;
    if (denom <= 0) denom = 1;
    return q / denom;
}

double peltierheat(int i)       // average of alpha T J in x-direction
{
    int j0 = jmin[i];
    int j1 = jmax[i];
    double q = 0;
    for (int j = j0; j <= j1; j++) {
        double c;
        int p = pos[i][j] & 3;
        if (p == 3) {
            c = (Jx[i-1][j] + Jx[i][j]) / 2;
        } else if (p == 1) {
            c = Jx[i-1][j];
        } else if (p == 2) {
            c = Jx[i][j];
        } else {
            c = 0;
        }
        double t = T[i][j];
        if (j == j0 || j == j1) {
            q += BOUNDFAC * seebeck(t) * t * c;
        } else {
            q += seebeck(t) * t * c;
        }
    }
    int denom = j1 - j0;
    if (denom <= 0) denom = 1;
    return q / denom;
}

double ettingshausen(int i)     // average of N T B cross J in x-direction
{
    int j0 = jmin[i];
    int j1 = jmax[i];
    if (j1 <= j0) return 0;
    double q = 0;
    for (int j = j0; j < j1; j++) {
        double t = (T[i][j] + T[i][j+1]) / 2;
        q -= nernst(t) * t * B * Jy[i][j];
    }
    return q / (j1 - j0);
}

double righiheat(int i)         // average of kappa M B cross nabla T in x-direction
{
    int j0 = jmin[i];
    int j1 = jmax[i];
    if (j1 <= j0) return 0;
    double q = 0;
    for (int j = j0; j < j1; j++) {
        double t = (T[i][j] + T[i][j+1]) / 2;
        double dt = (T[i][j+1] - T[i][j]) / h;
        q -= kappa(t) * righi(t) * B * dt;
    }
    return q / (j1 - j0);
}

#if 0
double efficiency(double w)
{
    double h1 = heatflow(imax);
    double p1 = peltierheat(imax);
    double e1 = ettingshausen(imax);
    double r1 = righiheat(imax);
    double a1 = h1 + p1 + e1 + r1;
    return fabs(w / a1);
}
#endif

void show_efficiency(double w)
{
    double h0 = heatflow(0);
    double p0 = peltierheat(0);
    double e0 = ettingshausen(0);
    double r0 = righiheat(0);
    double a0 = h0 + p0 + e0 + r0;
    double h1 = heatflow(imax);
    double p1 = peltierheat(imax);
    double e1 = ettingshausen(imax);
    double r1 = righiheat(imax);
    double a1 = h1 + p1 + e1 + r1;
    cout << "L: " << h0 << ' ' << p0 << ' ' << e0 << ' ' << r0 << ' ' << a0 << '\n';
    cout << "H: " << h1 << ' ' << p1 << ' ' << e1 << ' ' << r1 << ' ' << a1 << '\n';
    cout << "eff = " << fabs(w) << " / " << fabs(a1) << " = " << fabs(w / a1) << '\n';
}

void readparfile(char* parfile)
{
    ifstream in(parfile);
    if (!in) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> rho0 >> drho0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> seebeck0 >> dseebeck0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> kappa0 >> dkappa0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> hall0 >> dhall0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> nernst0 >> dnernst0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> righi0 >> drighi0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> B)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> length) || length <= 0) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> width) || width <= 0) error(__LINE__);

    int jT00, jT01, jT10, jT11, jV00, jV01, jV10, jV11;

    if (!nextline(in)) error(__LINE__);
    if (!(in >> jT00 >> jT01 >> T0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> jT10 >> jT11 >> T1)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> jV00 >> jV01 >> V0)) error(__LINE__);

    if (!nextline(in)) error(__LINE__);
    if (!(in >> jV10 >> jV11 >> V1)) error(__LINE__);
    V_not_fixed = V1 == 999;

    imax = -1;
    while (nextline(in)) {
        int i, j0, j1;
        if (!(in >> i >> j0 >> j1)
            || i > NX || j0 < 0 || j1 < j0 || j1 > NY)
            error(__LINE__);
        while (imax < i) {
            imax++;  jmin[imax] = j0;  jmax[imax] = j1;
        }
    }
    if (imax <= 0) error(__LINE__);
    h = length / imax;

    for (int i = 0; i <= imax; i++) {
        for (int j = jmin[i]; j <= jmax[i]; j++) {
            int a = (i-1 >= 0) && (j >= jmin[i-1]) && (j <= jmax[i-1]);
            int b = (i+1 <= imax) && (j >= jmin[i+1]) && (j <= jmax[i+1]);
            int c = (j-1 >= jmin[i]);
            int d = (j+1 <= jmax[i]);
            pos[i][j] = a | (b << 1) | (c << 2) | (d << 3);
        }
    }
    for (int j = jT00; j <= jT01; j++) pos[0][j] |= 32;
    for (int j = jT10; j <= jT11; j++) pos[imax][j] |= 32;
    for (int j = jV00; j <= jV01; j++) pos[0][j] |= 16;
    for (int j = jV10; j <= jV11; j++) pos[imax][j] |= 16;

    cout << "rho=" << rho0 << "+T*" << drho0 << '\n';
    cout << "seebeck=" << seebeck0 << "+T*" << dseebeck0 << '\n';
    cout << "kappa=" << kappa0 << "+T*" << dkappa0 << '\n';
    cout << "hall=" << hall0 << "+T*" << dhall0 << '\n';
    cout << "nernst=" << nernst0 << "+T*" << dnernst0 << '\n';
    cout << "righi=" << righi0 << "+T*" << drighi0 << '\n';

    cout << "B=" << B << "T, length=" << length << "m, width="
         << width << "m, T0=" << T0 << "K, T1=" << T1 << "K, V0="
         << V0 << "V, V1=";
    if (V_not_fixed)
        cout << "not fixed, ";
    else
        cout << V1 << "V, ";
    cout << "imax=" << imax << ", h=" << h << endl;

    /*
    for (int i = 0; i <= imax; i++) {
        for (int j = jmin[i]; j <= jmax[i]; j++)
            cout << ' ' << pos[i][j];
        cout << endl;
    }
    */
}

void init(char* initfile)
{
    if (initfile) {
        readinitfile(initfile);
    } else {
        for (int i = 0; i <= imax; i++) {
            double t = ((double) i / (double) imax) * (T1 - T0) + T0;
            for (int j = jmin[i]; j <= jmax[i]; j++) {
                T[i][j] = t;
                Jx[i][j] = Jy[i][j] = 0;
            }
        }
        if (V_not_fixed) {
            double v = V0;
            for (int j = jmin[0]; j <= jmax[0]; j++) V[0][j] = v;
            double dt = (T1 - T0) / imax;
            for (int i = 1; i <= imax; i++) {
                double t = (T[i-1][jmin[i-1]] + T[i][jmin[i]]) / 2;
                v -= dt * seebeck(t);
                for (int j = jmin[i]; j <= jmax[i]; j++) V[i][j] = v;
            }
        } else {
            for (int i = 0; i <= imax; i++) {
                double v = ((double) i / (double) imax) * (V1 - V0) + V0;
                for (int j = jmin[i]; j <= jmax[i]; j++)
                    V[i][j] = v;
            }
        }
    }
}

void solve(char* initfile)      // main loop
{
    init(initfile);
    int j0, j1;
    for (j0 = jmin[0]; j0 <= jmax[0]; j0++)
        if ((pos[0][j0] & 16) != 0) break;
    if (j0 > jmax[0]) j0 = (jmin[0] + jmax[0]) / 2;
    for (j1 = jmin[imax]; j1 <= jmax[imax]; j1++)
        if ((pos[imax][j1] & 16) != 0) break;
    if (j1 > jmax[imax]) j1 = (jmin[imax] + jmax[imax]) / 2;
    cout << "iter      emax_T    emax_V         V      W/m2 ";
    if (T1 != T0) {
        cout << "alpha_eff ";
        if (B != 0) cout << "    N_eff";
    } else if (V0 != V1) {
        if (B != 0) cout << "    R_eff";
        cout << "   rho_eff";
    }
    cout << endl;
    (void) update_J_V();
    double dvx = 0;
    double dvy = 0;
    double cur = 0;
    for (int iter = 1; iter <= 1000000; iter++) {
        if (quitting) break;
        double emax_T = 0;
        if (iter % 50 == 0) {
            emax_T = update_T();
            if (diverging()) {
                cerr << "Diverging at SOR = " << sor;
                sor = 0.9 * (sor + 1) - 1;
                cerr << ", New SOR = " << sor << endl;
                init(initfile);
                iter = 1;
            }
        }
        double emax_V = update_J_V();
        if (iter % SHOWSTEP == 0) {
            dvx = V[imax][j1] - V[0][j0];
            dvy = V[imax/2][jmax[imax/2]] - V[imax/2][jmin[imax/2]];
            cur = Javg(imax / 2);
            if (V_not_fixed) {
                double u = 0;
                for (int i = 1; i <= imax; i++) {
                    double t = (T[i-1][(jmin[i-1] + jmax[i-1]) / 2]
                                + T[i][(jmin[i] + jmax[i]) / 2]) / 2;
                    u += cur * rho(t) * h;
                    for (int j = jmin[i]; j <= jmax[i]; j++)
                        V[i][j] += u;
                }
            }
            savefmt dummy(cout);
            cout << setw(6) << iter;
            cout.precision(3);
            cout << ' ' << setw(9) << emax_T;
            cout << ' ' << setw(9) << emax_V;
            cout << ' ' << setw(9) << dvx;
            // cout.precision(6);
            // cout.setf(ios_base::scientific);
            cout << ' ' << setw(9) << dvx * cur;
            if (T1 != T0) {
                double alpha_eff = -dvx / (T1 - T0);
                cout << ' ' << setw(9) << alpha_eff;
                if (B != 0) {
                    double nernst_eff = -dvy * length / (B * (T1 - T0) * width);
                    cout << ' ' << setw(9) << nernst_eff;
                }
            } else if (V0 != V1) {
                if (B != 0) {
                    double hall_eff = -dvy / (B * cur * width);
                    cout << ' ' << setw(9) << hall_eff;
                }
                double rho_eff = (V0 - V1) / (cur * length);
                cout << ' ' << setw(9) << rho_eff;
            }
            cout << endl;
            if (emax_V < 1e-6 && emax_T < 1e-6) break;
        }
    }
    show_efficiency(dvx * cur);
    cout << " R21=" << (V[imax][jmin[imax]] - V[0][jmin[0]]) / (cur * width * 0.002079); // tentative
    cout << " R24=" << (V[imax][jmin[imax]] - V[0][jmax[0]]) / (cur * width * 0.002079); // tentative
    cout << " R31=" << (V[imax][jmax[imax]] - V[0][jmin[0]]) / (cur * width * 0.002079); // tentative
    cout << " R34=" << (V[imax][jmax[imax]] - V[0][jmax[0]]) / (cur * width * 0.002079); // tentative
    cout << endl;
}

void show_coefs()
{
    cout << "T   kappa      rho        hall       seebeck    nernst\n";
    for (int i = 273; i <= 373; i += 20) {
        savefmt dummy(cout);
        cout << setw(3) << i << ' '; 
        cout.precision(2);
        cout.setf(ios_base::scientific);
        cout << setw(10) << kappa(i) << ' '
             << setw(10) << rho(i) << ' '
             << setw(10) << hall(i) << ' '
             << setw(10) << seebeck(i) << ' '
             << setw(10) << nernst(i) << endl;
    }
}

int main(int argc, char* argv[])
{
    if (argc < 3 || argc > 4) {
        cerr <<
            "Usage: nernst parfile outfile [initfile]\n"
            "    parfile   ... name of parameter file\n"
            "    outfile   ... name of output file\n"
            "    initfile  ... name of initial-value file (optional)\n";
        exit(1);
    }
    char* parfile = argv[1];
    char* outfile = argv[2];
    char* initfile = (argc == 4) ? argv[3] : 0;

    readparfile(parfile);

    double t = cos(M_PI / imax) + cos(M_PI / (width / h));
    sor = 4 / (2 + sqrt(4 - t * t)) - 1;

    if (signal(SIGINT, myhandler) == SIG_ERR)
        cerr << "Warning: Unable to catch SIGINT" << endl;

    solve(initfile);

    writeoutfile(outfile);
}