001package models;
002
003import org.apache.commons.math3.analysis.UnivariateFunction;
004import org.apache.commons.math3.complex.Complex;
005
006/**
007 * An integrable model with scattering function composed of the product of several sinh-Gordon scattering functions.
008 * 
009 * These sinh-Gordon type factors can have different coupling constants, and the coupling constant can be complex,
010 * as long as the complex values appear in complex-conjugate pairs. An overall minus sign in the scattering function,
011 * like in the Ising model, can also be included. 
012 * 
013 */
014public class ProductModel extends IntegrableModel
015{
016    private CachedFunction fminIpi;
017    private int coshPower;
018        
019    /**
020     * Constructs a new product model corresponding to a given list of sinh-Gordon coupling constants,
021     * but no overall minus sign. Identical to {@link #ProductModel(Object[], boolean)} with the last parameter set to <code>false</code>.
022     */
023    public ProductModel(final Object[] coupling) {
024        this(coupling, false);
025    }
026    
027    /**
028     * Constructs a new product model corresponding to a given list of sinh-Gordon coupling constants.
029     * 
030     * These constants can be either real or complex; in the latter case, two complex-conjugate factors are added to the scattering function.
031     * An overall plus/minus sign can also be included.
032     * The overall scattering function is
033     * \[ S(\zeta) = \pm \Big( \prod_{j =1}^n \dfrac{\sinh \zeta - i \sin B_j \pi/2 }{\sinh \zeta + i \sin B_j \pi/2 }\Big)
034     *    \Big( \prod_{j =1}^n \dfrac{\sinh \zeta - i \sin C_j \pi/2 }{\sinh \zeta + i \sin C_j \pi/2 }
035     *           \dfrac{\sinh \zeta - i \sin \bar C_j \pi/2 }{\sinh \zeta + i \sin \bar C_j \pi/2 } \Big)
036     *  \] 
037     * where \(B_j\) are the real and \(C_j\) are the complex coupling constants. The real part of \(B_j,C_j\) is expected 
038     * to be in the interval \( [0,2] \).  
039     * 
040     * Technically, these coupling constants are passed as an array that should contain either <code>Double</code> objects or
041     * <code>org.apache.commons.math3.complex.Complex</code> objects, or a combination of both. The plus/minus sign is determined
042     * by the parameter <code>hasMinus</code>
043     * 
044     * @param coupling the list of coupling constants, as described above
045     * @param hasMinus whether an overall minus sign is included in the S function.
046     */
047    public ProductModel(final Object[] coupling, final boolean hasMinus) {
048        
049        int nCoupling = 0;
050        for (Object o: coupling) {
051                if (o instanceof Complex) {
052                        nCoupling += 2;
053                } else {
054                        nCoupling += 1;
055                }
056        }
057        
058        if (hasMinus) {
059                this.coshPower =  1 - 2 * ((nCoupling+1) / 2);
060 
061        } else {
062                this.coshPower = -2 * (nCoupling / 2);
063        }
064                
065        UnivariateFunction fminI = new UnivariateFunction() {
066            public double value(double theta) {
067                double prod = 1.0;
068                for (Object b : coupling) {
069                        if (b instanceof Double) {
070                                prod *= SinhTools.sinhIntegral(theta, (Double) b);
071                        } else if (b instanceof Complex) {
072                                prod *= SinhTools.sinhIntegralCCPair(theta, (Complex) b);
073                        } else {
074                                throw new IllegalArgumentException("coupling neither Double nor Complex: "+b);
075                        }
076                }
077                return prod;
078            }
079        };
080       
081        fminIpi = new CachedFunction(fminI, 1e-14);
082    }
083    
084    @Override
085    public double FminIpi(double theta) {
086        double val = fminIpi.value(Math.abs(theta));
087                val *= Math.pow(Math.cosh(theta/2.0), this.coshPower);
088        return val;
089    }    
090    
091}