///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Bench.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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// A program to solve the problem of a person sitting
// on a bench as described in the text.

import java.lang.*;
public class Bench {
  final static int n = 99, m = 2;
  public static void main(String argv[]) {
  double d[] = new double[n];
  double b[] = new double[n];
  double c[] = new double[n];
  double l = 3, l2 = l/2, h = l/(n+1), h2 = h*h;
  double x0 = 0.25, x2 = x0*x0, e0 = 1/Math.E;
  double rho = 3, g = 9.8, f0 = 200;
  double y = 1e9*Math.pow(0.03,3)*0.2/3;
  
 // Evaluate the coefficient matrix elements
    for (int i=0; i<n; ++i) {
      d[i] = -2;
      c[i] =  1;
      b[i] = -rho*g;
      double x = h*(i+1)-l2;
      if (Math.abs(x) < x0)
        b[i] -= f0*(Math.exp(-x*x/x2)-e0);
      b[i] *= h2/y;
    }

 // Obtain the solution of the curverture of the bench
    double u[] = tridiagonalLinearEq(d, c, c, b);

 // Output the result in every m time steps
    double x = h;
    double mh = m*h;
    for (int i=0; i<n; i+=m) {
      System.out.println(x + " " + 100*u[i]);
      x += mh;
    }
  }

// Method to solve the tridiagonal linear equation set.

  public static double[] tridiagonalLinearEq(double d[],
    double e[], double c[], double b[]) {
    int m = b.length;
    double w[] = new double[m];
    double y[] = new double[m];
    double z[] = new double[m];
    double v[] = new double[m-1];
    double u[] = new double[m-1];

 // Evaluate the elements in the LU decomposition
    w[0] = d[0];
    v[0]  = c[0];
    u[0]  = e[0]/w[0];
    for (int i=1; i<m-1; ++i) {
      w[i]  = d[i]-v[i-1]*u[i-1];
      v[i]  = c[i];
      u[i]  = e[i]/w[i];
    }
    w[m-1]  = d[m-1]-v[m-2]*u[m-2];

 // Forward substitution to obtain y
    y[0] = b[0]/w[0];
    for (int i=1; i<m; ++i)
      y[i] = (b[i]-v[i-1]*y[i-1])/w[i];

 // Backward substitution to obtain z
    z[m-1] = y[m-1];
    for (int i=m-2; i>=0; --i) {
      z[i] = y[i]-u[i]*z[i+1];
    }
    return z;
  }
}