///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Scattering.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 quantum scattering in one 
// dimension through the Numerov and steepest-descent 
// methods.

import java.lang.*;
import java.io.*;
public class Scattering {
  static final int nx = 100; 
  static final double hartree = 0.0116124;
  static final double bohr = 98.964, a = 125/bohr;
  static double e; 
  
  public static void main(String argv[]) throws
    IOException {
    double x[] = new double[2];

// Read in the particle energy
    System.out.println("Enter particle energy: ");
    InputStreamReader c
      = new InputStreamReader(System.in);
    BufferedReader b = new BufferedReader(c);
    String energy = b.readLine();
    e = Double.valueOf(energy).doubleValue();
    e /= hartree;
    x[0] =  0.1;
    x[1] = -0.1;
    double d = 0.1, del = 1e-6;
    steepestDescent(x, d, del);
    double r = x[0]*x[0]+x[1]*x[1];
    double t = 1-r;
    System.out.println("The reflectivity: " + r);
    System.out.println("The transmissivity: " + t);
  }

// Method to carry out the steepest descent search.
  public static void steepestDescent(double x[],
    double a, double del) {
    int n = x.length;
    double h = 1e-6;
    double g0 = g(x);
    double fi[] = new double[n];
    fi = f(x, h);
    double deriv = 0;
    for (int i=0; i<n; ++i) deriv += fi[i]*fi[i];
    deriv = Math.sqrt(deriv);
    double b = a/deriv;
    while (deriv > del) {
      for (int i=0; i<n; ++i) x[i] -= b*fi[i];
      h /= 2;
      fi = f(x, h);
      deriv = 0;
      for (int i=0; i<n; ++i) deriv += fi[i]*fi[i];
      deriv = Math.sqrt(deriv);
      b  = a/deriv;
      double g1 = g(x);
      if (g1 > g0) a /= 2;
      else g0 = g1;
    }
  }

// Method to provide function f=gradient g(x).
  public static double[] f(double x[], double h) {
    int n = x.length;
    double z[] = new double[n];
    double y[] = (double[]) x.clone();
    double g0 = g(x);
    for (int i=0; i<n; ++i) {
      y[i] += h;
      z[i] = (g(y)-g0)/h;
    }
    return z;
  }

// Method to provide function g(x).
  public static double g(double x[]) {
    double h = a/nx;
    double ar = x[0];
    double ai = x[1];
    double ur[] = new double[nx+1];
    double ui[] = new double[nx+1];
    double qf[] = new double[nx+1];
    double qb[] = new double[nx+1];
    double s[] = new double[nx+1];

    for (int i=0; i<=nx; ++i) {
      double xi = h*i; 
      s[i] = 0;
      qf[i] = 2*(e-v(xi));
      qb[nx-i] = qf[i];
    }

 // Perform forward integration
    double delta = 1e-6;
    double k = Math.sqrt(2*e);
    double ur0 = 1+ar;
    double ur1 = ur0+k*ai*h
      +h*h*(v(delta)-e)*ur0; 
    double ui0 = ai;
    double ui1 = ui0+k*(1-ar)*h
      +h*h*(v(delta)-e)*ui0; 
    ur = numerov(nx+1, h, ur0, ur1, qf, s);
    ui = numerov(nx+1, h, ui0, ui1, qf, s);

 // Perform backward integration
    ur0 = ur[nx];
    ur1 = ur0+k*ui[nx]*h
      +h*h*(v(a-delta)-e)*ur0; 
    ui0 = ui[nx];
    ui1 = ui0-k*ur[nx]*h
      +h*h*(v(a-delta)-e)*ui0; 
    ur = numerov(nx+1, h, ur0, ur1, qb, s);
    ui = numerov(nx+1, h, ui0, ui1, qb, s);

    return (1+ar-ur[nx])*(1+ar-ur[nx])
      +(ai-ui[nx])*(ai-ui[nx]);
  }

// Method to provide the potential V(x).
  public static double v(double x) {
    double v0 = 0.3/hartree;
    double x1 = 25/bohr;
    double x2 = 75/bohr;
    if (((x<x1)&&(x>0)) || ((x<a)&&(x>x2)))
      return v0;
    else return 0;
  }

// Method to perform the Numerov integration.
  public static double[] numerov(int m, double h,
    double u0, double u1, double q[], double s[]) {
    double u[] = new double[m];
    u[0] = u0;
    u[1] = u1;
    double g = h*h/12;
    for (int i=1; i<m-1; ++i) {
      double c0 = 1+g*q[i-1];
      double c1 = 2-10*g*q[i];
      double c2 = 1+g*q[i+1];
      double d  = g*(s[i+1]+s[i-1]+10*s[i]);
      u[i+1] = (c1*u[i]-c0*u[i-1]+d)/c2;
    }
    return u;
  }
}