/*
* Classes GouyChapmanStern
*
* Class GouyChapmanStern contains the methods for
* calculating a surface potential, surface charge,
* potential profiles and ionic profiles at an
* electrolyte - charged surface interface
* using the Gouy-Chapman or Gouy-Chapman-Stern models.
*
* Class GouyChapmanStern requires Class GCSMinim that
* implements an interface to the required minimisation method
*
* WRITTEN BY: Michael Thomas Flanagan
*
* DATE: 1 December 2004
* UPDATE: 26 February 2005, 5-7 July 2008
*
* DOCUMENTATION:
* See Michael Thomas Flanagan's JAVA library on-line web page:
* http://www.ee.ucl.ac.uk/~mflanaga/java/GouyChapmanStern.html
* http://www.ee.ucl.ac.uk/~mflanaga/java/
*
* Copyright (c) December 2004, February 2005
*
* 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.
*
***************************************************************************************/
package flanagan.physchem;
import java.util.ArrayList;
import flanagan.physprop.IonicRadii;
import flanagan.io.*;
import flanagan.math.*;
import flanagan.integration.Integration;
import flanagan.integration.IntegralFunction;
import javax.swing.JOptionPane;
// Class to evaluate the function needed to evaluate the potential at any distance x
// Equation 32a in the documentation, GouyChapmanStern.html
// used in the method, getPotentialAtX(), in the class GouyChapmanStern
class FunctionPatX implements IntegralFunction{
public int numOfIons = 0;
public double termOne = 0.0D;
public double expTerm = 0.0D;
public double[] bulkConcn = null;
public double[] charges = null;
public double function(double x){
double sigma = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigma += (this.bulkConcn[i]*this.termOne)*(Math.exp(-this.expTerm*this.charges[i]*x) - 1.0D);
}
return 1.0D/Math.sqrt(sigma);
}
}
// GouyChapmanStern Class
public class GouyChapmanStern{
private ArrayList<Object> vec = new ArrayList<Object>();
// vector holding:
// element 0: name of first ion
// element 1: initial concentration of the ion in the electrolyte (moles per cubic metre but entered as moles per cubic decimetre)
// element 2: radius of first ion (metres)
// element 3: charge of first ion
// element 4: association constant of ion with surface sites (cubic metres per mole but entered as cubic decimetres per mole)
// element 5: name of second ion
// etc
private boolean unpackArrayList = false; // = true when the above vector has been unpacked for the calculations
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[] initConcnM = null; // initial concentration of ions (Molar)
private double[] initConcn = null; // initial concentration of ions (mol per cubic metre)
private double[] siteConcn = null; // surface concentration of ions adsorbed
private double[] sternConcn = null; // surface concentration of ions in Stern layer but not specifically adsorbed
private double[] bulkConcn = null; // concentration of ions in bulk phase (mol per cubic metre)
private double electrolyteConcn = 0.0D; // electrolyte concentration (average if asymmetric)
private double ionicStrength = 0.0D; // ionic strength of the electrolyte
private int[] indexK = null; // ion index of ionic species with non zero assocConsts
private int nonZeroAssocK = 0; // number of ionic species with affinity for surface sites
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 tolNeutral = 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 double[] assocConstsM = null; // association constants of the ions with surface sites (M^-1)
private double[] assocConsts = null; // association constants of the ions with surface sites (cubic metres per mol)
private double surfaceSiteDensity = 0.0D; // surface site density - entered as number per square metre
// - converted to moles per square metre
private double freeSurfaceSiteDensity = 0.0D; // surface site density minus all occupied site densities
private boolean surfaceDensitySet = false; // = true if the surface density is enterd
private double epsilon = 0.0D; // relative electrical permittivity of electrolyte
private double epsilonStern = 0.0D; // relative electrical permittivity of the Stern layer
private boolean epsilonSet = false; // = true when epsilon has been set
private double temp = 25.0; // Temperature (degrees Celcius) [Enter temperature in degrees Celsius]
private double tempK = 25.0-Fmath.T_ABS; // Temperature (degrees Kelvin)
private boolean tempSet = false; // = true when temperature has been set
private double surfacePotential = 0.0D; // Surface potential relative to the bulk potential [volts]
private boolean psi0set = false; // = true when surface potential has been entered
private double diffPotential = 0.0D; // diffuse layer potential difference [volts]
private double sternPotential = 0.0D; // Stern potential difference [volts]
private double surfaceArea = 1.0D; // surface area in square metres
private boolean surfaceAreaSet = false; // = true if the surface area is enterd
private double volume = 1.0D; // electrolyte volume in cubic metres
private boolean volumeSet = false; // = true if the volume is enterd
private double sternCap = 0.0D; // Stern layer capacitance [F]
private double diffCap = 0.0D; // diffuse double layer capacitance [F]
private double totalCap = 0.0D; // total surface capacitance [F]
private double sternDelta = 0.0D; // Stern layer thickness [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 charge value not the same for all ions
private double averageCharge = 0; // number average absolute charge value of the ions
private double surfaceChargeDensity = 0.0D; // surface charge Density [C per square metre]
private double adsorbedChargeDensity = 0.0D; // adsorbed charge Density [C per square metre]
private double diffuseChargeDensity = 0.0D; // diffuse layer charge Density [C per square metre]
private boolean sigmaSet = false; // = true when surface charge density has been set
private double surfaceCharge = 0.0D; // surface charge [C]
private boolean chargeSet = false; // = true when surface charge set
private double recipKappa = 0.0D; // Debye length
private boolean sternOption = true; // = true: Guoy-Chapman-Stern theory used
// = false: Gouy-Chapman theory used
private double expTerm = 0.0D; // common group of constants: e/kT
private double expTermOver2 = 0.0D; // common group of constants: e/2kT
private double twoTerm = 0.0D; // common group of constants: 2.N.k.T.eps0.eps
private double eightTerm = 0.0D; // common group of constants: 8.N.k.T.eps0.eps
// Constructor
public GouyChapmanStern(){
}
// 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 ignoring the Stern layer
// i.e. using Gouy-Chapman theory NOT Gouy-Chapman-Stern theory
public void ignoreStern(){
this.sternOption = false;
}
// Method for reinstating the Stern layer
// i.e. using Gouy-Chapman-Stern theory NOT Gouy-Chapman-theory
public void includeStern(){
this.sternOption = true;
}
// Method to add an ionic species
// Concentrations - Molar, assocK - M^-1, radius - metres, charge - valency e.g. +1
public void setIon(String ion, double concn, double radius, int charge, double assocK){
this.vec.add(ion);
this.vec.add(new Double(concn));
this.vec.add(new Double(radius));
this.vec.add(new Integer(charge));
this.vec.add(new Double(assocK));
if(assocK!=0.0D)this.nonZeroAssocK++;
this.numOfIons++;
this.unpackArrayList = false;
}
// Method to add an ionic species
// Concentrations - Molar, radius - metres, charge - valency e.g. +1
// association constant = 0.0D
public void setIon(String ion, double concn, double radius, int charge){
this.vec.add(ion);
this.vec.add(new Double(concn));
this.vec.add(new Double(radius));
this.vec.add(new Integer(charge));
this.vec.add(new Double(0.0D));
this.numOfIons++;
this.unpackArrayList = false;
}
// Method to add an ionic species
// 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 concn, double assocK){
IonicRadii ir = new IonicRadii();
this.vec.add(ion);
this.vec.add(new Double(concn));
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.vec.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 as, e.g +2";
charg = Db.readInt(mess1+mess2);
}
this.vec.add(new Integer(charg));
this.vec.add(new Double(assocK));
if(assocK!=0.0D)this.nonZeroAssocK++;
this.numOfIons++;
this.unpackArrayList = false;
}
// Method to add an ionic species
// default radii and charge taken from class IonicRadii
// association constant = 0.0D
// Concentrations - Molar
public void setIon(String ion, double concn){
IonicRadii ir = new IonicRadii();
this.vec.add(ion);
this.vec.add(new Double(concn));
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.vec.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 as, e.g +2";
charg = Db.readInt(mess1+mess2);
}
this.vec.add(new Integer(charg));
this.vec.add(new Double(0.0D));
this.numOfIons++;
this.unpackArrayList = false;
}
// Method to set the surface adsorption site density, moles per square metre
public void setSurfaceSiteDensity(double density){
this.surfaceSiteDensity = density/Fmath.N_AVAGADRO;
this.surfaceDensitySet = true;
}
// Method to set the surface area
public void setSurfaceArea(double area){
this.surfaceArea = area;
this.surfaceAreaSet = true;
}
// Method to set the electrolyte volume
public void setVolume(double vol){
this.volume = vol;
this.volumeSet = true;
}
// Method to set relative permittivities
public void setRelPerm(double epsilon, double epsilonStern){
this.epsilon = epsilon;
this.epsilonStern = epsilonStern;
this.epsilonSet = true;
}
// Method to set relative permittivity
// No Stern layer permittivities included
public void setRelPerm(double epsilon){
this.epsilon = epsilon;
this.epsilonSet = true;
}
// Method to set temperature (enter as degrees Celsius)
public void setTemp(double temp){
this.tempK = temp - Fmath.T_ABS;
this.tempSet = true;
}
// Method to set Surface Charge Density (C per square metre)
public void setSurfaceChargeDensity(double sigma){
if(this.psi0set){
System.out.println("You have already entered a surface potential");
System.out.println("This class allows the calculation of a surface charge density for a given surface potential");
System.out.println("or the calculation of a surface potential for a given surface charge density");
System.out.println("The previously entered surface potential will now be ignored");
this.psi0set = false;
}
this.surfaceChargeDensity = sigma;
this.sigmaSet = true;
if(this.surfaceAreaSet){
this.surfaceCharge = sigma*this.surfaceArea;
this.chargeSet = true;
}
}
// Method to set Surface Charge(C) and surface area
public void setSurfaceCharge(double charge, double area){
if(this.psi0set){
System.out.println("You have already entered a surface potential");
System.out.println("This class allows the calculation of a surface charge density for a given surface potential");
System.out.println("or the calculation of a surface potential for a given surface charge density");
System.out.println("The previously entered surface potential will now be ignored");
this.psi0set = false;
}
this.surfaceCharge = charge;
this.chargeSet = true;
this.surfaceArea = area;
this.surfaceAreaSet = true;
this.surfaceChargeDensity = charge/this.surfaceArea;
this.sigmaSet = true;
}
// Method to set Surface Charge(C)
public void setSurfaceCharge(double charge){
if(this.psi0set){
System.out.println("You have already entered a surface potential");
System.out.println("This class allows the calculation of a surface charge density for a given surface potential");
System.out.println("or the calculation of a surface potential for a given surface charge density");
System.out.println("The previously entered surface potential will now be ignored");
this.psi0set = false;
}
this.surfaceCharge = charge;
this.chargeSet = true;
if(this.surfaceAreaSet){
this.surfaceChargeDensity = charge/this.surfaceArea;
this.sigmaSet = true;
}
}
// Method to set Surface Potential (Volts)
public void setSurfacePotential(double psi0){
if(this.sigmaSet){
System.out.println("You have already entered a surface charge density");
System.out.println("This class allows the calculation of a surface potential for a given surface charge density");
System.out.println("or the calculation of a surface charge density for a given surface potential");
System.out.println("The previously entered surface charge density will now be ignored");
this.sigmaSet = false;
}
this.surfacePotential = psi0;
this.psi0set = true;
}
// Method to get the diffuse double layer potential [volts]
public double getDiffuseLayerPotential(){
if(!this.sternOption){
System.out.println("Class: GouyChapmanStern\nMethod: getDiffuseLayerPotential\nThe Stern modification was not included");
System.out.println("The value of the diffuse layer potential has been set equal to the surface potential");
return getSurfacePotential();
}
if(this.psi0set && this.sigmaSet){
return this.diffPotential;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.diffPotential;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.diffPotential;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getDiffuseLayerPotential\nThe value of the diffuse layer potential has not been calculated\nzero returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return 0.0D;
}
}
}
}
// Method to get the Stern layer potential [volts]
public double getSternLayerPotential(){
if(!this.sternOption){
System.out.println("Class: GouyChapmanStern\nMethod: getSternLayerPotential\nThe Stern modification has not been included");
System.out.println("The value of zero has been returned");
return 0.0D;
}
if(this.psi0set && this.sigmaSet){
return this.sternPotential;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.sternPotential;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.sternPotential;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getSternLayerPotential\nThe value of the Stern layer potential has not been calculated\nzero returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return 0.0D;
}
}
}
}
// Method to get the Stern Capacitance per square metre
public double getSternCapPerSquareMetre(){
if(!this.sternOption){
System.out.println("Class: GouyChapmanStern\nMethod: getSternCapacitance\nThe Stern modification has not been included");
System.out.println("A value of infinity has been returned");
return Double.POSITIVE_INFINITY;
}
if(this.psi0set && this.sigmaSet){
return this.sternCap;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.sternCap;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.sternCap;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getSternCap\nThe value of the Stern capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the Stern Capacitance
public double getSternCapacitance(){
if(!this.sternOption){
System.out.println("Class: GouyChapmanStern\nMethod: getSternCapacitance\nThe Stern modification has not been included");
System.out.println("A value of infinity has been returned");
return Double.POSITIVE_INFINITY;
}
if(!this.surfaceAreaSet){
System.out.println("Class: GouyChapmanStern\nMethod: getSternCapacitance\nThe surface area has not bee included");
System.out.println("A value per square metre has been returned");
return this.sternCap;
}
if(this.psi0set && this.sigmaSet){
return this.sternCap*this.surfaceArea;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.sternCap*this.surfaceArea;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.sternCap*this.surfaceArea;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getSternCapacitance\nThe value of the Stern capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the diffuse layer Capacitance per square metre
public double getDiffuseLayerCapPerSquareMetre(){
if(this.psi0set && this.sigmaSet){
return this.diffCap;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.diffCap;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.diffCap;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getDiffuseLayerCapPerSquareMetre\nThe value of the diffuse layer capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the diffuse layer Capacitance
public double getDiffuseLayerCapacitance(){
if(!this.surfaceAreaSet){
System.out.println("Class: GouyChapmanStern\nMethod: getDiffuseLayerCapacitance\nThe surface area has not bee included");
System.out.println("A value per square metre has been returned");
return this.diffCap;
}
if(this.psi0set && this.sigmaSet){
return this.diffCap*this.surfaceArea;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.diffCap*this.surfaceArea;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.diffCap*this.surfaceArea;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getDiffuseLayerCap\nThe value of the diffuse layer capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the total capacitance per square metre
public double getTotalCapPerSquareMetre(){
if(this.psi0set && this.sigmaSet){
return this.totalCap;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.totalCap;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.totalCap;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getTotalCapPerSquareMetre\nThe value of the total capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the total capacitance
public double getTotalCapacitance(){
if(!this.surfaceAreaSet){
System.out.println("Class: GouyChapmanStern\nMethod: getTotalCapacitance\nThe surface area has not bee included");
System.out.println("A value per square metre has been returned");
return this.diffCap;
}
if(this.psi0set && this.sigmaSet){
return this.totalCap*this.surfaceArea;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.totalCap*this.surfaceArea;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.totalCap*this.surfaceArea;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getTotalCapacitance\nThe value of the total capacitance has not been calculated\ninfinity returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return Double.POSITIVE_INFINITY;
}
}
}
}
// Method to get the Stern layer thickness
public double getSternThickness(){
if(!this.sternOption){
System.out.println("Class: GouyChapmanStern\nMethod: getSternThickness");
System.out.println("The Stern modification has not been included");
System.out.println("A value of zero has been returned");
return 0.0D;
}
if(this.psi0set && this.sigmaSet){
return this.sternDelta;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.sternDelta;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.sternDelta;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getSternThickness\nThe value of the Stern thickness has not been calculated\nzero returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return 0.0D;
}
}
}
}
// Method to get the Debye length
public double getDebyeLength(){
if(!this.unpackArrayList)unpack();
return this.calcDebyeLength();
}
// Calculate Debye length
private double calcDebyeLength(){
if(!this.epsilonSet)throw new IllegalArgumentException("The relative permittivitie/s have not been entered");
if(!this.tempSet)System.out.println("The temperature has not been entered\na value of 25 degrees Celsius has been used");
double preterm = 2.0D*Fmath.square(Fmath.Q_ELECTRON)*Fmath.N_AVAGADRO/(Fmath.EPSILON_0*this.epsilon*Fmath.K_BOLTZMANN*this.tempK);
this.recipKappa = 0.0;
for(int i=0; i<this.numOfIons; i++){
this.recipKappa += this.bulkConcn[i]*charges[i]*charges[i];
}
this.recipKappa = 1.0D/Math.sqrt(this.recipKappa*preterm);
return this.recipKappa;
}
// Method to get the ionic strength
public double getIonicStrength(){
if(!this.unpackArrayList)unpack();
return this.ionicStrength;
}
// unpacks the ion storage vector and fills relevant arrays
private void unpack(){
if(this.numOfIons==0)throw new IllegalArgumentException("No ions have been entered");
// 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.siteConcn = new double[this.numOfIons];
this.sternConcn = new double[this.numOfIons];
this.initConcnM = new double[this.numOfIons];
this.initConcn = new double[this.numOfIons];
this.bulkConcn = new double[this.numOfIons];
this.radii = new double[this.numOfIons];
this.charges = new double[this.numOfIons];
this.assocConsts = new double[this.numOfIons];
this.assocConstsM = new double[this.numOfIons];
this.indexK = new int[this.nonZeroAssocK];
Double hold = null;
Integer holi = null;
int ii = 0;
this.chargeValue = 0;
this.chargeSame = true;
for(int i=0; i<numOfIons; i++){
this.ionNames[i] = (String)this.vec.get(0+i*5);
hold = (Double)this.vec.get(1+i*5);
this.initConcnM[i] = hold.doubleValue();
this.initConcn[i] = this.initConcnM[i]*1e3;
hold = (Double)this.vec.get(2+i*5);
this.radii[i] = hold.doubleValue();
holi = (Integer)this.vec.get(3+i*5);
this.charges[i] = holi.intValue();
hold = (Double)this.vec.get(4+i*5);
this.assocConstsM[i]= hold.doubleValue();
this.assocConsts[i] = this.assocConstsM[i]*1e-3;
if(this.assocConsts[i]>0.0D){
indexK[ii] = i;
ii++;
}
// 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;
}
// calculate number of anionic and cationic species
if(charges[i]<0){
numOfAnions++;
}
else{
numOfCations++;
}
if(this.surfaceSiteDensity==0.0D)this.nonZeroAssocK = 0;
}
// Calculate overall charge, average absolute charge,
// ionic strength and electrolye concentration (only average if asymmetric)
this.averageCharge = 0.0D;
this.ionicStrength = 0.0D;
double overallCharge = 0.0D;
double positives = 0.0D;
double negatives = 0.0D;
for(int i=0; i<numOfIons; i++){
if(charges[i]>0.0D){
positives += this.initConcn[i]*charges[i];
}
else{
negatives += this.initConcn[i]*charges[i];
}
overallCharge = positives + negatives;
}
if(Math.abs(overallCharge)>positives*this.tolNeutral){
String quest0 = "Class: GouyChapmanStern, method: unpack()\n";
String quest1 = "Total charge = " + overallCharge + " mol/dm, i.e. is not equal to zero\n";
String quest2 = "Positive charge = " + positives + " mol/dm\n";
String quest3 = "Do you wish to continue?";
String quest = quest0 + quest1 + quest2 + quest3;
int res = JOptionPane.showConfirmDialog(null, quest, "Neutrality check", JOptionPane.YES_NO_OPTION, JOptionPane.QUESTION_MESSAGE);
if(res==1)System.exit(0);
}
double numer = 0.0D;
double denom = 0.0D;
double anionConc = 0.0D;
double cationConc = 0.0D;
for(int i=0; i<numOfIons; i++){
this.ionicStrength += initConcn[i]*charges[i]*charges[i];
if(charges[i]>0){
cationConc += initConcn[i];
}
else{
anionConc += initConcn[i];
}
if(initConcn[i]>0.0D){
numer += initConcn[i]*Math.abs(charges[i]);
denom += initConcn[i];
}
}
this.ionicStrength = this.ionicStrength*1e-3/2.0D;
this.averageCharge = numer/denom;
this.electrolyteConcn = (anionConc + cationConc)/2.0D;
// initialise site and Stern concentrations
for(int i=0; i<this.numOfIons; i++){
bulkConcn[i] = initConcn[i];
siteConcn[i] = 0.0D;
sternConcn[i] = 0.0D;
}
// calculate commonly used combinations of constants
this.expTerm = -Fmath.Q_ELECTRON/(Fmath.K_BOLTZMANN*this.tempK); // |e|/kT
this.expTermOver2 = this.expTerm/2.0D; // |e|/2kT
this.twoTerm = 2.0D*(Fmath.N_AVAGADRO*Fmath.K_BOLTZMANN)*Fmath.EPSILON_0*this.epsilon*this.tempK; // 2.N.k.T.eps0.eps
this.eightTerm = 4.0D*this.twoTerm; // 8.N.k.T.eps0.eps
this.unpackArrayList = true;
}
// Calculate the surface charge density for the set surface potential
public double getSurfaceChargeDensity(){
if(this.sigmaSet && this.psi0set){
return this.surfaceChargeDensity;
}
else{
if(!this.psi0set)throw new IllegalArgumentException("No surface potential has been entered");
return getSurfaceChargeDensity(this.surfacePotential);
}
}
// Get the surface charge density for a given surface potential, psi
public double getSurfaceChargeDensity(double psi){
this.surfaceChargeDensity = this.calcSurfaceChargeDensity(psi);
this.sigmaSet = true;
if(this.surfaceAreaSet){
this.surfaceCharge = this.surfaceChargeDensity*this.surfaceArea;
this.chargeSet = true;
}
return this.surfaceChargeDensity;
}
// Calculate the surface charge density for a given surface potential, psi
private double calcSurfaceChargeDensity(double psi){
double surCharDen = 0.0D; // temporary values of surface charge density
if(!this.epsilonSet)throw new IllegalArgumentException("The relative permittivitie/s have not been entered");
if(!this.tempSet)System.out.println("The temperature has not been entered\na value of 25 degrees Celsius has been used");
if(!this.unpackArrayList)unpack();
if(this.sternOption){
if(!this.surfaceAreaSet)throw new IllegalArgumentException("The surface area has not been entered");
if(!this.volumeSet)throw new IllegalArgumentException("The electrolyte volume has not been entered");
// Gouy-Chapman-Stern
if(this.nonZeroAssocK==0){
// No specific adsorption
if(this.chargeSame){
// symmetric electrolyte
surCharDen = surfaceChargeDensity0(psi);
}
else{
// asymmetric electrolyte
surCharDen = surfaceChargeDensity1(psi);
}
}
else{
// Specific Adsorption
if(this.chargeSame){
// symmetric electrolyte
surCharDen = surfaceChargeDensity2(psi);
}
else{
// asymmetric electrolyte
surCharDen = surfaceChargeDensity3(psi);
}
}
}
else{
// Gouy-Chapman
if(this.chargeSame){
// symmetric electrolyte
surCharDen = Math.sqrt(this.eightTerm*this.electrolyteConcn)*Fmath.sinh(this.chargeValue*this.expTermOver2*psi);
}
else{
// asymmetric electrolyte
double sigmaSum = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaSum += bulkConcn[i]*(Math.exp(-this.expTerm*psi*charges[i])-1.0D);
}
surCharDen = Fmath.sign(psi)*Math.sqrt(this.twoTerm*sigmaSum);
}
}
this.totalCap = surCharDen/psi;
if(!this.sternOption){
this.diffPotential = psi;
this.sternCap = Double.POSITIVE_INFINITY;
this.sternPotential = 0.0D;
this.diffCap = this.totalCap;
}
else{
this.diffPotential = psi - surCharDen/this.sternCap;
this.sternPotential = psi - this.diffPotential;
this.diffCap = (surCharDen + this.adsorbedChargeDensity)/this.diffPotential;
}
return surCharDen;
}
// Calculate the surface charge for the set surface potential
public double getSurfaceCharge(){
return getSurfaceCharge(this.surfacePotential);
}
// Calculate the surface charge for a given surface potential, psi
public double getSurfaceCharge(double psi){
if(!this.surfaceAreaSet)throw new IllegalArgumentException("No surface area has been entered");
if(this.sigmaSet){
this.surfaceCharge = this.surfaceChargeDensity*this.surfaceArea;
}
else{
if(!this.psi0set)throw new IllegalArgumentException("No surface potential has been entered");
this.surfaceCharge = this.getSurfaceChargeDensity(psi)*this.surfaceArea;
}
return this.surfaceCharge;
}
// Calculate surface charge density for a symmetric electrolyte
// Stern modification without specific adsorption
private double surfaceChargeDensity0(double psi){
double surCharDen = 0.0D; // temporary values of surface charge density
// bisection method
double ionSum = 0.0D;
for(int i=0; i<this.numOfIons; i++){
if(this.charges[i]>0.0){
ionSum += bulkConcn[i];
}
}
double sigmaLow = 0.0D;
double sFuncLow = this.sigmaFunction0(sigmaLow, psi);
double sigmaHigh = Math.sqrt(this.eightTerm*ionSum)*Fmath.sinh(this.chargeValue*this.expTermOver2*psi);
double sFuncHigh = this.sigmaFunction0(sigmaHigh, psi);
if(sFuncHigh*sFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(sigmaHigh)*1e-6;
boolean test = true;
double sigmaMid = 0.0D;
double sFuncMid = 0.0D;
int nIter = 0;
while(test){
sigmaMid = (sigmaLow + sigmaHigh)/2.0D;
sFuncMid = this.sigmaFunction0(sigmaMid, psi);
if(Math.abs(sFuncMid)<=check){
surCharDen = sigmaMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfaceChargeDensity0\nnumber of iterations exceeded in bisection\ncurrent value of sigma returned");
surCharDen = sigmaMid;
test = false;
}
else{
if(sFuncMid*sFuncLow>0){
sigmaLow=sigmaMid;
sFuncLow=sFuncMid;
}
else{
sigmaHigh=sigmaMid;
sFuncHigh=sFuncMid;
}
}
}
}
return surCharDen;
}
// function to calculate sigma for surfaceChargeDensity0()
private double sigmaFunction0(double sigma, double psi){
// calculate psi(delta) Stern capacitance for current estimate of sigma
this.calcSurfacePotential(sigma);
return this.diffPotential - psi + sigma/this.sternCap;
}
// Calculate surface charge density for a asymmetric electrolyte
// Stern modification without specific adsorption
private double surfaceChargeDensity1(double psi){
double surCharDen = 0.0D; // temporary values of surface charge density
// bisection method
double sigmaLow = 0.0D;
double sFuncLow = this.sigmaFunction0(sigmaLow, psi);
double sigmaHigh = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaHigh += bulkConcn[i]*this.twoTerm*(Math.exp(-this.expTerm*charges[i]*psi) - 1.0D);
}
sigmaHigh = Fmath.sign(psi)*Math.sqrt(this.twoTerm*sigmaHigh);
double sFuncHigh = this.sigmaFunction0(sigmaHigh, psi);
if(sFuncHigh*sFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(sigmaHigh)*1e-6;
boolean test = true;
double sigmaMid = 0.0D;
double sFuncMid = 0.0D;
int nIter = 0;
while(test){
sigmaMid = (sigmaLow + sigmaHigh)/2.0D;
sFuncMid = this.sigmaFunction0(sigmaMid, psi);
if(Math.abs(sFuncMid)<=check){
surCharDen = sigmaMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfaceChargeDensity1\nnumber of iterations exceeded in outer bisection\ncurrent value of sigma returned");
surCharDen = sigmaMid;
test = false;
}
else{
if(sFuncLow*sFuncMid>0.0D){
sigmaLow = sigmaMid;
sFuncLow = sFuncMid;
}
else{
sigmaHigh = sigmaMid;
sFuncHigh = sFuncMid;
}
}
}
}
return surCharDen;
}
// Calculate the Stern thickness, adsorbed ion concentrations and Stern layer non-adsorbed ion concentrations
// for a given potential, psi
private double calcDelta(double psi){
double numer = 0.0D;
double denom = 0.0D;
double delta = 0.0;
for(int i=0; i<this.numOfIons; i++){
this.sternConcn[i] = this.bulkConcn[i]*Math.exp(-this.charges[i]*this.expTerm);
numer += this.sternConcn[i]*this.radii[i];
denom += this.sternConcn[i];
}
this.sternDelta = numer/denom;
return this.sternDelta;
}
// Calculate surface charge density for a symmetric electrolyte
// Stern modification with specific adsorption
private double surfaceChargeDensity2(double psi){
double surCharDen = 0.0D; // temporary values of surface charge density
// bisection method
double sigmaAdsPos = 0.0D;
double sigmaAdsNeg = 0.0D;
for(int i=0; i<this.nonZeroAssocK; i++){
if(this.charges[indexK[i]]>0){
sigmaAdsPos=surfaceSiteDensity;
}
else{
sigmaAdsNeg=-surfaceSiteDensity;
}
}
double sigmaLow = 0.0D;
double sigmaHigh = 0.0D;
double ionSum = 0.0D;
for(int i=0; i<this.numOfIons; i++){
ionSum += bulkConcn[i];
}
ionSum /=2.0D;
sigmaHigh = Math.sqrt(ionSum*this.eightTerm)*Fmath.sinh(this.expTermOver2*psi*this.chargeValue);
if(sigmaHigh>0.0D){
sigmaHigh += sigmaAdsPos;
sigmaLow -= sigmaAdsNeg;
}
else{
sigmaHigh -= sigmaAdsNeg;
sigmaLow += sigmaAdsPos;
}
System.out.println("pp0 "+ sigmaLow + " " + psi + " "+ionSum + " " +ionSum*this.eightTerm);
double sFuncLow = this.sigmaFunction2(sigmaLow, psi);
System.out.println("pp1 "+ sigmaHigh + " " + psi + " "+ionSum + " " +this.chargeValue+" "+Math.sqrt(ionSum*this.eightTerm) + " " + this.expTermOver2*psi*this.chargeValue + " " + Fmath.sinh(this.expTermOver2*psi*this.chargeValue));
double sFuncHigh = this.sigmaFunction2(sigmaHigh, psi);
if(sFuncHigh*sFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(sigmaHigh)*1e-6;
boolean test = true;
double sigmaMid = 0.0D;
double sFuncMid = 0.0D;
int nIter = 0;
while(test){
sigmaMid = (sigmaLow + sigmaHigh)/2.0D;
sFuncMid = this.sigmaFunction2(sigmaMid, psi);
if(Math.abs(sFuncMid)<=check){
surCharDen = sigmaMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfaceChargeDensity2\nnumber of iterations exceeded in outer bisection\ncurrent value of sigma returned");
surCharDen = sigmaMid;
test = false;
}
else{
if(sFuncLow*sFuncMid>0.0D){
sigmaLow = sigmaMid;
sFuncLow = sFuncMid;
}
else{
sigmaHigh = sigmaMid;
sFuncHigh = sFuncMid;
}
}
}
}
return surCharDen;
}
// function to calculate sigma for surfaceChargeDensity2()
private double sigmaFunction2(double sigma, double psi){
System.out.println("F2begin");
// bisection method for psi(delta)
double psiLow = -10*psi;
double pFuncLow = this.psiFunctionQ(psiLow, psi, sigma);
double psiHigh = 10*psi;
double pFuncHigh = this.psiFunctionQ(psiHigh, psi, sigma);
System.out.println("qq " + pFuncLow + " "+pFuncHigh+" "+ psiLow+" "+psiHigh);
if(pFuncHigh*pFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(psi)*1e-6;
boolean test = true;
double psiMid = 0.0D;
double pFuncMid = 0.0D;
int nIter = 0;
while(test){
psiMid = (psiLow + psiHigh)/2.0D;
pFuncMid = this.psiFunctionQ(psiMid, psi, sigma);
if(Math.abs(pFuncMid)<=check){
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfaceChargeDensity3\nnumber of iterations exceeded in inner bisection\ncurrent value of sigma returned");
test = false;
}
if(pFuncMid*pFuncHigh>0){
psiHigh = psiMid;
pFuncHigh = pFuncMid;
}
else{
psiLow = psiMid;
pFuncLow = pFuncMid;
}
}
}
double sigmaEst = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaEst += bulkConcn[i];
}
sigmaEst = Math.sqrt(this.eightTerm*sigmaEst/2.0D)*Fmath.sinh(this.expTermOver2*psi*this.chargeValue);
System.out.println("F2end");
return sigma + this.adsorbedChargeDensity - sigmaEst;
}
// Calculate surface charge density for a asymmetric electrolyte
// Stern modification with specific adsorption
private double surfaceChargeDensity3(double psi){
double surCharDen = 0.0D; // temporary values of surface charge density
// bisection method
double sigmaAdsPos = 0.0D;
double sigmaAdsNeg = 0.0D;
for(int i=0; i<this.nonZeroAssocK; i++){
if(this.charges[indexK[i]]>0){
sigmaAdsPos=surfaceSiteDensity;
}
else{
sigmaAdsNeg=-surfaceSiteDensity;
}
}
double sigmaLow = 0.0D;
double sigmaHigh = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaHigh += bulkConcn[i]*this.twoTerm*(Math.exp(-this.expTerm*charges[i]*psi) - 1.0D);
}
sigmaHigh = Fmath.sign(psi)*Math.sqrt(sigmaHigh);
if(sigmaHigh>0.0D){
sigmaHigh += sigmaAdsPos;
sigmaLow -= sigmaAdsNeg;
}
else{
sigmaHigh -= sigmaAdsNeg;
sigmaLow += sigmaAdsPos;
}
double sFuncLow = this.sigmaFunction3(sigmaLow, psi);
double sFuncHigh = this.sigmaFunction3(sigmaHigh, psi);
if(sFuncHigh*sFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(sigmaHigh)*1e-6;
boolean test = true;
double sigmaMid = 0.0D;
double sFuncMid = 0.0D;
int nIter = 0;
while(test){
sigmaMid = (sigmaLow + sigmaHigh)/2.0D;
sFuncMid = this.sigmaFunction3(sigmaMid, psi);
if(Math.abs(sFuncMid)<=check){
surCharDen = sigmaMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfaceChargeDensity3\nnumber of iterations exceeded in outer bisection\ncurrent value of sigma returned");
surCharDen = sigmaMid;
test = false;
}
else{
if(sFuncLow*sFuncMid>0.0D){
sigmaLow = sigmaMid;
sFuncLow = sFuncMid;
}
else{
sigmaHigh = sigmaMid;
sFuncHigh = sFuncMid;
}
}
}
}
return surCharDen;
}
// function to calculate sigma for surfaceChargeDensity2()
private double sigmaFunction3(double sigma, double psi){
// bisection method for psi(delta)
double psiLow = 0.0D;
double pFuncLow = this.psiFunctionQ(psiLow, psi, sigma);
double psiHigh = psi;
double pFuncHigh = this.psiFunctionQ(psiHigh, psi, sigma);
if(pFuncHigh*pFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(psi)*1e-6;
boolean test = true;
double psiMid = 0.0D;
double pFuncMid = 0.0D;
int nIter = 0;
while(test){
psiMid = (psiLow + psiHigh)/2.0D;
pFuncMid = this.psiFunctionQ(psiMid, psi, sigma);
if(Math.abs(pFuncMid)<=check){
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: sigmaFunction3\nnumber of iterations exceeded in inner bisection\ncurrent value of sigma returned");
test = false;
}
if(pFuncMid*pFuncHigh>0){
psiHigh = psiMid;
pFuncHigh = pFuncMid;
}
else{
psiLow = psiMid;
pFuncLow = pFuncMid;
}
}
}
double sigmaEst = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaEst += bulkConcn[i]*this.twoTerm*(Math.exp(-this.expTerm*charges[i]*psiMid) - 1.0D);
}
sigmaEst = Fmath.sign(psiMid)*Math.sqrt(sigmaEst);
return sigma + this.adsorbedChargeDensity - sigmaEst;
}
// Function to calculate psi(delta) when specific adsorption occurs
private double psiFunctionQ(double psiDelta, double psi0, double sigma){
this.sternDelta = this.calcDeltaQ(psiDelta);
this.diffPotential = psiDelta;
this.sternCap = this.epsilonStern*Fmath.EPSILON_0/this.sternDelta;
System.out.println("aaa " + psiDelta + " " + psi0 + " " + sigma/this.sternCap + " " + sigma);
return psiDelta - psi0 + sigma/this.sternCap;
}
// Calculate the Stern thickness, adsorbed ion concentrations and Stern layer non-adsorbed ion concentrations
// for a given potential, psi
// and speicific adsorption occurring
private double calcDeltaQ(double psi){
double numer = 0.0D;
double denom = 0.0D;
double convFac = this.surfaceArea/this.volume;
int ii = 0;
// calculate adsorbed ion concentrations
if(this.nonZeroAssocK==1){
ii = indexK[0];
double hold = this.assocConsts[ii]*Math.exp(-charges[ii]*psi*this.expTerm);
double aa = hold*convFac;
double bb = -(1.0D + this.initConcn[ii]*hold + this.surfaceSiteDensity*hold*convFac);
double cc = this.initConcn[ii]*this.surfaceSiteDensity*hold;
// solve quadratic equatiom
double root = bb*bb - 4.0D*aa*cc;
if(root<0.0D){
System.out.println("Class: GouyChapmanStern\nMethod: calcDeltaQ\nthe square root term (b2-4ac) of the quadratic = "+root);
System.out.println("this term was set to zero as the negative value MAY have arisen from rounding errors");
root = 0.0D;
}
double qq = -0.5*(bb + Fmath.sign(bb)*Math.sqrt(root));
double root1 = qq/aa;
double root2 = cc/qq;
double limit = this.surfaceSiteDensity*1.001D;
if(root1>=0.0D && root1<=limit){
if(root2<0.0D || root2>limit){
this.siteConcn[indexK[0]] = root1;
this.bulkConcn[indexK[0]] = this.initConcn[indexK[0]] - this.siteConcn[indexK[0]]*this.surfaceArea/this.volume;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: ionConcns");
System.out.println("error1: no physically meaningfull root");
System.out.println("root1 = " + root1 + " root2 = " + root2 + " limit = " + limit);
System.exit(0);
}
}
else{
if(root2>=0.0D && root2<=limit){
if(root1<0.0D || root1>limit){
this.siteConcn[indexK[0]] = root2;
this.bulkConcn[indexK[0]] = this.initConcn[indexK[0]] - this.siteConcn[indexK[0]]*this.surfaceArea/this.volume;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: ionConcns");
System.out.println("error2: no physically meaningfull root");
System.out.println("root1 = " + root1 + " root2 = " + root2 + " limit = " + limit);
System.exit(0);
}
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: ionConcns");
System.out.println("error3: no physically meaningfull root");
System.out.println("root1 = " + root1 + " root2 = " + root2 + " limit = " + limit);
System.exit(0);
}
}
}
else{
// More than one non-zero association constant
// Minimisation procedure
// Initial estimates
double[] vec = new double[this.nonZeroAssocK];
double[][] mat = new double[this.nonZeroAssocK][this.nonZeroAssocK];
// set up simultaneous equations
double expPsiTerm = -psi*this.expTerm;
for(int i=0; i<this.nonZeroAssocK; i++){
ii = indexK[i];
vec[i] = this.assocConsts[ii]*this.initConcn[ii]*this.surfaceSiteDensity*Math.exp(charges[ii]*expPsiTerm);
for(int j=0; j<this.nonZeroAssocK; j++){
mat[i][j] = this.assocConsts[ii]*this.initConcn[ii]*Math.exp(charges[ii]*expPsiTerm);
if(i==j)mat[i][j] += 1.0D;
}
}
// solve simultaneous equations to obtain complex concentrations
Matrix matrix = new Matrix(mat);
vec = matrix.solveLinearSet(vec);
for(int i=0; i<this.nonZeroAssocK; i++){
this.siteConcn[indexK[i]] = vec[i];
}
}
// Create an instances of Minimisation and the minimisation function, GCSminim
Minimisation min = new Minimisation();
GCSminim functA = new GCSminim();
// set minimisation function variables
functA.psiDelta = psi; // current value of psi(delta) [diffuseLayerPotential]
functA.tempK = this.tempK; // temperature
functA.surfaceSiteDensity = this.surfaceSiteDensity; // surface adsorption site density
functA.surfaceArea = this.surfaceArea; // surface ares
functA.volume = this.volume; // electrolyte volume
functA.nonZeroAssocK = this.nonZeroAssocK; // number of non-zero association constants
functA.assocK = this.assocConsts; // surface site association constants
functA.initConcn = this.initConcn; // initial ion concentrations
functA.charges = this.charges; // ion charges
functA.indexK = this.indexK; // indices of the non-zero associaton constants
// set initial estimates, step sizes
double[] start = new double[this.nonZeroAssocK];
double[] step = new double[this.nonZeroAssocK];
for(int i=0; i<this.nonZeroAssocK; i++){
start[i] = this.surfaceSiteDensity/this.nonZeroAssocK;
step[i] = start[i]*0.05D;
}
// covergence tolerance and maximum number of iteration allowed
double tolerance = this.surfaceSiteDensity*1e-8;
int maxIter = 100000;
// Perform minimisation
min.nelderMead(functA, start, step, tolerance, maxIter);
// Get best estimates of the site concentrations
double[] param = min.getParamValues();
for(int i=0; i<this.nonZeroAssocK; i++){
ii = indexK[i];
this.siteConcn[ii] = param[i];
this.bulkConcn[ii] = this.initConcn[ii] - param[i]*this.surfaceArea/this.volume;
}
// Calculte Stern delta and adsorbed charge density
this.adsorbedChargeDensity = 0.0D;
double factor1 = -Fmath.Q_ELECTRON*Fmath.N_AVAGADRO;
double factor2 = this.surfaceArea/this.volume;
for(int i=0; i<this.numOfIons; i++){
this.sternConcn[i] = this.bulkConcn[i]*Math.exp(-this.charges[i]*this.expTerm);
this.adsorbedChargeDensity += this.siteConcn[i]*charges[i]*factor1;
numer += (this.sternConcn[i] + this.siteConcn[i]*factor2)*this.radii[i];
denom += this.sternConcn[i] + this.siteConcn[i]*factor2;
}
double delta = numer/denom;
return delta;
}
// Method to get the adsorbed ion charge density [C per square metre]
public double getAdsorbedChargeDensity(){
if(!this.sternOption || this.nonZeroAssocK==0){
return 0.0D;
}
if(this.psi0set && this.sigmaSet){
return this.adsorbedChargeDensity;
}
else{
if(this.sigmaSet){
this.getSurfacePotential();
return this.sternPotential;
}
else{
if(this.psi0set){
this.getSurfaceChargeDensity();
return this.adsorbedChargeDensity;
}
else{
System.out.println("Class: GouyChapmanStern\nMethod: getAdsorbedChargeDensity\nThe value of the adsorbed ion charge density has not been calculated\nzero returned");
System.out.println("Neither a surface potential nor a surface charge density have been entered");
return 0.0D;
}
}
}
}
// Method to get the diffuse layer charge density [C per square metre]
public double getDiffuseChargeDensity(){
double ads = this.getAdsorbedChargeDensity();
this.diffuseChargeDensity = -(this.surfaceChargeDensity + ads);
return this.diffuseChargeDensity;
}
// Get the surface potential for a given surface charge density, sigma
public double getSurfacePotential(double sigma){
this.surfacePotential = calcSurfacePotential(sigma);
this.psi0set = true;
return this.surfacePotential;
}
// Calculate the surface potential for a given surface charge density, sigma
private double calcSurfacePotential(double sigma){
double surPot = 0.0D; // temporary values of the surface potential
if(!this.epsilonSet)throw new IllegalArgumentException("The relative permittivitie/s have not been entered");
if(!this.tempSet)System.out.println("The temperature has not been entered\na value of 25 degrees Celsius has been used");
if(this.psi0set && this.sigmaSet)return this.surfacePotential;
if(!this.unpackArrayList)unpack();
if(this.sternOption){
// Gouy-Chapman-Stern
if(this.nonZeroAssocK==0){
// No specific adsorption
if(this.chargeSame){
// symmetric electrolyte
this.diffPotential = Fmath.asinh(sigma/Math.sqrt(this.eightTerm*this.electrolyteConcn))/(this.chargeValue*this.expTermOver2);
}
else{
// asymmetric electrolyte
this.diffPotential = surfacePotential1(sigma);
}
this.sternCap = Fmath.EPSILON_0*this.epsilonStern/this.calcDelta(this.diffPotential);
surPot = this.diffPotential + sigma/this.sternCap;
this.totalCap = sigma/this.surfacePotential;
this.diffCap = sigma/this.diffPotential;
}
else{
// Specific absorption
if(this.chargeSame){
// symmetric electrolyte
surPot = surfacePotential2(sigma);
}
else{
// asymmetric electrolyte
surPot = surfacePotential3(sigma);
}
}
}
else{
// Gouy-Chapman
if(this.chargeSame){
// symmetric electrolyte
double preterm = Math.sqrt(this.eightTerm*this.electrolyteConcn);
this.surfacePotential = Fmath.asinh(this.surfaceChargeDensity/preterm)/(this.chargeValue*this.expTermOver2);
}
else{
// asymmetric electrolyte
surPot = surfacePotential4(sigma);
}
this.diffPotential = surPot;
this.sternPotential = 0.0D;
this.totalCap = sigma/surPot;
this.diffCap = sigma/this.diffPotential;
this.sternCap = Double.POSITIVE_INFINITY;
}
return surPot;
}
// Calculate the surface charge for the set surface potential
public double getSurfacePotential(){
if(!this.sigmaSet)throw new IllegalArgumentException("No surface charge density has been entered");
return getSurfacePotential(this.surfaceChargeDensity);
}
// Calculate surface potential for a asymmetric electrolyte
// Stern modification ignored
private double surfacePotential4(double sigma){
double surPot = 0.0D; // temporary values of the surface potential
// bisection method
double psiLow = 0.0D;
double pFuncLow = this.psiFunction4(psiLow, sigma);
double asinhDenom = Math.sqrt(this.eightTerm*this.electrolyteConcn);
double psiHigh = (10.0D/(this.averageCharge*this.expTerm))*Fmath.asinh(sigma/asinhDenom);
double pFuncHigh = this.psiFunction4(psiHigh, sigma);
if(pFuncHigh*pFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(psiHigh)*1e-6;
boolean test = true;
double psiMid = 0.0D;
double pFuncMid = 0.0D;
int nIter = 0;
while(test){
psiMid = (psiLow + psiHigh)/2.0D;
pFuncMid = this.psiFunction4(psiMid, sigma);
if(Math.abs(pFuncMid)<=check){
surPot = psiMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: getSurfacePotential\nnumber of iterations exceeded in outer bisection\ncurrent value of sigma returned");
surPot = psiMid;
test = false;
}
else{
if(pFuncLow*pFuncMid>0.0D){
psiLow = psiMid;
pFuncLow = pFuncMid;
}
else{
psiHigh = psiMid;
pFuncHigh = pFuncMid;
}
}
}
}
return surPot;
}
// function to calculate surfacePotential function for surfacePotential4()
private double psiFunction4(double psi, double sigma){
double sigmaEst = 0.0D;
for(int i=0; i<this.numOfIons; i++){
sigmaEst += bulkConcn[i]*this.twoTerm*(Math.exp(-this.expTerm*charges[i]*psi) - 1.0D);
}
sigmaEst = Fmath.sign(sigma)*Math.sqrt(sigmaEst);
return sigma - sigmaEst;
}
// Calculate surface potential for a given charge density for a asymmetric electrolyte
// Stern modification without specific adsorption
private double surfacePotential1(double sigma){
double difPot = 0.0D; // temporary values of the diffuse potential
// bisection method
double psiLow = 0.0D;
double pFuncLow = this.psiFunction1(psiLow, sigma);
double psiHigh = (5.0D/(this.expTerm*this.chargeValue))*Fmath.asinh(sigma/Math.sqrt(this.eightTerm*this.electrolyteConcn));
double pFuncHigh = this.psiFunction1(psiHigh, sigma);
if(pFuncHigh*pFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(psiHigh)*1e-6;
boolean test = true;
double psiMid = 0.0D;
double pFuncMid = 0.0D;
int nIter = 0;
while(test){
psiMid = (psiLow + psiHigh)/2.0D;
pFuncMid = this.psiFunction1(psiMid, sigma);
if(Math.abs(pFuncMid)<=check){
this.diffPotential = psiMid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: getSurfacePotential\nnumber of iterations exceeded in outer bisection\ncurrent value of sigma returned");
this.diffPotential = psiMid;
test = false;
}
else{
if(pFuncLow*pFuncMid>0.0D){
psiLow = psiMid;
pFuncLow = pFuncMid;
}
else{
psiHigh = psiMid;
pFuncHigh = pFuncMid;
}
}
}
}
return difPot;
}
// function to calculate diffPotential for surfacePotential1()
private double psiFunction1(double psi, double sigma){
double func = 0.0D;
for(int i=0; i<this.numOfIons; i++){
func += this.twoTerm*bulkConcn[i]*(Math.exp(-charges[i]*psi*this.expTerm) - 1.0D);
}
return sigma - Fmath.sign(sigma)*Math.sqrt(func);
}
// method to calculate the potential at a distance, x, from the surface
public double getPotentialAtX(double xDistance){
if(!this.psi0set && !this.sigmaSet)throw new IllegalArgumentException("Neither a surface potential nor a surface charge/density have been entered");
if(this.sigmaSet && !this.psi0set)this.getSurfacePotential();
if(this.psi0set && !this.sigmaSet)this.getSurfaceChargeDensity();
double potAtX = 0.0D;
if(xDistance==0.0D){
// distance, x, is zero metres from charged surface
potAtX = this.surfacePotential;
}
else{
if(this.sternOption){
// distance, x, coincides with the Stern interface
if(xDistance==this.sternDelta){
potAtX = this.diffPotential;
}
else{
if(xDistance<this.sternDelta){
// distance, x, is within the Stern layer
potAtX = this.surfacePotential - (xDistance/this.sternDelta)*(this.surfacePotential - this.diffPotential);
}
else{
// distance, x, is greater than the Stern layer thickness
potAtX = this.calcPotAtX(this.diffPotential, xDistance);
}
}
}
else{
// no Stern layer and distance, x, is greater than zero
potAtX = this.calcPotAtX(this.surfacePotential, xDistance);
}
}
return potAtX;
}
// calculates potential at x for getPotentialAtX()
private double calcPotAtX(double psiL, double xDistance){
double potAtX = 0.0D;
if(this.chargeSame){
//symmetric elctrolyte
double kappa = Math.sqrt(2.0D*Fmath.square(Fmath.Q_ELECTRON*this.chargeValue)*Fmath.N_AVAGADRO*this.electrolyteConcn/(Fmath.EPSILON_0*this.epsilon*Fmath.K_BOLTZMANN*this.tempK));
double expPart = Math.exp(this.expTerm*this.chargeValue*psiL/2.0D);
double gamma = (expPart - 1.0D)/(expPart + 1.0D);
double gammaExp = gamma*Math.exp(-kappa*xDistance);
potAtX = 2.0D*Math.log((1.0D + gammaExp)/(1.0D - gammaExp))/(this.expTerm*this.chargeValue);
}
else{
// asymmetric electrolye
// Create instance of integration function
FunctionPatX func = new FunctionPatX();
func.numOfIons = this.numOfIons;
func.termOne = 2.0D*Fmath.N_AVAGADRO*Fmath.K_BOLTZMANN*this.tempK/(Fmath.EPSILON_0*this.epsilon);
func.expTerm = this.expTerm ;
func.bulkConcn = this.bulkConcn;
func.charges = this.charges;
int nPointsGQ = 2000; // number of points in the numerical integration
// bisection procedure
double psiXlow = 0.0;
double pFuncLow = xDistance - Integration.trapezium(func, psiXlow, psiL, nPointsGQ);
double psiXhigh = psiL;
double pFuncHigh = xDistance - Integration.trapezium(func, psiXhigh, psiL, nPointsGQ);
if(pFuncHigh*pFuncLow>0.0D)throw new IllegalArgumentException("root not bounded");
double check = Math.abs(xDistance)*1e-2;
boolean test = true;
double psiXmid = 0.0D;
double pFuncMid = 0.0D;
int nIter = 0;
while(test){
psiXmid = (psiXlow + psiXhigh)/2.0D;
pFuncMid = xDistance - Integration.trapezium(func, psiXmid, psiL, nPointsGQ);
if(Math.abs(pFuncMid)<=check){
potAtX = psiXmid;
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: getPotentialAtX\nnumber of iterations exceeded in outer bisection\ncurrent value of psi(x) returned");
potAtX = psiXmid;
test = false;
}
else{
if(pFuncLow*pFuncMid>0.0D){
psiXlow = psiXmid;
pFuncLow = pFuncMid;
}
else{
psiXhigh = psiXmid;
pFuncHigh = pFuncMid;
}
}
}
}
}
return potAtX;
}
// Get concentrations at a distance, x, from the interface(M)
public double[] getConcnsAtX(double xDistance){
if(!this.psi0set && !this.sigmaSet)throw new IllegalArgumentException("Neither a surface potential nor a surface charge/density have been entered");
if(this.sigmaSet && !this.psi0set)this.getSurfacePotential();
if(this.psi0set && !this.sigmaSet)this.getSurfaceChargeDensity();
double[] conc = new double[this.numOfIons];
if(this.sternOption && xDistance<this.sternDelta){
for(int i=0; i<this.numOfIons; i++)conc[i] = 0.0D;
}
else{
double psi = this.getPotentialAtX(xDistance);
for(int i=0; i<this.numOfIons; i++)conc[i] = bulkConcn[i]*Math.exp(-this.expTerm*this.charges[i]*psi);
}
return conc;
}
// Get initial concentrations (M)
public double[] getInitConcns(){
if(!this.psi0set && !this.sigmaSet)unpack();
double[] conc = this.initConcn.clone();
for(int i=0; i<this.numOfIons; i++)conc[i] *= 1e-3;
return conc;
}
// Get equilibrium bulk concentrations (M)
public double[] getBulkConcns(){
if(!this.psi0set && !this.sigmaSet)throw new IllegalArgumentException("Neither a surface potential nor a surface charge/density have been entered");
if(this.sigmaSet && !this.psi0set)this.getSurfacePotential();
if(this.psi0set && !this.sigmaSet)this.getSurfaceChargeDensity();
double[] conc = this.bulkConcn.clone();
for(int i=0; i<this.numOfIons; i++)conc[i] *= 1e-3;
return conc;
}
// Get adsorbed ion concentrations (mol m-2)
public double[] getSiteConcns(){
if(!this.psi0set && !this.sigmaSet)throw new IllegalArgumentException("Neither a surface potential nor a surface charge/density have been entered");
if(this.sigmaSet && !this.psi0set)this.getSurfacePotential();
if(this.psi0set && !this.sigmaSet)this.getSurfaceChargeDensity();
return this.siteConcn;
}
// Calculate the surface potential
// Gouy-Chapman-Stern - Specific absorption - symmetric electrolyte
private double surfacePotential2(double sigma){
double surPot = 0.0D; // temporary values of the surface potential
// INITIAL ESTIMATES
// Ignore specific absorption
double diffPot = Fmath.asinh(sigma/Math.sqrt(this.eightTerm*this.electrolyteConcn))/(this.chargeValue*this.expTermOver2);
// Recalculate total charge with adsorption
// double totCharge = getSurfaceChargeDensity(diffPot);
double pLow = 0.0D;
double pHigh = 0.0D;
if(diffPot>0.0D){
pHigh = 2.0D*diffPot;
}
else{
pLow = 2.0D*diffPot;
}
// call bisection method
surPot = surfacePotentialBisection(pLow, pHigh, sigma, 2);
return surPot;
}
// Bisection method for surfacePotential2 and surfacePotential3
// n = 2 for surfacePotential2 and n = 3 for surfacePotential3
private double surfacePotentialBisection(double pLow, double pHigh, double sigma, int n){
double surPot = 0.0D; // temporary values of the surface potential
// Bisection method converging on sigma calculated = true sigma
boolean testC = true;
int testCiter = 0;
int testCmax = 10;
double pDiff = pHigh-pLow;
double cLow = cFunction(pLow, sigma);
double cHigh = cFunction(pHigh, sigma);
while(testC){
if(pHigh*pLow>0.0D){
testCiter++;
if(testCiter>testCmax)throw new IllegalArgumentException("root not bounded after " + testCiter + "expansions");
pLow -= pDiff;
pHigh += pDiff;
cLow = cFunction(pLow, sigma);
cHigh = cFunction(pHigh, sigma);
}
else{
testC=false;
}
}
double check = Math.abs(sigma)*1e-6;
boolean test = true;
double cMid = 0.0D;
int nIter = 0;
while(test){
surPot = (pLow + pHigh)/2.0D;
cMid = this.cFunction(surPot, sigma);
if(Math.abs(cMid)<=check){
test = false;
}
else{
nIter++;
if(nIter>10000){
System.out.println("Class: GouyChapmanStern\nMethod: surfacePotential"+n +"\nnumber of iterations exceeded in bisection\ncurrent value of sigma returned");
test = false;
}
if(cMid*cHigh>0){
pHigh = surPot;
cHigh = cMid;
}
else{
pLow = surPot;
cLow = cMid;
}
}
}
return surPot;
}
// Function to calculate difference between set value of the charge density and the estimated value for a current estimate of the surface potential
private double cFunction(double psi, double sigma){
double sigmaEst = this.calcSurfaceChargeDensity(psi);
return sigmaEst - sigma;
}
// Gouy-Chapman-Stern - Specific absorption - asymmetric electrolyte
private double surfacePotential3(double sigma){
double surPot = 0.0D; // temporary values of the surface potential
// INITIAL ESTIMATES
// Ignore specific absorption
double diffPot = Fmath.asinh(sigma/Math.sqrt(this.eightTerm*this.electrolyteConcn))/(this.chargeValue*this.expTermOver2);
// Recalculate total charge with adsorption
// double totCharge = getSurfaceChargeDensity(diffPot);
double pLow = 0.0D;
double pHigh = 0.0D;
if(diffPot>0.0D){
pHigh = 2.0D*diffPot;
}
else{
pLow = 2.0D*diffPot;
}
// call bisection method
surPot = surfacePotentialBisection(pLow, pHigh, sigma, 3);
return surPot;
}
}