*   Classes Donnan

*

*   Class Donnan contains the primary methods for

*   calculating a Donnan Potential between two

*   partitions of different electrical permittivity

*   (dielectric constant) between which any number of

*   ionic species may be partitioned and, in one

*   partition only, may bind to a neutral ionophore.

*

*   Class Donnan requires Class DonnanMinim and Class DonnanConcn

*   that implement interfaces to the required minimisation methods

*

*   WRITTEN BY: Michael Thomas Flanagan

*

*   DATE:           November 2004

*   LAST UPDATE:    1 December 2004, 5-7 July 2008

*

*   DOCUMENTATION:

*   See Michael Thomas Flanagan's JAVA library on-line web page:

*   http://www.ee.ucl.ac.uk/~mflanaga/java/Donnan.html

*   http://www.ee.ucl.ac.uk/~mflanaga/java/

*

*   Copyright (c) November 2004

*

*   PERMISSION TO COPY:

*   Permission to use, copy and modify this software and its documentation for

*   NON-COMMERCIAL purposes is granted, without fee, provided that an acknowledgement

*   to the author, Michael Thomas Flanagan at www.ee.ucl.ac.uk/~mflanaga, appears in all copies.

*

*   Dr Michael Thomas Flanagan makes no representations about the suitability

*   or fitness of the software for any or for a particular purpose.

*   Michael Thomas Flanagan shall not be liable for any damages suffered

*   as a result of using, modifying or distributing this software or its derivatives.

*

