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;
    }