package simulate;

/**
 * Meter for tabulation of the atomic radial distribution function (RDF)
 * Should work ok in ensembles with fluctuating volume and particle numbers
 * Not suited for tabulation of RDFs of mixtures or heteroatomic molecules
 */
public class MeterRDF extends MeterFunction {
    
    private AtomPair.Iterator iterator;
    private double[] vShell;
    private double delr;
    
    public MeterRDF() {
	    super();
	    setXLabel("r");
	    setLabel("rdf");
    }
    
    /**
     * Sets phase, evaluates abscissa values to maximum of half system edge-length, and constructs atom iterator
     */
	public void setPhase(Phase p) {
	    super.setPhase(p);
	    iterator = phase.iteratorFactory().makeAtomPairIteratorAll();
	    setX(0.0, 0.5*phase.boundary().dimensions().component(0),nPoints);
	}
	
	/**
	 * Computes RDF for the current configuration
	 *    For future development: It may be possible to extend to particular atom pairs by changing iterator and using a different normalization
	 */
	public double[] currentValue() {
	    iterator.reset();                           //prepare iterator of atom pairs
	    for(int i=0; i<nPoints; i++) {y[i] = 0.0;}  //zero histogram
	    while(iterator.hasNext()) {                 //iterate over all pairs
	        double r = Math.sqrt(iterator.next().r2()); //compute pair separation
	        if(r < xMax) {
	            int index = (int)(r/delr);          //determine histogram index
	            y[index]+=2;                        //add once for each atom
	        }
	    }
	    int n = phase.atomCount();                  //compute normalization: divide by
	    double norm = n*n/phase.volume();           //n, and density*(volume of shell)
	    for(int i=0; i<nPoints; i++) {y[i] /= (norm*vShell[i]);}
	    return y;
	}
	
	/**
	 * Sets the radial values at which RDF is tabulated
	 */
    public void setX(double min, double max, int n) {
	    super.setX(min, max, n);
	    if(phase == null) return;
	    //Compute normalization constants for RDF, including shell volumes for ideal-gas particle numbers
	    vShell = new double[n];         
	    Space space = Simulation.space;
	    double dx2 = 0.5*(xMax - xMin)/(double)nPoints;
	    for(int i=0; i<nPoints-1; i++) {
	        vShell[i] = space.sphereVolume(x[i]+dx2)-space.sphereVolume(x[i]-dx2);
	    }
	    delr = xMax/(double)(nPoints-1);
	}
	    
}