***************************************************************************************/

    private ArrayList<Object> arrayl = new ArrayList<Object>();      //  vector holding:

                                            //  element 0:  name of first ion

                                            //  element 1:  initial concentration of first ion in partition A (enter as M)

                                            //  element 2:  initial concentration of first ion in partition B (enter as M)

                                            //  element 3:  association constant with ionophore (enter as M^-1)

                                            //  element 4:  radius of first ion (metres)

                                            //  element 5:  charge of first ion

                                            //  element 6:  partition coefficient if provided by the user

                                            //              -1 if not; partition coefficient to be calculated later

                                            //  element 7:  name of second ion

                                            //  etc

    private int numOfIons = 0;              //  number of ionic species

    private int numOfAnions = 0;            //  number of anionic species

    private int numOfCations = 0;           //  number of cationic species

    private String[] ionNames = null;       //  names of the ions

    private double[] concnA = null;         //  concentration of ions in partition A (moles per cubic metre)

    private double[] concnB = null;         //  concentration of ions in partition B (moles per cubic metre)

    private double[] molesT = null;         //  total moles of an ion

    private double[] complex = null;        //  concentration of complex in partition B (moles per cubic metre)

    private double[] excessConcnA = null;   //  excess concentration of ions in partition A in interfacial charge (moles per cubic metre)

    private double[] excessConcnB = null;   //  excess concentration of ions in partition B in interfacial charge (moles per cubic metre)

    private double[] excessComplex = null;  //  excess concentration of complex in partition B in interfacial charge (moles per cubic metre)

    private int[] indexC = null;            //  ion index of ionic species with non zero total concentrations

    private int nonZeroConcns = 0;          //  number of ionic species with a non-zero total concentration

    private double[] assocConsts = null;    //  association constants of the ions with the ionophore (cubic metres per mole)

    private int[] indexK = null;            //  ion index of ionic species with non zero assocConsts

    private int nonZeroAssocK = 0;          //  number of ionic species with affinity for the ionophore

    private double[] radii = null;          //  radii of ions

    private boolean radiusType = true;      //  if = true - hydrated radii are taken from Class IonicRadii

                                            //  if = false - bare radii are taken from Class IonicRadii

    private double[] charges = null;        //  charge of ions

    private double tol = 1e-6;              //  fractional tolerance in checking for overall charge neutrality

                                            //  when multiplied by total concentration of species of lowest concentration

                                            //  gives limit for considering overall neurtality achieved

    private String ionophore = "ionophore"; //  name ionophore for output texts

    private double ionophoreConcn = 0.0D;   //  ionophore concentration - entered as molar, converted to moles per cubic metre

    private double freeIonophoreConcn = 0.0D;   // ionophore concentration - all complex concentrations

    private double ionophoreRad = 0.0D;     //  ionophore radius (metres)

    private boolean ionophoreSet = false;   //  = true when ionophore concentration and radius have been set

    private double volumeA = 0.0D;          //  volume of partition A (cubic metres)

    private double volumeB = 0.0D;          //  volume of partition B (cubic metres)

    private double interfacialArea = 0.0D;  //  onterfacial area between Partitions A and B(square metres)

    private boolean volumeSet = false;      //  = true when volumeA and volumeB have been set

    private double epsilonA = 0.0D;         //  relative electrical permittivity of partition A

    private double epsilonB = 0.0D;         //  relative electrical permittivity of partition B

    private double epsilonSternA = 0.0D;    //  relative electrical permittivity of partition A Stern layer

    private double epsilonSternB = 0.0D;    //  relative electrical permittivity of partition B SternLayer

    private boolean epsilonSet = false;     //  = true when epsilonA and epsilonB have been set

    private double temp = 25.0-Fmath.T_ABS; //  Temperature (degrees Kelvin) [Enter temperature in degrees Celsius]

    private boolean tempSet = false;        //  = true when temperature has been set

    private double[] deltaMu0 = null;       //  difference in Born charging energy (partitionB - partitionA)

    private double[] partCoeff = null;      //  partition coefficients in absence of inter-partition, e.g. Donnan, potential

    private double[] partCoeffPot = null;   //  partition coefficients in presence of inter-partition, e.g. Donnan, potential

    private boolean[] indexPC = null;       //  partition coefficient (pc) index

                                            //  = true if pc provided by user

                                            //  = false if pc calculated from Born charging

    private double donnanPotential = 0.0D;  //  Donanan potential with respect to partition A (psiB - psiA) [volts]

    private double diffPotentialA = 0.0D;   //  Double layer potential difference - compartment A [volts]

    private double diffPotentialB = 0.0D;   //  Double layer potential difference - compartment B [volts]

    private double sternPotential = 0.0D;   //  Stern potential difference  [volts]

    private double estimate = 0.0D;         //  Initial Estimate of the Donanan potential with respect to partition A (psiB - psiA) [volts]

    private double step = 0.0D;             //  Initial step size in the estimation of the Donanan potential with respect to partition A (psiB - psiA) [volts]

    private double tolerance = 1.0e-20;     //  tolerance in exiting estimation of the Donanan potential with respect to partition A (psiB - psiA)

    private int nMaxIter = 10000;           //  maximum number of iterations allowed in estimation of the Donanan potential with respect to partition A (psiB - psiA)

    private int numIterations = 0;          //  number of iterations taken

    private double minimum = 1.0e300;       //  value of function to be minimised at the minimum

    private double sternCap = 0.0D;         //  Stern layer capacitance [F]

    private double diffCapA = 0.0D;         //  diffuse double layer capacitance in partition A [F]

    private double diffCapB = 0.0D;         //  diffuse double layer capacitance in partition B [F]

    private double donnanCap = 0.0D;        //  total interfacial capacitance [F]

    private double sternDeltaA = 0.0D;      //  Stern layer thickness in partition A [m]

    private double sternDeltaB = 0.0D;      //  Stern layer thickness in partition B [m]

    private double chargeValue = 0;         //  Absolute value of the charge valency if all are the same, i.e. symmetric electrolytes of the same valency

    private boolean chargeSame = true;      //  = false if chargeValue not the same for all ions

    private double interfacialChargeDensity = 0.0D;// interface Charge Density [C per square metre]

    private double interfacialCharge = 0.0D;       // interface Charge [C]

    private boolean includeIc = true;       //  = true - interface charge included in the calculation of the Donnan potential

                                            //  = false - interface charge ignored in the calculation of the Donnan Poptential

    private double[] ratioA = null;         //  ratio of excess to bulk concentrations - compartment A

    private double[] ratioB = null;         //  ratio of excess to bulk concentrations - compartment B

    private double[] ratioC = null;         //  ratio of excess to bulk concentrations - complex

    private double recipKappaA = 0.0D;      //  Debye length - compartment A

    private double recipKappaB = 0.0D;      //  Debye length - compartment B

    // Constructor

    public Donnan(){

    }

    // Method for setting ionic radii taken from Class IonicRadii to hydrated radii

    // This is the default option

    public void setHydratedRadii(){

        this.radiusType = true;

    }

    // Method for setting ionic radii taken from Class IonicRadii to bare radii

    public void setBareRadii(){

        this.radiusType = false;

    }

    // Method for setting Donnan Potential calculatiom method to ignore interfacial charge

    public void ignoreInterfaceCharge(){

        this.includeIc = false;

    }

     // Method for setting Donnan Potential calculatiom method to include interfacial charge

     // This is the default option

    public void includeInterfaceCharge(){

        this.includeIc = true;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient to be calculated from Born charging equation

    // Concentrations - Molar,  assocK - M^-1, radius - metres, charge - valency e.g. +1

    public void setIon(String ion, double concnA, double concnB, double assocK, double radius, int charge){

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(assocK));

        if(assocK!=0.0D)this.nonZeroAssocK++;

        this.arrayl.add(new Double(radius));

        this.arrayl.add(new Integer(charge));

        this.arrayl.add(new Double(-1.0D));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient supplied by user

    // Concentrations - Molar,  assocK - M^-1, radius - metres, charge - valency e.g. +1

    public void setIon(double partCoeff, String ion, double concnA, double concnB, double assocK, double radius, int charge){

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(assocK));

        if(assocK!=0.0D)this.nonZeroAssocK++;

        this.arrayl.add(new Double(radius));

        this.arrayl.add(new Integer(charge));

        this.arrayl.add(new Double(partCoeff));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient to be calculated from Born charging equation

    // Concentrations - Molar,  radius - metres, charge - valency e.g. +1

    // association constant = 0.0D

    public void setIon(String ion, double concnA, double concnB, double radius, int charge){

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(0.0D));

        this.arrayl.add(new Double(radius));

        this.arrayl.add(new Integer(charge));

        this.arrayl.add(new Double(-1.0D));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient to be supplied by user

    // Concentrations - Molar,  radius - metres, charge - valency e.g. +1

    // association constant = 0.0D

    public void setIon(double partCoeff, String ion, double concnA, double concnB, double radius, int charge){

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(0.0D));

        this.arrayl.add(new Double(radius));

        this.arrayl.add(new Integer(charge));

        this.arrayl.add(new Double (partCoeff));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient to be calculated from Born charging equation

    // default radii and charge taken from class IonicRadii

    // if radii not in Ionic Radii, Donnan potential calculated with interface charge neglected

    // Concentrations - Molar

    public void setIon(String ion, double concnA, double concnB, double assocK){

        IonicRadii ir = new IonicRadii();

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(assocK));

        if(assocK!=0.0D)this.nonZeroAssocK++;

        double rad = 0.0D;

        if(this.radiusType){

            rad = ir.hydratedRadius(ion);

        }

        else{

            rad = ir.radius(ion);

        }

        if(rad==0.0D){

            String mess1 = ion + " radius is not in the IonicRadii list\n";

            String mess2 = "Please enter radius in metres\n";

            rad = Db.readDouble(mess1+mess2);

        }

        this.arrayl.add(new Double(rad));

        int charg = 0;

        charg = ir.charge(ion);

        if(charg==0){

            String mess1 = ion + " charge is not in the IonicRadii list\n";

            String mess2 = "Please enter charge, e.g +2";

            charg = Db.readInt(mess1+mess2);

        }

        this.arrayl.add(new Integer(charg));

        this.arrayl.add(new Double(-1.0D));

        this.numOfIons++;

    }

   // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficient to be supplied by user

    // default radii and charge taken from class IonicRadii

    // Concentrations - Molar

    public void setIon(double partCoeff, String ion, double concnA, double concnB, double assocK){

        IonicRadii ir = new IonicRadii();

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(assocK));

        if(assocK!=0.0D)this.nonZeroAssocK++;

        double rad = 0.0D;

        if(this.includeIc){

            if(this.radiusType){

                rad = ir.hydratedRadius(ion);

            }

            else{

                rad = ir.radius(ion);

            }

            if(rad==0.0D){

                String mess1 = ion + " radius is not in the IonicRadii list\n";

                String mess2 = "Please enter radius in metres\n";

                String mess3 = "Enter 0.0 if you wish interfacial charge to be neglected";

                rad = Db.readDouble(mess1+mess2+mess3);

                if(rad==0.0D)this.includeIc = false;

            }

        }

        this.arrayl.add(new Double(rad));

        int charg = 0;

        charg = ir.charge(ion);

        if(charg==0){

            String mess1 = ion + " charge is not in the IonicRadii list\n";

            String mess2 = "Please enter charge, e.g +2";

            charg = Db.readInt(mess1+mess2);

        }

        this.arrayl.add(new Integer(charg));

        this.arrayl.add(new Double(partCoeff));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficients to be calculated from Born charging equation

    // default radii and charge taken from class IonicRadii

    // association constant = 0.0D

    // Concentrations - Molar

    public void setIon(String ion, double concnA, double concnB){

        IonicRadii ir = new IonicRadii();

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(0.0D));

        double rad = 0.0D;

        if(this.radiusType){

            rad = ir.hydratedRadius(ion);

        }

        else{

            rad = ir.radius(ion);

        }

        if(rad==0.0D){

            String mess1 = ion + " radius is not in the IonicRadii list\n";

            String mess2 = "Please enter radius in metres\n";

            rad = Db.readDouble(mess1+mess2);

            if(rad==0.0D)this.includeIc = false;

        }

        this.arrayl.add(new Double(rad));

        int charg = 0;

        charg = ir.charge(ion);

        if(charg==0){

            String mess1 = ion + " charge is not in the IonicRadii list\n";

            String mess2 = "Please enter charge, e.g +2";

            charg = Db.readInt(mess1+mess2);

        }

        this.arrayl.add(new Integer(charg));

        this.arrayl.add(new Double(-1.0D));

        this.numOfIons++;

    }

    // Method to add an ionic species to the  Donnan Equilibrium

    // Partition coefficients to be supplied by the user

    // default radii and charge taken from class IonicRadii

    // association constant = 0.0D

    // Concentrations - Molar

    public void setIon(double partCoeff, String ion, double concnA, double concnB){

        IonicRadii ir = new IonicRadii();

        this.arrayl.add(ion);

        this.arrayl.add(new Double(concnA));

        this.arrayl.add(new Double(concnB));

        if(concnA>0.0D || concnB>0.0)this.nonZeroConcns++;

        this.arrayl.add(new Double(0.0D));

        double rad = 0.0D;

        if(this.includeIc){

            if(this.radiusType){

                rad = ir.hydratedRadius(ion);

            }

            else{

                rad = ir.radius(ion);

            }

            if(rad==0.0D){

                String mess1 = ion + " radius is not in the IonicRadii list\n";

                String mess2 = "Please enter radius in metres\n";

                String mess3 = "Enter 0.0 if you wish interfacial charge to be neglected";

                rad = Db.readDouble(mess1+mess2+mess3);

                if(rad==0.0D)this.includeIc = false;

            }

        }

        this.arrayl.add(new Double(rad));

        int charg = 0;

        charg = ir.charge(ion);

        if(charg==0){

            String mess1 = ion + " charge is not in the IonicRadii list\n";

            String mess2 = "Please enter charge, e.g +2";

            charg = Db.readInt(mess1+mess2);

        }

        this.arrayl.add(new Integer(charg));

        this.arrayl.add(new Double(partCoeff));

        this.numOfIons++;

    }

    // Method to add an ionophore to the partition B

    // Concentration = Molar, radius = metres

    public void setIonophore(double concn, double radius){

        this.ionophoreConcn = concn*1e3;

        this.ionophoreRad = radius;

        this.ionophoreSet = true;

    }

    // Method to add an ionophore to the partition B

    // Concentration = Molar, radius = metres

    // Ionophore name can be added for output text

    public void setIonophore(String ionophore, double concn, double radius){

        this.ionophore = ionophore;

        this.ionophoreConcn = concn*1e3;

        this.ionophoreRad = radius;

        this.ionophoreSet = true;

    }

    // Method to add an ionophore to the partition B

    // No radius added - calculation of Donnan potential will neglect interface charge

    // Concentration = Molar

    // Ionophore name can be added for output text

    public void setIonophore(String ionophore, double concn){

        this.ionophore = ionophore;

        this.ionophoreConcn = concn*1e3;

        this.includeIc = false;

        this.ionophoreSet = true;

    }

    // Method to add an ionophore to the partition B

    // No radius added - calculation of Donnan potential will neglect interface charge

    // Concentration = Molar

    public void setIonophore(double concn){

        this.ionophoreConcn = concn*1e3;

        this.includeIc = false;

        this.ionophoreSet = true;

    }

    // Method to set partition volumes (m^3) and interfacial area

    public void setVolumes(double volA, double volB, double area){

        this.volumeA = volA;

        this.volumeB = volB;

        this.interfacialArea = area;

        this.volumeSet = true;

    }

    // Method to set partition volumes (m^3)

    // No interfacial area entered - calculation of Donnan potential will neglect interface charge

    public void setVolumes(double volA, double volB){

        this.volumeA = volA;

        this.volumeB = volB;

        this.includeIc = false;

        this.volumeSet = true;

    }

    // Method to set partition relative permittivities

    public void setRelPerm(double epsA, double epsB, double epsSternA, double epsSternB){

        this.epsilonA = epsA;

        this.epsilonB = epsB;

        this.epsilonSternA = epsSternA;

        this.epsilonSternB = epsSternB;

        this.epsilonSet = true;

    }

    // Method to set partition relative permittivities

    // No Stern layer permittivities included - interface charge cannot be calculated

    public void setRelPerm(double epsA, double epsB){

        this.epsilonA = epsA;

        this.epsilonB = epsB;

        this.includeIc = false;

        this.epsilonSet = true;

    }

    // Method to set temperature (enter as degrees Celsius)

    public void setTemp(double temp){

        this.temp = temp - Fmath.T_ABS;

        this.tempSet = true;

    }

    // Method to set initial estimate of the Donnan potential (V)

    public void setEstimate(double pot){

        this.estimate = pot;

    }

    // Method to set initial step size of the Donnan potential Estimations (V)

    public void setStep(double step){

        this.step = step;

    }

    // Method to set tolerance in exiting the Donnan potential calculation

    public void setTolerance(double tol){

        this.tolerance = tol;

    }

    // Method to set maximum number of iterations in calculating the Donnan potential calculation

    public void setMaxIterations(int nMax){

        this.nMaxIter = nMax;

    }

    // Method to get the Donnan potential [volts]

    public double getDonnanPotential(){

        return this.donnanPotential;

    }

    // Method to get the diffuse double layer potential in partition A [volts]

    public double getDiffuseLayerPotentialA(){

        return this.diffPotentialA;

    }

    // Method to get the diffuse double layer potential in partition B [volts]

    public double getDiffuseLayerPotentialB(){

        return this.diffPotentialB;

    }

    // Method to get the Stern layer potential [volts]

    public double getSternLayerPotential(){

        return this.sternPotential;

    }

    // Method to get the bulk concentrations of the ionic species, in partition A, at the minimum [M]

    public double[] getConcnA(){

        double[] concn = this.concnA.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the bulk concentrations of the ionic species, in partition B, at the minimum [M]

    public double[] getConcnB(){

        double[] concn = this.concnB.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the bulk concentrations of the complex species, in partition B, at the minimum [M]

    public double[] getComplex(){

        double[] concn = this.complex.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the excess concentrations of the ionic species, in partition A, at the minimum [M]

    public double[] getExcessConcnA(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getExcessConcnA\nThe values of the excess concentrations have not been calculated\nzeros returned");

        }

        double[] concn = this.excessConcnA.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the excessConcentrations of the ionic species, in partition B, at the minimum [M]

    public double[] getExcessConcnB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getExcessConcnA\nThe values of the excess concentrations have not been calculated\nzeros returned");

        }

        double[] concn = this.excessConcnB.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the excess concentrations of the complex species, in partition B, at the minimum [M]

    public double[] getExcessComplex(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getExcessConcnA\nThe values of the excess concentrations have not been calculated\nzeros returned");

        }

        double[] concn = this.excessComplex.clone();

        for(int i=0; i<this.numOfIons; i++)concn[i] *= 1e-3;

        return concn;

    }

    // Method to get the ratio of excess concentrations over bulk concentrations of the ionic in partition A

    public double[] getRatioA(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getRatioA\nThe values of the excess to bulk concentrations have not been calculated\nzeros returned");

        }

      return this.ratioA;

    }

    // Method to get the ratio of excess concentrations over bulk concentrations of the ionic but not complex species, in partition B

    public double[] getRatioB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getRatioB\nThe values of the excess to bulk concentrations have not been calculated\nzeros returned");

        }

      return this.ratioB;

    }

    // Method to get the ratio of excess concentrations over bulk concentrations of the complex species

    public double[] getRatioComplex(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getRatioComplex\nThe values of the excess to bulk concentrations have not been calculated\nzeros returned");

        }

      return this.ratioC;

    }

    // Method to get the partition coefficients at equilibrium

    public double[] getPartitionCoefficients(){

      return this.partCoeffPot;

    }

    // Method to get the partition coefficients in the absence of a Donnan Potential

    public double[] getPartitionCoefficientsZero(){

      return this.partCoeff;

    }

    // Method to get the difference in Born Charging energy between partitions A and B

    public double[] getDeltaMu0(){

        return this.deltaMu0;

    }

     // Method to get the total interfacial charge

    public double getInterfaceCharge(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getInterfaceCharge\nThe value of the interface charge has not been calculated\nzero returned");

        }

        return this.interfacialCharge;

    }

    // Method to get the interfacial charge density

    public double getInterfaceChargeDensity(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getInterfaceChargeDensity\nThe value of the interface charge density has not been calculated\nzero returned");

        }

        return this.interfacialCharge;

    }

    // Method to get the Stern Capacitance

    public double getSternCapacitance(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getSternCapacitance\nThe value of the Stern capacitance has not been calculated\nzero returned");

        }

        return this.sternCap*this.interfacialArea;

    }

    // Method to get the diffuse layer Capacitance in Partition A

    public double getDiffuseLayerCapacitanceA(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getDiffuseLayerCapacitanceA\nThe values of the diffuse layer capacitances have not been calculated\nzero returned");

        }

        return this.diffCapA*this.interfacialArea;

    }

    // Method to get the diffuse layer Capacitance in Partition B

    public double getDiffuseLayerCapacitanceB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getDiffuseLayerCapacitanceB\nThe values of the diffuse layer capacitances have not been calculated\nzero returned");

        }

        return this.diffCapB*this.interfacialArea;

    }

    // Method to get the Donnan Capacitance

    public double getDonnanCapacitanceB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getDonnanCapacitance\nThe value of the Donnan capacitance has not been calculated\nzero returned");

        }

        return this.donnanCap*this.interfacialArea;

    }

    // Method to get the Stern layer thickness in Partition A

    public double getSternThicknessA(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getSternThicknessA\nThe values of the Stern layer thicknesses have not been calculated\nzero returned");

        }

        return this.sternDeltaA;

    }

    // Method to get the Stern layer thickness in Partition B

    public double getSternThicknessB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getSternThicknessB\nThe values of the Stern layer thicknesses have not been calculated\nzero returned");

        }

        return this.sternDeltaB;

    }

    // Method to get the Debye length in Partition A

    public double getDebyeLengthA(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getDebyeLengthA\nThe values of the Debye lengths have not been calculated\nzero returned");

        }

        return this.recipKappaA;

    }

    // Method to get the Debye length in Partition A

    public double getDebyeLengthB(){

        if(!this.includeIc){

            System.out.println("Class: Donnan\nMethod: getDebyeLengthB\nThe values of the Debye lengths have not been calculated\nzero returned");

        }

        return this.recipKappaB;

    }

    // Method to get the minimised function (square of charge in B) value at the minimum

    public double getMinimum(){

        return this.minimum;

    }

    // Method to calculate the Donnan potential (V)

    public double calcPotential(){

        // check all information needed is present and check overall charge neutrality

        unpack();

        // Repack to eliminate zero concentration ions

        double numOfIonsHold = this.numOfIons;

        double[] assocConstshold = null;

        double[] radiihold = null;

        double[] chargeshold = null;

        double[] deltaMu0hold = null;

        double[] partCoeffhold = null;

        if(this.nonZeroConcns<this.numOfIons){

            assocConstshold = this.assocConsts.clone();

            radiihold = this.radii.clone();

            chargeshold = this.charges.clone();

            deltaMu0hold = this.deltaMu0.clone();

            partCoeffhold = this.partCoeff.clone();

            boolean test = true;

            int jj = 0;

            while(test){

                if(indexC[jj]==0){

                    for(int k=jj+1; k<this.numOfIons; k++){

                        this.concnA[k-1] = this.concnA[k];

                        this.concnB[k-1] = this.concnB[k];

                        this.complex[k-1] = this.complex[k];

                        this.molesT[k-1] = this.molesT[k];

                        this.assocConsts[k-1] = this.assocConsts[k];

                        this.radii[k-1] = this.radii[k];

                        this.charges[k-1] = this.charges[k];

                        this.deltaMu0[k-1] = this.deltaMu0[k];

                        this.partCoeff[k-1] = this.partCoeff[k];

                    }

                    this.numOfIons--;

                }

                else{

                    jj++;

                }

                if(this.numOfIons==this.nonZeroConcns)test=false;

            }

        }

        // Obtain Donnan potential by miminimising

        // the square of the net charge in partition B

        // Check if interface charge is to be include

        // if it is - perform calculation first without interface charge

        // to obtain workable initial estimates for calculation with interface charge included

        boolean includeIcHold = false;

        if(this.includeIc){

            includeIcHold = true;

            this.includeIc = false;

        }

        //Create instance of Minimisation

        Minimisation minPot = new Minimisation();

        // Create instace of class holding function to be minimised

        DonnanMinim functD = new DonnanMinim(this.numOfIons);

        // Initialise function functD variable

        functD.numOfIons        = this.numOfIons;

        functD.concnA           = this.concnA;

        functD.concnB           = this.concnB;

        functD.molesT           = this.molesT;

        functD.complex          = this.complex;

        functD.excessConcnA     = this.excessConcnA;

        functD.excessConcnB     = this.excessConcnB;

        functD.excessComplex    = this.excessComplex;

        functD.assocConsts      = this.assocConsts;

        functD.indexK           = this.indexK;

        functD.nonZeroAssocK    = this.nonZeroAssocK;

        functD.radii            = this.radii;

        functD.charges          = this.charges;

        functD.ionophoreConcn   = this.ionophoreConcn;

        functD.ionophoreRad     = this.ionophoreRad;

        functD.volumeA          = this.volumeA;

        functD.volumeB          = this.volumeB;

        functD.interfacialArea  = this.interfacialArea;

        functD.epsilonA         = this.epsilonA;

        functD.epsilonB         = this.epsilonB;

        functD.epsilonSternA    = this.epsilonSternA;

        functD.epsilonSternB    = this.epsilonSternB;

        functD.temp             = this.temp;

        functD.partCoeff        = this.partCoeff;

        functD.partCoeffPot     = this.partCoeffPot;

        functD.sternCap         = this.sternCap;

        functD.sternDeltaA      = this.sternDeltaA;

        functD.sternDeltaB      = this.sternDeltaB;

        functD.chargeValue      = this.chargeValue;

        functD.chargeSame       = this.chargeSame;

        functD.interfacialCharge = this.interfacialCharge;

        functD.interfacialChargeDensity = this.interfacialChargeDensity;

        functD.includeIc        = this.includeIc;

        // initial estimate

        double[] start = {this.estimate};

        // initial step size

        double[] step = {this.step};

        // Nelder and Mead minimisation procedure

        minPot.nelderMead(functD, start, step, this.tolerance, this.nMaxIter);

        // get values of the Donnan potential at minimum

        double[] param = minPot.getParamValues();

        this.donnanPotential = param[0];

        // Repeat with interface charge included if this is required

        if(includeIcHold){

            this.includeIc = true;

            // initial estimate

            start[0] = this.donnanPotential;

            // initial step size

            step[0] = this.step;

            // Nelder and Mead minimisation procedure

            minPot.nelderMead(functD, start, step, this.tolerance, this.nMaxIter);

            // get values of the Donnan potential at minimum

            param = minPot.getParamValues();

            this.donnanPotential = param[0];

        }

        // Calculate concentartions at the minimum

        ionConcns(this.donnanPotential);

        // Repack if zero concentration ions were present

        if(this.nonZeroConcns!=numOfIonsHold){

            boolean test = true;

            int jj = 0;

            while(test){

                if(indexC[jj]==0){

                    for(int k=jj; k<this.numOfIons; k++){

                        this.concnA[k+1] = this.concnA[k];

                        this.concnB[k+1] = this.concnB[k];

                        this.complex[k+1] = this.complex[k];

                        this.excessConcnA[k+1] = this.excessConcnA[k];

                        this.excessConcnB[k+1] = this.excessConcnB[k];

                        this.excessComplex[k+1] = this.excessComplex[k];

                        this.molesT[k+1] = this.molesT[k];

                        this.assocConsts[k+1] = this.assocConsts[k];

                        this.radii[k+1] = this.radii[k];

                        this.charges[k+1] = this.charges[k];

                        this.deltaMu0[k+1] = this.deltaMu0[k];

                        this.partCoeff[k+1] = this.partCoeff[k];

                        this.partCoeffPot[k+1] = this.partCoeffPot[k];

                    }

                    this.numOfIons++;

                    this.concnA[jj] = 0.0D;

                    this.concnB[jj] = 0.0D;

                    this.complex[jj] = 0.0D;

                    this.excessConcnA[jj] = 0.0D;

                    this.excessConcnB[jj] = 0.0D;

                    this.excessComplex[jj] = 0.0D;

                    this.molesT[jj] = 0.0D;

                    this.assocConsts[jj] = assocConstshold[jj];

                    this.radii[jj] = radiihold[jj];

                    this.charges[jj] = chargeshold[jj];

                    this.deltaMu0[jj] = deltaMu0hold[jj];

                    this.partCoeff[jj] = partCoeffhold[jj];

                }

                else{

                    jj++;

                }

                if(this.numOfIons==this.nonZeroConcns)test=false;

            }

        }

        // get the minimum value

        this.minimum = minPot.getMinimum();

        // get the number of iterations

        this.numIterations = minPot.getNiter();

        if(this.includeIc){

            // Calculate excess to bulk ratios

            for(int i=0; i<this.numOfIons; i++){

                this.ratioA[i] = this.excessConcnA[i]/this.concnA[i];

                this.ratioB[i] = this.excessConcnB[i]/this.concnB[i];

                this.ratioC[i] = this.excessComplex[i]/this.complex[i];

            }

            // Calculate overall and diffuse layer capacitances

            this.diffCapA = this.interfacialCharge/this.diffPotentialA;

            this.diffCapB = this.interfacialCharge/this.diffPotentialB;

            this.donnanCap = this.interfacialCharge/this.donnanPotential;

            // Calculate Debye lengths

            double preterm = 2.0D*Fmath.square(Fmath.Q_ELECTRON)*Fmath.N_AVAGADRO/(Fmath.EPSILON_0*Fmath.K_BOLTZMANN*this.temp);

            double pretermA = preterm/this.epsilonA;

            double pretermB = preterm/this.epsilonB;

            this.recipKappaA = 0.0;

            this.recipKappaB = 0.0;

            for(int i=0; i<this.numOfIons; i++){

                this.recipKappaA += this.concnA[i]*charges[i]*charges[i];

                this.recipKappaB += (this.concnB[i] + this.complex[i])*charges[i]*charges[i];

            }

            this.recipKappaA = 1.0D/Math.sqrt(this.recipKappaA*pretermA);

            this.recipKappaB = 1.0D/Math.sqrt(this.recipKappaB*pretermB);

            // Calculate deltaMu0s for user supplied partition coefficients

            // and scale all deltaMu0s to per mole

            for(int ii=0; ii<this.numOfIons; ii++){

                if(this.indexPC[ii]){

                    this.deltaMu0[ii] = Math.log(this.partCoeff[ii])*(Fmath.N_AVAGADRO*Fmath.K_BOLTZMANN)*this.temp;

                }

                else{

                    this.deltaMu0[ii] *= Fmath.N_AVAGADRO;

                }

            }

            // Calculate partition coefficients at equilibrium

            for(int ii=0; ii<this.numOfIons; ii++){

                this.partCoeffPot[ii] = this.partCoeff[ii]*Math.exp((-this.donnanPotential*this.charges[ii]*Fmath.Q_ELECTRON)/(Fmath.K_BOLTZMANN*this.temp));

            }

        }

        return this.donnanPotential;

    }

    // Unpacks the ion storage ArrayList and fills the appropriate arrays

    // Checks all relevant data has been entered

    private void unpack(){

        if(!this.volumeSet)throw new IllegalArgumentException("The volumes of the partitions have not been set");

        if(this.numOfIons==0)throw new IllegalArgumentException("No ions have been entered");

        if(this.nonZeroConcns==0)throw new IllegalArgumentException("No non-zero ionic concentrations have been entered");

        if(!this.epsilonSet)throw new IllegalArgumentException("The relative permittivities of the partitions have not been set");

        if(!this.tempSet)System.out.println("The temperature has not been entered\na value of 25 degrees Celsius has been used");

        if(!this.ionophoreSet){

            System.out.println("The ionophore has not been entered\na concentration value of zero has been used");

            if(this.includeIc){

                this.includeIc=false;

                System.out.println("and the interface charge option has been set to neglect the interface charge");

            }

        }

        // change concentrations to moles per cubic metre

        // fill primitive data arrays

        // calculate total moles

        // check if the electrolyte is charge symmetric

        this.ionNames = new String[this.numOfIons];

        this.concnA = new double[this.numOfIons];

        this.concnB = new double[this.numOfIons];

        this.molesT = new double[this.numOfIons];

        this.complex = new double[this.numOfIons];

        this.excessConcnA = new double[this.numOfIons];

        this.excessConcnB = new double[this.numOfIons];

        this.excessComplex = new double[this.numOfIons];

        this.ratioA = new double[this.numOfIons];

        this.ratioB = new double[this.numOfIons];

        this.ratioC = new double[this.numOfIons];

        this.assocConsts = new double[this.numOfIons];

        this.radii = new double[this.numOfIons];

        this.charges = new double[this.numOfIons];

        this.deltaMu0 = new double[this.numOfIons];

        this.partCoeff = new double[this.numOfIons];

        this.partCoeffPot = new double[this.numOfIons];

        this.indexK = new int[this.nonZeroAssocK];

        this.indexC = new int[this.numOfIons];

        this.indexPC = new boolean[this.numOfIons];

        Double hold = null;

        Integer holi = null;

        int ii = 0;

        this.chargeValue = 0;

        this.chargeSame = true;

        for(int i=0; i<numOfIons; i++){

            // ion name

            this.ionNames[i]= (String)this.arrayl.get(0+i*7);

            // concentration in compartment A

            hold            = (Double)this.arrayl.get(1+i*7);

            this.concnA[i]  = hold.doubleValue()*1e3;

            hold            = (Double)this.arrayl.get(2+i*7);

            // concentration in compartment B

            this.concnB[i]  = hold.doubleValue()*1e3;

            // total moles of ion i

            this.molesT[i]  = this.concnA[i]*volumeA + this.concnB[i]*volumeB;

            if(this.molesT[i]>0.0D){

                indexC[i] = 1;

            }

            else{

                indexC[i] = 0;

            }

            // association constant

            hold            = (Double)this.arrayl.get(3+i*7);

            this.assocConsts[i] = hold.doubleValue()*1e-3;

            if(this.assocConsts[i]>0.0D){

                indexK[ii] = i;

                ii++;

            }

            // ion radius

            hold            = (Double)this.arrayl.get(4+i*7);

            this.radii[i]   = hold.doubleValue();

            // ion charge

            holi            = (Integer)this.arrayl.get(5+i*7);

            this.charges[i] = holi.intValue();

            // running check for all ions having same absolute charge

            if(i==0){

                this.chargeValue = Math.abs(this.charges[0]);

            }

            else{

                if(Math.abs(this.charges[i])!=this.chargeValue)this.chargeSame=false;

            }

            // partition coefficient

            hold        = (Double)this.arrayl.get(6+i*7);