/*
 * PhaseSrkEos.java
 *
 * Created on 3. juni 2000, 14:38
 */
package thermo.phase;

import thermo.component.*;
import thermo.mixingRule.*;
import MathLib.JamaLinAlg.*;

/**
 *
 * @author  Even Solbraa
 * @version
 */
public class PhaseSrkCPA extends PhaseSrkEos implements PhaseCPAInterface {

    private static final long serialVersionUID = 1000;

    public CPAMixing cpaSelect = new CPAMixing();
    public CPAMixingInterface cpamix;
    boolean updateXold = true;
    double dFdVdXdXdVtotal = 0.0;
    double dFCPAdXdXdTtotal = 0.0, dFCPAdTdT = 0.0;
    double hcpatot = 1.0, hcpatotdT = 0.0, hcpatotdTdT = 0.0;
    private double gcpav = 0.0;
    double lngcpa = 0.0, lngcpav = 0.0, gcpavv = 0.0, gcpavvv = 0.0, gcpa = 0.0;
    int cpaon = 1;
    int[][][] selfAccociationScheme = null;
    int[][][][] crossAccociationScheme = null;
    int totalNumberOfAccociationSites = 0;
    double[][] ksi = null, mk, u = null;
    double[][] Hij, KlkdV, KlkdT, Klk = null;
    int[] moleculeNumber = new int[totalNumberOfAccociationSites];
    int[] assSiteNumber = new int[totalNumberOfAccociationSites];

    /** Creates new PhaseSrkEos */
    public PhaseSrkCPA() {
        super();

    }

    public Object clone() {
        PhaseSrkCPA clonedPhase = null;
        try {
            clonedPhase = (PhaseSrkCPA) super.clone();
        } catch (Exception e) {
            logger.error("Cloning failed.", e);
        }
        // clonedPhase.cpaSelect = (CPAMixing) cpaSelect.clone();
        //clonedPhase.cpamix = (CPAMixingInterface) cpamix.clone();

        return clonedPhase;
    }

    public void setMixingRule(int type) {
        super.setMixingRule(type);
        cpamix = cpaSelect.getMixingRule(1, this);
    }

    public void init(double totalNumberOfMoles, int numberOfComponents, int type, int phase, double beta) { // type = 0 start init type =1 gi nye betingelser
        if (type == 0) {
            selfAccociationScheme = new int[numberOfComponents][0][0];
            crossAccociationScheme = new int[numberOfComponents][numberOfComponents][0][0];
            for (int i = 0; i < numberOfComponents; i++) {
                totalNumberOfAccociationSites += getComponent(i).getNumberOfAssociationSites();
                selfAccociationScheme[i] = cpaSelect.setAssociationScheme(i, this);
                for (int j = 0; j < numberOfComponents; j++) {
                    crossAccociationScheme[i][j] = cpaSelect.setCrossAssociationScheme(i, j, this);
                }
            }
            mk = new double[totalNumberOfAccociationSites][1];
            Klk = new double[totalNumberOfAccociationSites][totalNumberOfAccociationSites];
            KlkdV = new double[totalNumberOfAccociationSites][totalNumberOfAccociationSites];
            KlkdT = new double[totalNumberOfAccociationSites][totalNumberOfAccociationSites];
            Hij = new double[totalNumberOfAccociationSites][totalNumberOfAccociationSites];
            ksi = new double[totalNumberOfAccociationSites][1];
            u = new double[totalNumberOfAccociationSites][1];
            moleculeNumber = new int[totalNumberOfAccociationSites];
            assSiteNumber = new int[totalNumberOfAccociationSites];
        }
        do {
            super.init(totalNumberOfMoles, numberOfComponents, type, phase, beta);
        } while (!solveX());

        if (type != 0 && isConstantPhaseVolume()) {
            try {
                setMolarVolume(molarVolume(pressure, temperature, getA() / numberOfMolesInPhase / numberOfMolesInPhase, getB() / numberOfMolesInPhase, phase));
            } catch (Exception e) {
                e.printStackTrace();
            }
        }

        //     initCPAMatrix();

        if (type > 1) {
            hcpatotdT = calc_hCPAdT();
            //hcpatotdTdT = calc_hCPAdTdT();
            System.out.println("dFCPAdT" + dFCPAdT());
            //calcXsitedT();
            System.out.println("DFdVdT " + dFCPAdTdV() + " dFDtDT" + dFCPAdTdT());
            //         hcpatotdT = calc_hCPAdT();


            //          calcXsitedV();
            //  calcXsitedT();
        }
//        calcXsitedT();
///        initCPAMatrix();
    }

