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