/* 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.virial;

import etomica.box.Box;
import etomica.math.SpecialFunctions;
import etomica.molecule.IMoleculeList;
import etomica.potential.IPotential;
import etomica.potential.Potential2Spherical;

/**
 * Mayer function class that returns temperature expansion bond of the desired
 * order.
 * 
 * @author Andrew Schultz
 */
public class MayerSphericalPTAtt implements MayerFunction {

	/**
	 * Constructor takes reference potential, attractive component and order in temperature.
	 */
	public MayerSphericalPTAtt(Potential2Spherical potentialRef, Potential2Spherical potentialAtt, int order) {
		this.potentialRef = potentialRef;
		this.potentialAtt = potentialAtt;
		this.order = order;
		this.factOrder = SpecialFunctions.factorial(order);
	}

	public double f(IMoleculeList pair, double r2, double beta) {
		double uRef = potentialRef.u(r2);
		if (Double.isInfinite(uRef)) {
			return 0;
		}
		// double dfdkT = Math.exp(-beta*uRef)*Math.pow(-potentialAtt.u(r2)*beta, order)/SpecialFunctions.factorial(order);
		double dfdkT = Math.exp(-beta*uRef)/factOrder;
		if (order > 0) {
		    double betaUA = -potentialAtt.u(r2)*beta;
		    switch (order) {
		        case 1: 
		            dfdkT *= betaUA;
		            break;
		        case 2: 
		            dfdkT *= betaUA*betaUA;
		            break;
                case 3: 
                    dfdkT *= betaUA*betaUA*betaUA;
                    break;
                case 4:
                    double betaUA2 = betaUA*betaUA;
                    dfdkT *= betaUA2*betaUA2;
                    break;
                case 5:
                    betaUA2 = betaUA*betaUA;
                    dfdkT *= betaUA2*betaUA2*betaUA;
                    break;
                case 6:
                    double betaUA3 = betaUA*betaUA*betaUA;
                    dfdkT *= betaUA3*betaUA3;
                    break;
                case 7:
                    betaUA3 = betaUA*betaUA*betaUA;
                    dfdkT *= betaUA3*betaUA3*betaUA;
                    break;
                case 8:
                    betaUA2 = betaUA*betaUA;
                    double betaUA4 = betaUA2*betaUA2;
                    dfdkT *= betaUA4*betaUA4;
                    break;
                case 9:
                    betaUA3 = betaUA*betaUA*betaUA;
                    dfdkT *= betaUA3*betaUA3*betaUA3;
                    break;
                default:
                    throw new RuntimeException("oops, please implement order "+order);
		    }
		}
		
		
		if (Double.isNaN(dfdkT)) {
			throw new RuntimeException ("dfdT is NaN");
		}
		return dfdkT;
	}
	
	public void setBox(Box newBox) {
	    potentialRef.setBox(newBox);
	    potentialAtt.setBox(newBox);
	}

    public IPotential getPotential() {
        return potentialRef;
    }

    protected final Potential2Spherical potentialRef, potentialAtt;
	protected final int order;
	protected final long factOrder;
}

@@ -788,19 +788,20 @@
 otential
-Att
+Full
 , int or
@@ -859,19 +859,20 @@
 otential
-Att
+Full
  = poten
@@ -875,19 +875,20 @@
 otential
-Att
+Full
 ;%0A%09%09this
@@ -1298,20 +1298,17 @@
  0) %7B%0A%09%09
-    
+%09
 double b
@@ -1316,16 +1316,17 @@
 taUA = -
+(
 potentia
@@ -1330,26 +1330,37 @@
 tial
-Att
+Full
 .u(r2)
-*
+ - uRef) * 
 beta;%0A%09%09
     
@@ -1355,20 +1355,17 @@
 beta;%0A%09%09
+%09
-    
 switch (
@@ -1379,39 +1379,23 @@
 %7B%0A%09%09
-        
+%09%09
 case 1:
- 
 %0A%09%09
-            
+%09%09%09
 dfdk
@@ -1405,36 +1405,27 @@
 = betaUA;%0A%09%09
-            
+%09%09%09
 break;%0A%09%09   
@@ -1425,39 +1425,23 @@
 ;%0A%09%09
-        
+%09%09
 case 2:
- 
 %0A%09%09
-            
+%09%09%09
 dfdk
@@ -1443,33 +1443,35 @@
 %09dfdkT *= betaUA
-*
+ * 
 betaUA;%0A%09%09      
@@ -1468,72 +1468,35 @@
 ;%0A%09%09
-            break;%0A                case 3: %0A                    
+%09%09%09break;%0A%09%09%09%09case 3:%0A%09%09%09%09%09
 dfdk
@@ -1506,32 +1506,36 @@
 = betaUA
-*
+ * 
 betaUA
-*
+ * 
 betaUA;%0A
         
@@ -1526,36 +1526,21 @@
 betaUA;%0A
-                    
+%09%09%09%09%09
 break;%0A 
@@ -2730,17 +2730,16 @@
 dkT;%0A%09%7D%0A
-%09
 %0A%09public
@@ -2766,20 +2766,17 @@
 Box) %7B%0A%09
-    
+%09
 potentia
@@ -2797,20 +2797,17 @@
 wBox);%0A%09
-    
+%09
 potentia
@@ -2807,19 +2807,20 @@
 otential
-Att
+Full
 .setBox(
@@ -2832,20 +2832,17 @@
 x);%0A%09%7D%0A%0A
-    
+%09
 public I
@@ -2868,24 +2868,18 @@
 ial() %7B%0A
-        
+%09%09
 return p
@@ -2895,19 +2895,13 @@
 ef;%0A
-    %7D%0A%0A    
+%09%7D%0A%0A%09
 prot
@@ -2955,19 +2955,20 @@
 otential
-Att
+Full
 ;%0A%09prote

handle decompositions other than HS/SQW