package com.srbenoit.modeling.grid;

import com.srbenoit.math.grapher.Graphable;
import com.srbenoit.math.grapher.Grapher;

/**
 * A class that precomputes soft-sphere forces and energies at some number of discrete points then
 * returns an energy or force magnitude based on a square of the ratio of actual distance to
 * equilibrium distance. For any distance ratio greater than 1, the class returns zero energy and
 * force.
 *
 * <p>u = dist / EQ_DIST
 *
 * <p>E = WELL_DEPTH [ (1 / u)^12 - 1 ]
 *
 * <p>EQ_DIST * |F| = 12 (WELL_DEPTH / u) (1 / u)^12 where force is directed away from the particle
 * (positive forces are repulsive).
 */
public class FastSoftSphere implements Graphable {

    /** the number of discrete points to use in the representation of the potential */
    private final static int NUM_POINTS = 500;

    /** the default domain for graphing the function */
    private final double[] domain;

    /** the default range for graphing the function */
    private final double[] range;

    /** the precomputed energies */
    private final double[] energies;

    /** the precomputed forces */
    private final double[] forces;

    /**
     * Constructs a new <code>FastSoftSphere</code>.
     *
     * @param  wellDepth  the potential well depth.
     */
    public FastSoftSphere(final double wellDepth) {

        double uVal;
        double square;
        double sixth;
        double twelfth;
        int last;

        this.energies = new double[NUM_POINTS];
        this.forces = new double[NUM_POINTS];

        // Compute energies
        for (int i = 0; i < NUM_POINTS; i++) {

            // if dist ranges from 0 to EQ_DIST, u ranges from 0 to 1
            uVal = (i + 0.5) / NUM_POINTS;

            square = 1 / uVal;
            sixth = square * square * square;
            twelfth = sixth * sixth;

            this.energies[i] = wellDepth * (square - 1);
            this.forces[i] = 12 * wellDepth * square / uVal;
        }

        // Now adjust so the last term is zero
        last = NUM_POINTS - 1;
        for (int i = 0; i < NUM_POINTS; i++) {
            this.forces[i] -= this.forces[last];
        }

        this.domain = new double[] { 0, 1.5 };
        this.range = new double[] { -wellDepth, 30 * wellDepth };
    }

    /**
     * Computes the energy given a ratio of distance / equilibrium distance.
     *
     * @param   uVal  the square of the ratio of actual distance to equilibrium distance
     * @return  the energy
     */
    public double energy(final double uVal) {

        int index;

        index = (int) (uVal * NUM_POINTS);

        if (index < 0) {
            index = 0;
        } else if (index >= NUM_POINTS) {
            index = NUM_POINTS - 1;
        }

        return this.energies[index];
    }

    /**
     * Computes the force given a ratio of distance / equilibrium distance.
     *
     * @param   uVal  the square of the ratio of actual distance to equilibrium distance
     * @return  the force
     */
    public double forceTimesEqDist(final double uVal) {

        int index;

        index = (int) (uVal * NUM_POINTS);

        if (index < 0) {
            index = 0;
        } else if (index >= NUM_POINTS) {
            index = NUM_POINTS - 1;
        }

        return this.forces[index];
    }

    /**
     * Gets the number of dimensions of the graph (2 for a function of a single value).
     *
     * @return  the number of dimensions of the graph
     */
    public int graphDimensions() {

        return 2;
    }

    /**
     * Gets a default domain over which the graph will show the main features of the function. This
     * domain should be computed based on known attributes of the function.
     *
     * @return  a two-double array containing the left and right endpoints of the default domain
     */
    public double[] defaultDomain() {

        return this.domain.clone();
    }

    /**
     * Gets a default range over which the graph will show the main features of the function. This
     * range should be computed based on known attributes of the function.
     *
     * @return  a two-double array containing the lower and upper limits of the default range
     */
    public double[] defaultRange() {

        return this.range.clone();
    }

    /**
     * Computes the graph values at a coordinate or coordinates. The number of coordinates needed
     * is one less than the dimension.
     *
     * @param   coordinates  the list of coordinates
     * @return  the graph values at that location
     */
    public double valueAt(final double... coordinates) {

// return energy(coordinates[0]);
        return forceTimesEqDist(coordinates[0]);
    }

    /**
     * Main method to graph the force vs. distance.
     *
     * @param  args  command-line arguments
     */
    public static void main(final String... args) {

        FastSoftSphere obj;
        Grapher grapher;

        obj = new FastSoftSphere(1);

        grapher = new Grapher(500, 500);
        grapher.graph(obj);
        grapher.showInFrame("Fast Soft Sphere");
    }
}