    public void initCPAMatrix() {
        int temp = 0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                moleculeNumber[temp + j] = i;
                assSiteNumber[temp + j] = j;
                ksi[temp + j][0] = ((ComponentSrkCPA) getComponent(i)).getXsite()[j];
                u[temp + j][0] = Math.log(ksi[temp + j][0]) - ksi[temp + j][0] + 1.0;
                mk[temp + j][0] = getComponent(i).getNumberOfMolesInPhase();
            }
            temp = +getComponent(i).getNumberOfAssociationSites();
        }

        for (int i = 0; i < totalNumberOfAccociationSites; i++) {
            for (int j = 0; j < totalNumberOfAccociationSites; j++) {
                Klk[i][j] = mk[i][0] * mk[j][0] / getTotalVolume() * cpamix.calcDelta(assSiteNumber[i], assSiteNumber[j], moleculeNumber[i], moleculeNumber[j], this, getTemperature(), getPressure(), numberOfComponents);
                KlkdV[i][j] = Klk[i][j] * (getGcpav() - 1.0 / getTotalVolume());
                KlkdT[i][j] = Klk[i][j] * cpamix.calcDeltadT(assSiteNumber[i], assSiteNumber[j], moleculeNumber[i], moleculeNumber[j], this, getTemperature(), getPressure(), numberOfComponents) / getGcpa();
            }
        }

        Matrix KlkMatrix = new Matrix(KlkdT);
        Matrix ksiMatrix = new Matrix(ksi);
        //RealMatrixImpl KlkMatrix = new RealMatrixImpl(Klk);
        //RealMatrixImpl ksiMatrix = new RealMatrixImpl(ksi);

        Matrix tempMatrix = ksiMatrix.transpose().times(KlkMatrix.times(ksiMatrix)).times(-0.5);
        // System.out.println("matrix QdV");
        tempMatrix.print(10, 10);

    }

    public void addcomponent(String componentName, double moles, double molesInPhase, int compNumber) {
        super.addcomponent(componentName, moles, molesInPhase, compNumber);
        componentArray[compNumber] = new ComponentSrkCPA(componentName, moles, molesInPhase, compNumber);
    }

    public double getF() {
        return super.getF() + cpaon * FCPA();
    }

    public double dFdT() {
        //System.out.println("Fdt "+ dFCPAdT());
        return super.dFdT() + cpaon * dFCPAdT();
    }

    public double dFdTdV() {
        return super.dFdTdV() + cpaon * dFCPAdTdV();
    }

    public double dFdV() {
        double dv2 = dFCPAdV();
        return super.dFdV() + cpaon * dv2;
    }

    @Override
    public double dFdVdV() {
        return super.dFdVdV() + cpaon * dFCPAdVdV();
    }

    public double dFdVdVdV() {
        return super.dFdVdVdV() + cpaon * dFCPAdVdVdV();
    }

    public double dFdTdT() {
        return super.dFdTdT() + cpaon * dFCPAdTdT();
    }

    public double FCPA() {
        double tot = 0.0;
        double ans = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            tot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                double xai = ((ComponentSrkCPA) getComponent(i)).getXsite()[j];
                tot += (Math.log(xai) - 1.0 / 2.0 * xai + 1.0 / 2.0);
            }
            ans += getComponent(i).getNumberOfMolesInPhase() * tot;
        }
        return ans;
    }

    public double dFCPAdV() {
        return 1.0 / (2.0 * getTotalVolume()) * (1.0 - getTotalVolume() * getGcpav()) * hcpatot;
    }

    public double dFCPAdVdV() {
        return -1.0 / getTotalVolume() * dFCPAdV() + hcpatot / (2.0 * getTotalVolume()) * (-getGcpav() - getTotalVolume() * gcpavv) + getdFdVdXdXdVtotal();
    }

    public double dFCPAdVdVdV() {
        return -1.0 / getTotalVolume() * dFCPAdVdV() + 1.0 / Math.pow(getTotalVolume(), 2.0) * dFCPAdV() - hcpatot / (2.0 * Math.pow(getTotalVolume(), 2.0)) * (-getGcpav() - getTotalVolume() * gcpavv) + hcpatot / (2.0 * getTotalVolume()) * (-gcpavv - getTotalVolume() * gcpavvv - gcpavv);
    }

    public double dFCPAdT() {

        double tot = 0.0;
        double ans = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            tot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                double xai = ((ComponentSrkCPA) getComponent(i)).getXsite()[j];
                double xaidT = ((ComponentSrkCPA) getComponent(i)).getXsitedT()[j];
                tot += 1.0 / xai * xaidT - 0.5 * xaidT;// - 1.0 / 2.0 * xai + 1.0 / 2.0);
            }
            ans += getComponent(i).getNumberOfMolesInPhase() * tot;
        }
        // System.out.println("dFCPAdT " + ans);
        return ans;

        /*
        //  System.out.println
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
        for (int j = 0; j < numberOfComponents; j++) {
        double sum1 = 0;
        for (int k = 0; k < getComponent(i).getNumberOfAssociationSites(); k++) {
        for (int l = 0; l < getComponent(j).getNumberOfAssociationSites(); l++) {
        double temp = crossAccociationScheme[i][j][k][l];
        if (i == j) {
        temp = selfAccociationScheme[i][k][l];
        }
        sum1 += (((ComponentSrkCPA) getComponent(i)).getXsite()[k] * ((ComponentSrkCPA) getComponent(j)).getXsite()[l]) * temp * cpamix.calcDeltadT(k, l, i, j, this, getTemperature(), getPressure(), numberOfComponents);
        }
        }
        tot += getComponent(i).getNumberOfMolesInPhase() * getComponent(j).getNumberOfMolesInPhase() * sum1;
        }
        }

        return -1.0 * tot / 2.0 / getVolume();
         * */
    }

      public double dFCPAdTdT() {

        double tot = 0.0;
        double ans = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            tot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                double xai = ((ComponentSrkCPA) getComponent(i)).getXsite()[j];
                double xaidT = ((ComponentSrkCPA) getComponent(i)).getXsitedTdT()[j];
                tot += 1.0 / xai * xaidT - 0.5 * xaidT;// - 1.0 / 2.0 * xai + 1.0 / 2.0);
            }
            ans += getComponent(i).getNumberOfMolesInPhase() * tot;
        }
        // System.out.println("dFCPAdT " + ans);
        return ans;

    }

    public double dFCPAdTdT2() {
        //  System.out.println
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < numberOfComponents; j++) {
                double sum1 = 0;
                for (int k = 0; k < getComponent(i).getNumberOfAssociationSites(); k++) {
                    for (int l = 0; l < getComponent(j).getNumberOfAssociationSites(); l++) {
                        double temp = crossAccociationScheme[i][j][k][l];
                        double krondeltIJ = 0, krondeltKL = 0;
                        if (i == j) {
                            krondeltIJ = 1.0;
                            temp = selfAccociationScheme[i][k][l];
                        }
                        if (k == l) {
                            krondeltKL = 1.0;
                        }

                        sum1 += (((ComponentSrkCPA) getComponent(i)).getXsite()[k] * ((ComponentSrkCPA) getComponent(j)).getXsite()[l]) * temp * cpamix.calcDeltadTdT(k, l, i, j, this, getTemperature(), getPressure(), numberOfComponents);

                        double d1 = (((ComponentSrkCPA) getComponent(j)).getXsite()[l]) * temp * cpamix.calcDeltadT(k, l, i, j, this, getTemperature(), getPressure(), numberOfComponents);
                        double dQdXdX = -getComponent(i).getNumberOfMolesInPhase() / Math.pow(((ComponentSrkCPA) getComponent(i)).getXsite()[k], 2.0) * krondeltIJ * krondeltKL - getComponent(i).getNumberOfMolesInPhase() * getComponent(j).getNumberOfMolesInPhase() * temp * cpamix.calcDelta(k, l, i, j, this, getTemperature(), getPressure(), numberOfComponents);
                        double d2 = (((ComponentSrkCPA) getComponent(i)).getXsite()[k]);
                        double ans = d1 * dQdXdX * d2;
                    }
                }
                tot += getComponent(i).getNumberOfMolesInPhase() * getComponent(j).getNumberOfMolesInPhase() * sum1;
            }
        }

        return -1.0 * tot / 2.0 / getVolume();
    }

    public double dFCPAdTdV() {
        return 1.0 / (2.0 * getTotalVolume()) * (1.0 - getTotalVolume() * getGcpav()) * hcpatotdT;
    }

    public double dFCPAdTdV2() {
        //  System.out.println
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < numberOfComponents; j++) {
                double sum1 = 0;
                for (int k = 0; k < getComponent(i).getNumberOfAssociationSites(); k++) {
                    for (int l = 0; l < getComponent(j).getNumberOfAssociationSites(); l++) {
                        double temp = crossAccociationScheme[i][j][k][l];
                        if (i == j) {
                            temp = selfAccociationScheme[i][k][l];
                        }
                        sum1 += (((ComponentSrkCPA) getComponent(i)).getXsite()[k] * ((ComponentSrkCPA) getComponent(j)).getXsite()[l]) * temp * cpamix.calcDeltadTdV(k, l, i, j, this, getTemperature(), getPressure(), numberOfComponents);
                    }
                }
                tot += getComponent(i).getNumberOfMolesInPhase() * getComponent(j).getNumberOfMolesInPhase() * sum1;
            }
        }
        //System.out.println("Fcpa " + (-1.0 * tot / 2.0 / getVolume()));
        System.out.println("dFdTdV" + -1.0 * tot / 2.0 / getVolume());
        return -1.0 * tot / 2.0 / getVolume();
    }

    public double calc_hCPA() {
        double htot = 0.0;
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            htot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                htot += (1.0 - ((ComponentSrkCPA) getComponent(i)).getXsite()[j]);
            }
            tot += getComponent(i).getNumberOfMolesInPhase() * htot;
        }
        //System.out.println("tot " +tot );
        return tot;
    }

    public double calc_hCPAdT() {
        double htot = 0.0;
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            htot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                htot += -((ComponentSrkCPA) getComponent(i)).getXsitedT()[j];
            }
            tot += getComponent(i).getNumberOfMolesInPhase() * htot;
        }
        //System.out.println("tot " +tot );
        return tot;
    }

     public double calc_hCPAdTdT() {
        double htot = 0.0;
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            htot = 0.0;
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                htot += -((ComponentSrkCPA) getComponent(i)).getXsitedTdT()[j];
            }
            tot += getComponent(i).getNumberOfMolesInPhase() * htot;
        }
        //System.out.println("tot " +tot );
        return tot;
    }

    public double calc_hCPAdTdT2() {
        double htot = 0.0;
        double tot = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int k = 0; k < numberOfComponents; k++) {

                htot = 0.0;
                for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                    for (int l = 0; l < getComponent(k).getNumberOfAssociationSites(); l++) {
                        htot += ((ComponentSrkCPA) getComponent(i)).getXsite()[j] * ((ComponentSrkCPA) getComponent(k)).getXsite()[l] * getCpamix().calcDeltadTdT(j, l, i, k, this, temperature, pressure, numberOfComponents);
                    }
                }

                tot += getComponent(i).getNumberOfMolesInPhase() * getComponent(k).getNumberOfMolesInPhase() * htot;
            }
        }
        //System.out.println("tot " +tot );
        return tot / getTotalVolume();
    }

    public double calc_g() {
        double g = (2.0 - getb() / 4.0 / getMolarVolume()) / (2.0 * Math.pow(1.0 - getb() / 4.0 / getMolarVolume(), 3.0));
        return g;
    }

    public double calc_lngV() {
        double gv = 0.0, gv2 = 0.0;
        gv = -2.0 * getB() * (10.0 * getTotalVolume() - getB()) / getTotalVolume() / ((8.0 * getTotalVolume() - getB()) * (4.0 * getTotalVolume() - getB()));

        gv2 = 1.0 / (2.0 - getB() / (4.0 * getTotalVolume())) * getB() / (4.0 * Math.pow(getTotalVolume(), 2.0)) - 3.0 / (1.0 - getB() / (4.0 * getTotalVolume())) * getB() / (4.0 * Math.pow(getTotalVolume(), 2.0));
//
//        System.out.println("err gv " + (100.0-gv/gv2*100));
        //-2.0*getB()*(10.0*getTotalVolume()-getB())/getTotalVolume()/((8.0*getTotalVolume()-getB())*(4.0*getTotalVolume()-getB()));
//         System.out.println("gv " + gv);

        return gv2;
    }

    public double calc_lngVV() {
        double gvv = 0.0;
        gvv = 2.0 * (640.0 * Math.pow(getTotalVolume(), 3.0) - 216.0 * getB() * getTotalVolume() * getTotalVolume() + 24.0 * Math.pow(getB(), 2.0) * getTotalVolume() - Math.pow(getB(), 3.0)) * getB() / (getTotalVolume() * getTotalVolume()) / Math.pow(8.0 * getTotalVolume() - getB(), 2.0) / Math.pow(4.0 * getTotalVolume() - getB(), 2.0);
        return gvv;
    }

    public double calc_lngVVV() {
        double gvvv = 0.0;
        gvvv = 4.0 * (Math.pow(getB(), 5.0) + 17664.0 * Math.pow(getTotalVolume(), 4.0) * getB() - 4192.0 * Math.pow(getTotalVolume(), 3.0) * Math.pow(getB(), 2.0) + 528.0 * Math.pow(getB(), 3.0) * getTotalVolume() * getTotalVolume() - 36.0 * getTotalVolume() * Math.pow(getB(), 4.0) - 30720.0 * Math.pow(getTotalVolume(), 5.0)) * getB() / (Math.pow(getTotalVolume(), 3.0)) / Math.pow(-8.0 * getTotalVolume() + getB(), 3.0) / Math.pow(-4.0 * getTotalVolume() + getB(), 3.0);
        return gvvv;
    }

    public void setXsiteOld() {
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                ((ComponentCPAInterface) getComponent(i)).setXsiteOld(j, ((ComponentCPAInterface) getComponent(i)).getXsite()[j]);
            }
        }
    }

    public void setXsitedV(double dV) {
        dFdVdXdXdVtotal = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                double XdV = (((ComponentCPAInterface) getComponent(i)).getXsite()[j] - ((ComponentCPAInterface) getComponent(i)).getXsiteOld()[j]) / dV;
                //((ComponentCPAInterface) getComponent(i)).setXsitedV(j, XdV);
                dFdVdXdXdVtotal += XdV * ((ComponentCPAInterface) getComponent(i)).dFCPAdVdXi(j, this);
                //System.out.println("xidv " + XdV);
            }
        }
        //  System.out.println("xidv " + XdV);
    }

    public void calcXsitedV() {
        double dXdV = 0, temp = 0, temp2 = 0.0, temp22 = 0.0, temp3 = 0.0, temp4 = 0.0, temp6 = 0.0, temp7 = 0.0, temp8 = 0.0;
        double temp10 = 0.0, temp11 = 0.0, temp12 = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            temp4 = 0.0;
            temp8 = 0.0;
            temp12 = 0.0;
            for (int k = 0; k < getComponent(i).getNumberOfAssociationSites(); k++) {
                temp3 = 0.0;
                temp7 = 0.0;
                temp11 = 0.0;
                for (int j = 0; j < numberOfComponents; j++) {
                    temp2 = 0.0;
                    temp22 = 0.0;
                    temp6 = 0.0;
                    temp10 = 0.0;
                    for (int l = 0; l < getComponent(j).getNumberOfAssociationSites(); l++) {
                        temp2 += ((ComponentCPAInterface) getComponent(i)).getXsite()[k] * ((ComponentCPAInterface) getComponent(j)).getXsite()[l] * cpamix.calcDeltadV(k, l, i, j, this, temperature, pressure, numberOfComponents);
                        temp22 += ((ComponentCPAInterface) getComponent(i)).getXsite()[k] * ((ComponentCPAInterface) getComponent(j)).getXsite()[l] * cpamix.calcDelta(k, l, i, j, this, temperature, pressure, numberOfComponents) / getTotalVolume();

                        temp6 += ((ComponentCPAInterface) getComponent(i)).getXsite()[k] * ((ComponentCPAInterface) getComponent(j)).getXsite()[l] * cpamix.calcDeltadT(k, l, i, j, this, temperature, pressure, numberOfComponents);///getGcpa();

                        temp10 += ((ComponentCPAInterface) getComponent(i)).getXsite()[k] * ((ComponentCPAInterface) getComponent(j)).getXsite()[l] * cpamix.calcDeltadTdT(k, l, i, j, this, temperature, pressure, numberOfComponents);///getGcpa();
                    }
                    temp3 += (temp22 - temp2) * getComponent(j).getNumberOfMolesInPhase();
                    temp7 += temp6 * getComponent(j).getNumberOfMolesInPhase();
                    temp11 += temp10 * getComponent(j).getNumberOfMolesInPhase();
                }
                temp4 = temp3 * getComponent(i).getNumberOfMolesInPhase();
                temp8 = -temp7 * getComponent(i).getNumberOfMolesInPhase();
                temp12 = temp11 * getComponent(i).getNumberOfMolesInPhase();
                // ((ComponentCPAInterface) getComponent(i)).setXsitedV(k, temp4 / (getComponent(i).getNumberOfMolesInPhase() * (1.0 / ((ComponentCPAInterface) getComponent(i)).getXsite()[k] - 1.0 / 2.0)) / (2.0 * getTotalVolume()));
                //System.out.println("X " + ((ComponentCPAInterface) getComponent(i)).getXsite()[k] + " Xdv" + ((ComponentCPAInterface) getComponent(i)).getXsitedV()[k] + " volume " + getTotalVolume());
                ((ComponentCPAInterface) getComponent(i)).setXsitedT(k, temp8 / (getComponent(i).getNumberOfMolesInPhase() * (1.0 / ((ComponentCPAInterface) getComponent(i)).getXsite()[k] - 1.0 / 2.0)) / (2.0 * getTotalVolume()));
                ((ComponentCPAInterface) getComponent(i)).setXsitedTdT(k,temp12 / (getComponent(i).getNumberOfMolesInPhase() * (1.0 / ((ComponentCPAInterface) getComponent(i)).getXsite()[k] - 1.0 / 2.0)) / (2.0 * getTotalVolume()));
                //  System.out.println("X " + ((ComponentCPAInterface) getComponent(i)).getXsite()[k] + " Xdt" + ((ComponentCPAInterface) getComponent(i)).getXsitedT()[k] + " volume " + getTotalVolume());

            }
        }
    }

    public void calcXsitedV2() {
        double dXdV = 0, temp = 0, temp2 = 0.0, temp3 = 0.0, temp4 = 0.0;

        for (int i = 0; i < numberOfComponents; i++) {
            temp3 = 0.0;
            for (int j = 0; j < numberOfComponents; j++) {
                temp2 = 0.0;
                for (int k = 0; k < getComponent(i).getNumberOfAssociationSites(); k++) {
                    for (int l = 0; l < getComponent(j).getNumberOfAssociationSites(); l++) {
                        temp2 += ((ComponentCPAInterface) getComponent(i)).getXsite()[k] * ((ComponentCPAInterface) getComponent(j)).getXsite()[l] * cpamix.calcDeltadV(k, l, i, j, this, temperature, pressure, numberOfComponents);
                    }
                }
                temp3 += temp2 * getComponent(j).getNumberOfMolesInPhase();
                // nned to be filled in Michelsen paper p- 172
            }
            temp4 += temp3 * getComponent(i).getNumberOfMolesInPhase();
        }
        //System.out.println("sum " +temp4);//
        temp4 = (getHcpatot() - temp4) / 2.0 / getTotalVolume();

        for (int i = 0; i < numberOfComponents; i++) {
            for (int pp = 0; pp < getComponent(i).getNumberOfAssociationSites(); pp++) {
                ((ComponentCPAInterface) getComponent(i)).setXsitedV(pp, temp4 / (getComponent(i).getNumberOfMolesInPhase() * (1.0 / ((ComponentCPAInterface) getComponent(i)).getXsite()[pp] - 1.0 / 2.0)));
                //      System.out.println("X " + ((ComponentCPAInterface) getComponent(i)).getXsite()[pp] + " Xdv" + ((ComponentCPAInterface) getComponent(i)).getXsitedV()[pp] + " volume " + getTotalVolume());
            }
        }
        dFdVdXdXdVtotal = 0.0;
        for (int i = 0; i < numberOfComponents; i++) {
            for (int j = 0; j < getComponent(i).getNumberOfAssociationSites(); j++) {
                //((ComponentCPAInterface) getComponent(i)).setXsitedV(j, XdV);
                dFdVdXdXdVtotal += ((ComponentCPAInterface) getComponent(i)).getXsitedV()[j] * ((ComponentCPAInterface) getComponent(i)).dFCPAdVdXi(j, this);
            }
        }
    }

    public void calcXsitedTanalyt() { // analytically from page 172 michelsen paper
    }

    public void calcXsitedT() { // numericaly
        double dt = 1e-5, XdT = 0.0, dFCPAdTold=0.0;
        setTemperature(temperature - dt);
         do {
           // super.init(this.getNumberOfMolesInPhase(), numberOfComponents, this.getInitType(), this.getPhaseType(), this.getBeta());
        } while (!solveX());
        setXsiteOld();
        dFCPAdTold = dFCPAdT();
        updateXold = false;
        setTemperature(temperature + 2*dt);
        do {
          //  super.init(this.getNumberOfMolesInPhase(), numberOfComponents, this.getInitType(), this.getPhaseType(), this.getBeta());
        } while (!solveX());

        dFCPAdXdXdTtotal =
                0.0;
        dFCPAdTdT =
                (dFCPAdT()-dFCPAdTold)/(2*dt);
        System.out.println("dFdTdT" + dFCPAdTdT);
        for (int i = 0; i <
                numberOfComponents; i++) {
            for (int j = 0; j <
                    getComponent(i).getNumberOfAssociationSites(); j++) {
                XdT = (((ComponentCPAInterface) getComponent(i)).getXsite()[j] - ((ComponentCPAInterface) getComponent(i)).getXsiteOld()[j]) / dt;
                ((ComponentCPAInterface) getComponent(i)).setXsitedT(j, XdT);
               // System.out.println("X " +((ComponentCPAInterface) getComponent(i)).getXsite()[j]);
                //System.out.println("Xold " +((ComponentCPAInterface) getComponent(i)).getXsiteOld()[j]);
                //      System.out.println("dXdt " + XdT);
                dFCPAdXdXdTtotal +=
                        XdT * ((ComponentCPAInterface) getComponent(i)).dFCPAdXi(j, this);
            }

        }
        for (int i = 0; i <
                numberOfComponents; i++) {
            for (int j = 0; j <
                    getComponent(i).getNumberOfAssociationSites(); j++) {
                for (int k = 0; k <
                        numberOfComponents; k++) {
                    for (int j2 = 0; j2 <
                            getComponent(k).getNumberOfAssociationSites(); j2++) {
                        dFCPAdTdT += ((ComponentCPAInterface) getComponent(i)).dFCPAdXidXj(j, j2, k, this) * ((ComponentCPAInterface) getComponent(i)).getXsitedT()[j] * ((ComponentCPAInterface) getComponent(k)).getXsitedT()[j2];
                    }

                }
            }
        }
        setTemperature(temperature - dt);
        updateXold = true;
        do {
            super.init(this.getNumberOfMolesInPhase(), numberOfComponents, this.getInitType(), this.getPhaseType(), this.getBeta());
        } while (!solveX());

    }

    public double getdFdVdXdXdVtotal() {
        return dFdVdXdXdVtotal;
    }

    public boolean solveX() {
        double err = .0;
        int iter = 0;


        do {
            iter++;
            err =
                    0.0;
            for (int i = 0; i <
                    numberOfComponents; i++) {
                for (int j = 0; j <
                        getComponent(i).getNumberOfAssociationSites(); j++) {
                    double old = ((ComponentSrkCPA) getComponent(i)).getXsite()[j];
                    double neeval = getCpamix().calcXi(selfAccociationScheme, crossAccociationScheme, j, i, this, temperature, pressure, numberOfComponents);
                    ((ComponentCPAInterface) getComponent(i)).setXsite(j, neeval);
                    err +=
                            Math.abs((old - neeval) / neeval);
                }
            }
            //System.out.println("err " + err);
        } while (Math.abs(err) > 1e-12 && iter < 100);
        //System.out.println("iter " +iter);
        if (Math.abs(err) < 1e-10) {
            return true;
        } else {
            return false;
        }

    }

    /** Getter for property hcpatot.
     * @return Value of property hcpatot.
     */
    public double getHcpatot() {
        return hcpatot;
    }

    /** Setter for property hcpatot.
     * @param hcpatot New value of property hcpatot.
     */
    public void setHcpatot(double hcpatot) {
        this.hcpatot = hcpatot;
    }

    public double getGcpa() {
        return gcpa;
    }

    public double molarVolume(double pressure, double temperature, double A, double B, int phasetype) throws util.exception.IsNaNException, util.exception.TooManyIterationsException {

        double BonV = phasetype == 0 ? 2.0 / (2.0 + temperature / getPseudoCriticalTemperature()) : pressure * getB() / (numberOfMolesInPhase * temperature * R);

        if (BonV < 0) {
            BonV = 1.0e-8;
        }

        if (BonV > 1.0) {
            BonV = 1.0 - 1.0e-8;
        }

        double BonVold = BonV;
        double Btemp = 0, Dtemp = 0, h = 0, dh = 0, gvvv = 0, fvvv = 0, dhh = 0;
        double d1 = 0, d2 = 0;
        Btemp =
                getB();
        Dtemp =
                getA();
        if (Btemp < 0) {
            System.out.println("b negative in volume calc");
        }

        setMolarVolume(1.0 / BonV * Btemp / numberOfMolesInPhase);
        int iterations = 0;
        do {

            gcpa = calc_g();
            lngcpa =
                    Math.log(gcpa);
            setGcpav(calc_lngV());
            gcpavv =
                    calc_lngVV();
            gcpavvv =
                    calc_lngVVV();

            solveX();

            hcpatot =
                    calc_hCPA();


            calcXsitedV();
            iterations++;

            BonVold =
                    BonV;
            h =
                    BonV - Btemp / numberOfMolesInPhase * dFdV() - pressure * Btemp / (numberOfMolesInPhase * R * temperature);
            dh =
                    1.0 + Btemp / Math.pow(BonV, 2.0) * (Btemp / numberOfMolesInPhase * dFdVdV());
            dhh =
                    -2.0 * Btemp / Math.pow(BonV, 3.0) * (Btemp / numberOfMolesInPhase * dFdVdV()) - Math.pow(Btemp, 2.0) / Math.pow(BonV, 4.0) * (Btemp / numberOfMolesInPhase * dFdVdVdV());

            d1 =
                    -h / dh;
            d2 =
                    -dh / dhh;

            if (Math.abs(d1 / d2) <= 1.0) {
                BonV += d1 * (1.0 + 0.5 * d1 / d2);
            } else if (d1 / d2 < -1) {
                BonV += d1 * (1.0 + 0.5 * -1.0);
            } else if (d1 / d2 > 1) {
                BonV += d2;
                double hnew = h + d2 * -h / d1;
                if (Math.abs(hnew) > Math.abs(h)) {
                    BonV = phasetype == 1 ? 2.0 / (2.0 + temperature / getPseudoCriticalTemperature()) : pressure * getB() / (numberOfMolesInPhase * temperature * R);
                }

            }

            if (BonV > 1) {
                BonV = 1.0 - 1.0e-8;
                BonVold =
                        10;
            }

            if (BonV < 0) {
                BonV = 1.0e-8;
                BonVold =
                        10;
            }

            setMolarVolume(1.0 / BonV * Btemp / numberOfMolesInPhase);
            Z =
                    pressure * getMolarVolume() / (R * temperature);

            // System.out.println("Z" + Z);
        } while (Math.abs((BonV - BonVold) / BonV) > 1.0e-10 && iterations < 1001);
        //     System.out.println("Z" + Z + " iterations " + iterations + " h " + h);
        //System.out.println("pressure " + Z*R*temperature/getMolarVolume());
        //if(iterations>=100) throw new util.exception.TooManyIterationsException();
        //System.out.println("error in volume " + (-pressure+R*temperature/getMolarVolume()-R*temperature*dFdV()));// + " firstterm " + (R*temperature/molarVolume) + " second " + R*temperature*dFdV());
        if (Double.isNaN(getMolarVolume())) {
            throw new util.exception.IsNaNException();
            //System.out.println("BonV: " + BonV + " "+"  itert: " +   iterations +" " +h + " " +dh + " B " + Btemp + "  D " + Dtemp + " gv" + gV() + " fv " + fv() + " fvv" + fVV());
        }

        return getMolarVolume();
    }

    public double molarVolume2(double pressure, double temperature, double A, double B, int phase) throws util.exception.IsNaNException, util.exception.TooManyIterationsException {

        Z = phase == 0 ? 1.0 : 1.0e-5;
        setMolarVolume(Z * R * temperature / pressure);
        // super.molarVolume(pressure,temperature, A, B, phase);
        int iterations = 0;
        double err = 0.0, dErrdV = 0.0;
        double deltaV = 0;

        do {
            A = calcA(this, temperature, pressure, numberOfComponents);
            B =
                    calcB(this, temperature, pressure, numberOfComponents);



            double dFdV = dFdV(), dFdVdV = dFdVdV(), dFdVdVdV = dFdVdVdV();
            double factor1 = 1.0e0, factor2 = 1.0e0;
            err =
                    -R * temperature * dFdV + R * temperature / getMolarVolume() - pressure;

            //   System.out.println("pressure " + -R * temperature * dFdV + " " + R * temperature / getMolarVolume());
            //-pressure;
            dErrdV =
                    -R * temperature * dFdVdV - R * temperature * numberOfMolesInPhase / Math.pow(getVolume(), 2.0);

            // System.out.println("errdV " + dErrdV);
            // System.out.println("err " + err);

            deltaV =
                    -err / dErrdV;

            setMolarVolume(getMolarVolume() + deltaV / numberOfMolesInPhase);

            Z =
                    pressure * getMolarVolume() / (R * temperature);
            if (Z < 0) {
                Z = 1e-6;
                setMolarVolume(Z * R * temperature / pressure);
            }

//     System.out.println("Z " + Z);
        } while (Math.abs(err) > 1.0e-8 || iterations < 100);
        //System.out.println("Z " + Z);
        return getMolarVolume();
    }

    public double getGcpav() {
        return gcpav;
    }

    public void setGcpav(double gcpav) {
        this.gcpav = gcpav;
    }

    public CPAMixingInterface getCpamix() {
        return cpamix;
    }

    public double calcPressure() {
        gcpa = calc_g();
        lngcpa =
                Math.log(gcpa);
        setGcpav(calc_lngV());
        gcpavv =
                calc_lngVV();
        gcpavvv =
                calc_lngVVV();
        solveX();

        hcpatot =
                calc_hCPA();
        return super.calcPressure();
    }

    public int getCrossAssosiationScheme(int comp1, int comp2, int site1, int site2) {
        if (comp1 == comp2) {
            return selfAccociationScheme[comp1][site1][site2];
        } else {
            return crossAccociationScheme[comp1][comp2][site1][site2];
        }
    }
}
