///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Groundwater.java                                   //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of studying the 2-dimensional groundwater
// dynamics through the relaxation method.

import java.lang.*;
public class Groundwater {
  final static int nx = 100, ny = 50, ni = 5;
  public static void main(String argv[]) {
    double sigma0 = 1, a = -0.04, phi0 = 200, b = -20;
    double lx = 1000, hx = lx/nx, ly = 500, hy =ly/ny;
    double phi[][] = new double[nx+1][ny+1];
    double sigma[][] = new double[nx+1][ny+1];
    double f[][] = new double[nx+1][ny+1];
    double p = 0.5;

 // Set up boundary values and a trial solution 
    for (int i=0; i<=nx; ++i) {
      double x = i*hx;
      for (int j=0; j<=ny; ++j) {
        double y = j*hy;
        sigma[i][j]  = sigma0+a*y;
        phi[i][j] = phi0+b*Math.cos(Math.PI*x/lx)*y/ly; 
        f[i][j] = 0;
      }
    }
    for (int step=0; step<ni; ++step) {

   // Ensure boundary conditions by 4-point formula
      for (int j=0; j<ny; ++j) {
         phi[0][j] = (4*phi[1][j]-phi[2][j])/3;
         phi[nx][j] = (4*phi[nx-1][j]-phi[nx-2][j])/3;
      }
      relax2d(p, hx, hy, phi, sigma, f);
    }

 // Output the result
    for (int i=0; i<=nx; ++i) {
      double x = i*hx;
      for (int j=0; j<=ny; ++j) {
        double y = j*hy;
        System.out.println(x + " " + y + " "
          + phi[i][j]);
      }
    }
  }

// Method to perform a relaxation step in 2D.

  public static void relax2d(double p, double hx,
    double hy, double u[][], double d[][],
    double s[][]) {
    double h2 = hx*hx, a = h2/(hy*hy),
      b = 1/(4*(1+a)), ab = a*b, q = 1-p;
    for (int i=1; i<nx; ++i) {
      for (int j = 1; j <ny; ++j) {
        double xp = b*(d[i+1][j]/d[i][j]+1);
        double xm = b*(d[i-1][j]/d[i][j]+1);
        double yp = ab*(d[i][j+1]/d[i][j]+1);
        double ym = ab*(d[i][j-1]/d[i][j]+1);
        u[i][j] = q*u[i][j]+p*(xp*u[i+1][j]
                 +xm*u[i-1][j]+yp*u[i][j+1]
                 +ym*u[i][j-1]+h2*s[i][j]);
      }
    }
  }
}