/* This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
package etomica.cavity;

import etomica.atom.IAtom;
import etomica.atom.IAtomKinetic;
import etomica.atom.IAtomList;
import etomica.potential.P2HardSphere;
import etomica.space.Space;

public class P2HardSphereCavity extends P2HardSphere {

    protected double pairWell;
    protected IAtomKinetic pairedAtom1, pairedAtom2;
    protected int idxSum, idxProduct;

    public enum CollisionType {INTERNAL_BOUNCE, ESCAPE, EXTERNAL_BOUNCE, CAPTURE}

    protected CollisionType lastCollisionType;

    public P2HardSphereCavity(Space space) {
        this(space, 0);
    }

    public P2HardSphereCavity(Space space, double pairWell) {
        super(space, 1, false);
        this.pairWell = pairWell;
    }

    public IAtom getPairedAtom1() {
        return pairedAtom1;
    }

    public IAtom getPairedAtom2() {
        return pairedAtom2;
    }

    public double collisionTime(IAtomList pair, double falseTime) {
        int idx0 = pair.get(0).getLeafIndex();
        int idx1 = pair.get(1).getLeafIndex();
        boolean paired = idx0 * idx1 == idxProduct && idx0 + idx1 == idxSum;
        if (!paired) {
            // this pair is not overlapping, treat as normal
            return super.collisionTime(pair, falseTime);
        }
        // this pair is overlapping.  look for escape
        IAtomKinetic atom0 = (IAtomKinetic) pair.get(0);
        IAtomKinetic atom1 = (IAtomKinetic) pair.get(1);
        dv.Ev1Mv2(atom1.getVelocity(), atom0.getVelocity());

        dr.Ev1Mv2(atom1.getPosition(), atom0.getPosition());
        dr.PEa1Tv1(falseTime, dv);
        boundary.nearestImage(dr);

        double bij = dr.dot(dv);
        double v2 = dv.squared();
        double discriminant = bij * bij - v2 * (dr.squared() - sig2);
        if (discriminant < 0) {
            throw new RuntimeException("paired atoms not overlapping?");
        }

        double time = (-bij + Math.sqrt(discriminant)) / v2;
        return time + falseTime;
    }

    /**
     * Implements collision dynamics and updates lastCollisionVirial
     */
    public void bump(IAtomList pair, double falseTime) {
        int idx0 = pair.get(0).getLeafIndex();
        int idx1 = pair.get(1).getLeafIndex();
        boolean paired = idx0 * idx1 == idxProduct && idx0 + idx1 == idxSum;
        if (idxSum > 0 && (!paired || pairWell == Double.POSITIVE_INFINITY)) {
            lastCollisionType = paired ? CollisionType.INTERNAL_BOUNCE : CollisionType.EXTERNAL_BOUNCE;
            super.bump(pair, falseTime);
            return;
        }

        IAtomKinetic atom0 = (IAtomKinetic) pair.get(0);
        IAtomKinetic atom1 = (IAtomKinetic) pair.get(1);
        dv.Ev1Mv2(atom1.getVelocity(), atom0.getVelocity());

        dr.Ev1Mv2(atom1.getPosition(), atom0.getPosition());
        dr.PEa1Tv1(falseTime, dv);
        boundary.nearestImage(dr);

        double r2 = dr.squared();
        double bij = dr.dot(dv);

        if (paired != (bij > 0.0)) {         // Separating
            throw new RuntimeException("paired: " + paired + " and bij " + bij);
        }

        double rm0 = atom0.getType().rm();
        double rm1 = atom1.getType().rm();
        double reducedMass = 2.0 / (rm0 + rm1);
        double ke = bij * bij * reducedMass / (2.0 * r2);

        if (bij > 0) {
            if (ke < pairWell) {     // Not enough kinetic energy to escape
                lastCollisionVirial = reducedMass * bij;
                lastCollisionType = CollisionType.INTERNAL_BOUNCE;
            } else {                 // Escape
                lastCollisionVirial = 0;
                pairedAtom1 = pairedAtom2 = null;
                idxSum = idxProduct = 0;
                lastCollisionType = CollisionType.ESCAPE;
            }
        } else if (ke > -pairWell) {   // Approach/capture
            // core collision -- capture
            lastCollisionVirial = 0;
            pairedAtom1 = atom0;
            pairedAtom2 = atom1;
            idxSum = idx0 + idx1;
            idxProduct = idx0 * idx1;
            lastCollisionType = CollisionType.CAPTURE;
        } else {
            // not enough kinetic energy to overcome shoulder
            lastCollisionVirial = reducedMass * bij;
            lastCollisionType = CollisionType.EXTERNAL_BOUNCE;
        }

        lastCollisionVirialr2 = lastCollisionVirial / r2;
        //dv is now change in velocity due to collision
        dv.Ea1Tv1(lastCollisionVirialr2, dr);
        atom0.getVelocity().PEa1Tv1(rm0, dv);
        atom1.getVelocity().PEa1Tv1(-rm1, dv);
        atom0.getPosition().PEa1Tv1(-falseTime * rm0, dv);
        atom1.getPosition().PEa1Tv1(falseTime * rm1, dv);
    }

    public CollisionType getLastCollisionType() {
        return lastCollisionType;
    }
}

@@ -381,16 +381,52 @@
 e.Space;
+%0Aimport etomica.util.random.IRandom;
 %0A%0Apublic
@@ -475,16 +475,52 @@
 here %7B%0A%0A
+    protected final IRandom random;%0A
     prot
@@ -829,16 +829,22 @@
 s(space,
+ null,
  0);%0A   
@@ -889,16 +889,32 @@
 e space,
+ IRandom random,
  double 
@@ -953,24 +953,54 @@
  1, false);%0A
+        this.random = random;%0A
         this
@@ -3509,111 +3509,96 @@
 1);%0A
-        double ke = bij * bij * reducedMass / (2.0 * r2);%0A%0A        if (bij %3E 0) %7B%0A            if (ke %3C 
+%0A        if (bij %3E 0) %7B%0A            if (pairWell %3E 0 && random.nextDouble() %3C Math.exp(-
 pair
@@ -3598,24 +3598,25 @@
 p(-pairWell)
+)
  %7B     // No
@@ -4047,14 +4047,55 @@
 if (
-ke %3E -
+pairWell %3E 0 %7C%7C random.nextDouble() %3C Math.exp(
 pair
@@ -4099,16 +4099,17 @@
 airWell)
+)
  %7B   // 

determine escape/capture without KE (use MC-type probabilities)