import java.io.*;
import java.lang.*;

public class H2_1D_analyticSolution {
  public static void main( String[] args ) {

    int i, index;
    double beta_a, kappa, x, B; 
    double MaxX;
    double aa = 0.0, bb = 0.0, pp = 0.0, Faa = 0.0, Fbb = 0.0, Fpp = 0.0;
    double normalize;

//  Hardwired data for now
    int  MaxFixPoint   = 100;
    int  MaxBissection = 100;
    int Nmesh = 500; // must be even
    double Z1 = 2.0;
    double a = 0.3;
    double TOL = 0.000000000001;

    double [] wave, wave1, wave2, beta; 
    double [] density, density1, density2, density1iso, density2iso, densitySUMiso, densitySUM;
    double [] waveSYME, waveANTI;
    double [] wave1iso, wave2iso;

    MaxX = a*5; // adds 500 percent at the right of the atomic position

    wave  = new double [2*Nmesh+1];
    wave1 = new double [2*Nmesh+1];
    wave2 = new double [2*Nmesh+1];
    waveSYME = new double [2*Nmesh+1];
    waveANTI = new double [2*Nmesh+1];

    wave1iso = new double [2*Nmesh+1];
    wave2iso = new double [2*Nmesh+1];

    density  = new double [2*Nmesh+1];
    density1 = new double [2*Nmesh+1];
    density2 = new double [2*Nmesh+1];
    densitySUM = new double [2*Nmesh+1];

    density1iso = new double [2*Nmesh+1];
    density2iso = new double [2*Nmesh+1];
    densitySUMiso = new double [2*Nmesh+1];

    beta = new double [2*Nmesh+1];

    FileOutputStream TotalOut, FileData_Total; 
    FileOutputStream PartitionOut, FileData_Partition;
    FileOutputStream IsolatedOut, FileData_Isolated;
    PrintStream Total, Partition, Isolated; // declare a print stream object

//  simple clear screen
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");

    System.out.println("This program generates the example by Morrel Cohen");
    System.out.println("''Partition Theory: A Very Simple Illustration''");
    System.out.println("J. Phys. Chem. A 2007, 111, 12447-12453");
    System.out.println("Program written by carrier@cs.umn.edu\n");


// PART 1: 2 hydrogenoids system  Section 2. of article
//
    // fix point algorithm
    kappa = 0.9; // initial guess
    for (i=0; i < MaxFixPoint; i++) {
      kappa = 2.0*Z1/(1.0 + Math.tanh(kappa*a));
    }
    System.out.printf( "kappa=%2.12f", kappa );
    System.out.printf( "    Residual=%2.12f\n", 
                      kappa -2.0*Z1/(1.0 + Math.tanh(kappa*a)) );

    System.out.printf( "Eigenvalue [-kappa^2/2] =%2.12f\n", -kappa*kappa/2 );

    B = kappa / ( 1.0 + kappa*a/Math.cosh(kappa*a)/Math.cosh(kappa*a) 
                      + Math.tanh(kappa*a)); 
    B = Math.sqrt(B);
    System.out.printf( "B=%2.12f\n", B );

    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      x = i*MaxX/Nmesh;
      if (Math.abs(x) >= a ) {
         wave[index] = B * Math.exp(a - Math.abs(x));
      } else {
         wave[index] = B * Math.cosh(kappa*x)/Math.cosh(kappa*a);
      }
    }

    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      density[index] = 2.0*wave[index]*wave[index];
    }


// PART 2: partitioned system  Section 4. and 5. of article
//
    // fix point iteration does not converge with that function
    // use bissection method instead
    // First, find end points by hand
    aa = 0.0; // assume positive root
    Faa = betaFunc(aa, kappa, a, Z1);
    for (i=0; i<100; i++) {
      bb = aa + 1.0; // initial guess
      Fbb = betaFunc(bb, kappa, a, Z1);
      if (Faa*Fbb < 0) break;
    }

    Faa = betaFunc(aa, kappa, a, Z1);
    Fbb = betaFunc(bb, kappa, a, Z1);
    if (Faa*Fbb > 0) {
      System.out.printf("Could not find bounds for the bissection; STOP\n");
      System.out.printf( "F(%2.12f)=%2.12f\n", aa, Faa);
      System.out.printf( "F(%2.12f)=%2.12f\n", bb, Fbb);
      System.exit(0);
    }

    for (i=0; i < MaxBissection; i++) {

      pp = aa + (bb - aa)/2.0;
      Fpp = betaFunc(pp, kappa, a, Z1); 
      if (Math.abs(Fpp) < TOL || Math.abs((bb - aa)/2.0) < TOL) {
        System.out.printf( "F(%2.12f)=%2.12f found after %4d iterations\n", 
                                pp,     Fpp,              i               );
        break; // found it; stop the loop
      }
      Faa = betaFunc(aa, kappa, a, Z1);
      if (Faa*Fpp > 0) {
       aa = pp; 
      } else {
       bb = pp; 
      };
    }
    beta_a = pp;

//  Now calculate beta(x) on the same mesh as the 2 hydrogenoids
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      x = i*MaxX/Nmesh;
      if (x > a ) {
         beta[index] =  beta_a; // that is equation (5.6)
      }
      if (x < -a ) {
         beta[index] = -beta_a;
      }
      if (0 <= Math.abs(x) && Math.abs(x) <= a ) {
        beta[index] = Math.tanh(kappa*x) /
                      Math.tanh(kappa*a) * beta_a; 
      }

    }
//  now determine the corresponding partition densities
//  all the index do correspond, therefore no need to use x anymore
//  x and i are 1 to 1.
// first normalize the wave functions to the density calculated in PART 1
// That is the total density, that is a constraint here CAREFUL
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      waveSYME[index] = Math.sqrt(density[index])*Math.cos(beta[index]);
      waveANTI[index] = Math.sqrt(density[index])*Math.sin(beta[index]);
    }

// use definition of wave1 and wave2 eq. 3.4
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      wave1[index] = 1.0/Math.sqrt(2.0)*(waveSYME[index] + waveANTI[index]);
      wave2[index] = 1.0/Math.sqrt(2.0)*(waveSYME[index] - waveANTI[index]);
    }
   
// use definition 3.2 to get the partitioned densities
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      density1[index] = wave1[index]*wave1[index];
      density2[index] = wave2[index]*wave2[index];
      densitySUM[index] = density1[index] + density2[index];
    }

// normalize to the partition densities the isolated density
    normalize = 1.0;
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      if (density1[index] > normalize) normalize = density1[index];
    }
    normalize /= Z1;

// PART 3: solution of isolated atom centered at +a and -a
    for (i=-Nmesh; i < Nmesh; i++) {
      index = i + Nmesh;
      x = i*MaxX/Nmesh;
      wave1iso[index] = Math.sqrt(Z1)*Math.exp(-Z1*Math.abs(x-a));
      density1iso[index] = wave1iso[index]*wave1iso[index];

      wave2iso[index] = Math.sqrt(Z1)*Math.exp(-Z1*Math.abs(x+a));
      density2iso[index] = wave2iso[index]*wave2iso[index];
      
      densitySUMiso[index] = density1iso[index] + density2iso[index];
    }

   
// prints to file out
    try
    {
      TotalOut = new FileOutputStream("TotalDensity.out"); 
      Total = new PrintStream( TotalOut ); 

      PartitionOut = new FileOutputStream("PartitionDensity.out"); 
      Partition = new PrintStream( PartitionOut ); 

      IsolatedOut = new FileOutputStream("IsolatedDensity.out"); 
      Isolated = new PrintStream( IsolatedOut ); 

      FileData_Total = new FileOutputStream("density.BIN");
      DataOutputStream Data_out = new DataOutputStream( FileData_Total );

      for (i=-Nmesh; i < Nmesh; i++) {
        index = i + Nmesh;
        x = i*MaxX/Nmesh;
        Total.printf ("%3.6f %3.6f\n", x, density[index]);
        Partition.printf ("%3.6f %3.6f %3.6f %3.6f\n", 
        x, density1[index], density2[index], densitySUM[index]);
        Isolated.printf ("%3.6f %3.6f %3.6f %3.6f      %3.6f   %3.6f   %3.6f\n", 
        x,           density1iso[index],           density2iso[index],          densitySUMiso[index],
           normalize*density1iso[index], normalize*density2iso[index],normalize*densitySUMiso[index]);
      }

      Data_out.writeInt(1);
  
      Total.close();
      Data_out.close();
    }
    catch (Exception e)
    {
      System.err.println ("Error writing to file");
    }


  }

  public static double betaFunc(double beta_a, double kappa, double a, double Z) {
    double beta;

    beta = beta_a - Z/(2.0*kappa)*(Math.sinh(2.0*kappa*a)*Math.cos(2.0*beta_a));

    return beta;
  }
}