void setValue(float value) {
			if (this.value < value)
				return;
			if (logCube)
				Logger.info("voxel.setValue " + x + " " + y + " " + z + ptXyz + ": " + value + " was " + this.value);
			this.value = value;
		}
	}

	int setVoxelRange(int index, float min, float max, float ptsPerAngstrom, int gridMax) {
		float range = max - min;
		int nGrid;
		if (resolution != Float.MAX_VALUE) {
			ptsPerAngstrom = resolution;
			nGrid = (int) (range * ptsPerAngstrom);
		}
		else {
			nGrid = (int) (range * ptsPerAngstrom);
		}
		if (nGrid > gridMax) {
			if ((dataType & HAS_MAXGRID) > 0) {
				if (resolution != Float.MAX_VALUE)
					Logger.info("Maximum number of voxels for index=" + index);
				nGrid = gridMax;
			}
			else if (resolution == Float.MAX_VALUE) {
				nGrid = gridMax;
			}
		}
		ptsPerAngstrom = nGrid / range;
		float d = volumetricVectorLengths[index] = 1f / ptsPerAngstrom;
		voxelCounts[index] = nGrid + ((dataType & IS_SOLVENTTYPE) != 0 ? 3 : 0);

		switch (index) {
		case 0:
			volumetricVectors[0].set(d, 0, 0);
			volumetricOrigin.x = min;
			break;
		case 1:
			volumetricVectors[1].set(0, d, 0);
			volumetricOrigin.y = min;
			break;
		case 2:
			volumetricVectors[2].set(0, 0, d);
			volumetricOrigin.z = min;
			if (isEccentric)
				eccentricityMatrix.transform(volumetricOrigin);
			if (center.x != Float.MAX_VALUE)
				volumetricOrigin.add(center);
		}
		if (isEccentric)
			eccentricityMatrix.transform(volumetricVectors[index]);
		unitVolumetricVectors[index].normalize(volumetricVectors[index]);
		return voxelCounts[index];
	}

	String jvxlGetVolumeHeader(int nAtoms) {
		String str = (-nAtoms) + " " + (volumetricOrigin.x / ANGSTROMS_PER_BOHR) + " "
				+ (volumetricOrigin.y / ANGSTROMS_PER_BOHR) + " " + (volumetricOrigin.z / ANGSTROMS_PER_BOHR) + "\n";
		for (int i = 0; i < 3; i++)
			str += voxelCounts[i] + " " + (volumetricVectors[i].x / ANGSTROMS_PER_BOHR) + " "
					+ (volumetricVectors[i].y / ANGSTROMS_PER_BOHR) + " "
					+ (volumetricVectors[i].z / ANGSTROMS_PER_BOHR) + "\n";
		return str;
	}

	// just for small factorials
	static float factorial(int n) {
		if (n == 0)
			return 1;
		return n * factorial(n - 1);
	}

	float[] fact = new float[20];

	float getDefaultResolution() {
		return Float.MAX_VALUE;
		// for popup menu?
		/*
		 * maybe.... int nAtoms = viewer.getAtomCount(); float res = 0; if (nAtoms < 200) res = 4; else if (nAtoms <
		 * 2000) res = 1; return res;
		 */
	}

	// ///// spheres and ellipsoids //////

	int sphere_gridMax = 20;
	float sphere_ptsPerAngstrom = 10f;
	float sphere_radiusAngstroms;

	void setupSphere() {
		if (center.x == Float.MAX_VALUE)
			center.set(0, 0, 0);
		float radius = sphere_radiusAngstroms * 1.1f * eccentricityScale;
		for (int i = 0; i < 3; i++)
			setVoxelRange(i, -radius, radius, sphere_ptsPerAngstrom, sphere_gridMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("SPHERE \nres="
				+ sphere_ptsPerAngstrom
				+ " rad="
				+ sphere_radiusAngstroms
				+ (isAnisotropic ? " anisotropy=(" + anisotropy[0] + "," + anisotropy[1] + "," + anisotropy[2] + ")"
						: "") + "\n");
		jvxlFileHeader.append(jvxlGetVolumeHeader(2));
		atomCount = 0;
		negativeAtomCount = false;
	}

	float getSphereValue(int x, int y, int z) {
		voxelPtToXYZ(x, y, z, ptPsi);
		getCalcPoint(ptPsi);
		return sphere_radiusAngstroms - (float) Math.sqrt(ptPsi.x * ptPsi.x + ptPsi.y * ptPsi.y + ptPsi.z * ptPsi.z);
	}

	// ///// hydrogen-like Schroedinger orbitals ///////

	int psi_gridMax = 40;
	float psi_ptsPerAngstrom = 5f;
	float psi_radiusAngstroms;

	double[] rfactor = new double[10];
	double[] pfactor = new double[10];
	int lastFactorial = -1;

	final static double A0 = 0.52918f; // x10^-10 meters

	final static double ROOT2 = 1.414214;
	int psi_n = 2;
	int psi_l = 1;
	int psi_m = 1;
	float psi_Znuc = 1; // hydrogen

	void setupOrbital() {
		psi_radiusAngstroms = autoScaleOrbital();
		if (center.x == Float.MAX_VALUE)
			center.set(0, 0, 0);
		for (int i = 0; i < 3; i++)
			setVoxelRange(i, -psi_radiusAngstroms, psi_radiusAngstroms, psi_ptsPerAngstrom, psi_gridMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("hydrogen-like orbital \nn="
				+ psi_n
				+ ", l="
				+ psi_l
				+ ", m="
				+ psi_m
				+ " Znuc="
				+ psi_Znuc
				+ " res="
				+ psi_ptsPerAngstrom
				+ " rad="
				+ psi_radiusAngstroms
				+ (isAnisotropic ? " anisotropy=(" + anisotropy[0] + "," + anisotropy[1] + "," + anisotropy[2] + ")"
						: "") + "\n");
		jvxlFileHeader.append(jvxlGetVolumeHeader(2));
		atomCount = 0;
		negativeAtomCount = false;
		calcFactors(psi_n, psi_l, psi_m);
	}

	float autoScaleOrbital() {
		float w = (psi_n * (psi_n + 3) - 5f) / psi_Znuc;
		if (w < 1)
			w = 1;
		if (psi_n < 3)
			w += 1;
		float aMax = 0;
		if (!isAnisotropic)
			return w;
		for (int i = 3; --i >= 0;)
			if (anisotropy[i] > aMax)
				aMax = anisotropy[i];
		return w * aMax;
	}

	void calcFactors(int n, int el, int m) {
		int abm = Math.abs(m);
		if (lastFactorial < n + el) {
			for (int i = lastFactorial + 1; i <= n + el; i++)
				fact[i] = factorial(i);
			lastFactorial = n + el;
		}
		double Nnl = Math.pow(2 * psi_Znuc / n / A0, 1.5)
				* Math.sqrt(fact[n - el - 1] / 2 / n / Math.pow(fact[n + el], 3));
		double Lnl = fact[n + el] * fact[n + el];
		double Plm = Math.pow(2, -el) * fact[el] * fact[el + abm]
				* Math.sqrt((2 * el + 1) * fact[el - abm] / 2 / fact[el + abm]);

		for (int p = 0; p <= n - el - 1; p++)
			rfactor[p] = Nnl * Lnl / fact[p] / fact[n - el - p - 1] / fact[2 * el + p + 1];
		for (int p = abm; p <= el; p++)
			pfactor[p] = Math.pow(-1, el - p) * Plm / fact[p] / fact[el + abm - p] / fact[el - p] / fact[p - abm];
	}

	final Point3f ptPsi = new Point3f();

	float getPsi(int x, int y, int z) {
		voxelPtToXYZ(x, y, z, ptPsi);
		getCalcPoint(ptPsi);
		return (float) hydrogenAtomPsiAt(ptPsi, psi_n, psi_l, psi_m);
	}

	double hydrogenAtomPsiAt(Point3f pt, int n, int el, int m) {
		// ref: http://www.stolaf.edu/people/hansonr/imt/concept/schroed.pdf
		int abm = Math.abs(m);
		double x2y2 = pt.x * pt.x + pt.y * pt.y;
		double r2 = x2y2 + pt.z * pt.z;
		double r = Math.sqrt(r2);
		double rho = 2d * psi_Znuc * r / n / A0;
		double ph, th, cth, sth;
		double theta_lm = 0;
		double phi_m = 0;
		double sum = 0;
		for (int p = 0; p <= n - el - 1; p++)
			sum += Math.pow(-rho, p) * rfactor[p];
		double rnl = Math.exp(-rho / 2) * Math.pow(rho, el) * sum;
		ph = Math.atan2(pt.y, pt.x);
		th = Math.atan2(Math.sqrt(x2y2), pt.z);
		cth = Math.cos(th);
		sth = Math.sin(th);
		sum = 0;
		for (int p = abm; p <= el; p++)
			sum += Math.pow(1 + cth, p - abm) * Math.pow(1 - cth, el - p) * pfactor[p];
		theta_lm = Math.abs(Math.pow(sth, abm)) * sum;
		if (m == 0)
			phi_m = 1;
		else if (m > 0)
			phi_m = Math.cos(m * ph) * ROOT2;
		else phi_m = Math.sin(-m * ph) * ROOT2;
		if (Math.abs(phi_m) < 0.0000000001)
			phi_m = 0;
		return rnl * theta_lm * phi_m;
	}

	// ////// lobes ////////

	float lobe_sizeAngstroms = 1.0f;
	int lobe_gridMax = 21;
	int lobe_ptsPerAngstrom = 10;

	void setupLobe() {
		psi_n = 3;
		psi_l = 2;
		psi_m = 0;
		// using Hens Borkent's 3dz2 lobe idea
		psi_Znuc = 15;
		if (center.x == Float.MAX_VALUE)
			center.set(0, 0, 0);
		float radius = lobe_sizeAngstroms * 1.1f * eccentricityRatio * eccentricityScale;
		if (eccentricityScale > 0 && eccentricityScale < 1)
			radius /= eccentricityScale;
		setVoxelRange(0, -radius, radius, lobe_ptsPerAngstrom, lobe_gridMax);
		setVoxelRange(1, -radius, radius, lobe_ptsPerAngstrom, lobe_gridMax);
		setVoxelRange(2, 0, radius / eccentricityRatio, lobe_ptsPerAngstrom, lobe_gridMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("lobe \nn="
				+ psi_n
				+ ", l="
				+ psi_l
				+ ", m="
				+ psi_m
				+ " Znuc="
				+ psi_Znuc
				+ " res="
				+ lobe_ptsPerAngstrom
				+ " rad="
				+ radius
				+ (isAnisotropic ? " anisotropy=(" + anisotropy[0] + "," + anisotropy[1] + "," + anisotropy[2] + ")"
						: "") + "\n");
		jvxlFileHeader.append(jvxlGetVolumeHeader(2));
		atomCount = 0;
		negativeAtomCount = false;
		calcFactors(psi_n, psi_l, psi_m);
	}

	float getLobeValue(int x, int y, int z) {
		voxelPtToXYZ(x, y, z, ptPsi);
		getCalcPoint(ptPsi);
		float value = (float) hydrogenAtomPsiAt(ptPsi, psi_n, psi_l, psi_m);
		if (value < 0)
			value = 0;
		return value;
	}

	// ///// ab initio/semiempirical quantum mechanical orbitals ///////

	int qm_gridMax = QuantumCalculation.MAX_GRID;
	float qm_ptsPerAngstrom = 10f;
	float qm_marginAngstroms = 1f; // may have to adjust this
	int qm_nAtoms;
	int qm_moNumber = Integer.MAX_VALUE;

	Atom[] qm_atoms;

	void addMOTitleInfo(int iLine, Hashtable mo) {
		String line = title[iLine];
		int pt = line.indexOf("%");
		if (line.length() == 0 || pt < 0)
			return;
		boolean replaced = false;
		for (int i = pt; i < line.length() - 1; i++) {
			if (line.charAt(i) == '%') {
				String info = "";
				switch (line.charAt(i + 1)) {
				case 'F':
					info = viewer.getFileName();
					break;
				case 'I':
					info += qm_moNumber;
					break;
				case 'N':
					info += qmOrbitalCount;
					break;
				case 'M':
					info += viewer.getModelNumber(viewer.getDisplayModelIndex());
					break;
				case 'E':
					info += mo.get("energy");
					break;
				case 'U':
					if (mo.containsKey("energyUnits"))
						info += moData.get("energyUnits");
					break;
				case 'S':
					if (mo.containsKey("symmetry"))
						info += mo.get("symmetry");
					break;
				case 'O':
					if (mo.containsKey("occupancy"))
						info += mo.get("occupancy");
					break;
				}
				replaced |= (info.length() > 0);
				line = line.substring(0, i) + info + line.substring(i + 2);
				i += info.length();
			}
		}
		line = (replaced || line.charAt(0) != '?' ? line : "");
		title[iLine] = (line.length() > 1 && line.charAt(0) == '?' ? line.substring(1) : line);
	}

	void setupQMOrbital() {
		Atom[] atoms = frame.atoms;
		Point3f xyzMin = new Point3f(Float.MAX_VALUE, Float.MAX_VALUE, Float.MAX_VALUE);
		Point3f xyzMax = new Point3f(-Float.MAX_VALUE, -Float.MAX_VALUE, -Float.MAX_VALUE);
		int modelIndex = viewer.getDisplayModelIndex();
		int iAtom = 0;
		int nSelected = 0;
		int nAtoms = viewer.getAtomCount();
		iUseBitSets = true;
		for (int i = 0; i < nAtoms; i++)
			if (atoms[i].modelIndex == modelIndex) {
				++iAtom;
				if (bsSelected == null || bsSelected.get(i))
					++nSelected;
			}
		qm_nAtoms = iAtom;
		if (nSelected > 0)
			Logger.info(nSelected + " of " + qm_nAtoms + " atoms will be used in the orbital calculation");
		if (qm_nAtoms > 0)
			qm_atoms = new Atom[qm_nAtoms];
		iAtom = 0;
		for (int i = 0; i < nAtoms; i++) {
			Atom atom = atoms[i];
			if (atom.modelIndex != modelIndex)
				continue;
			Point3f pt = new Point3f(atom);
			if (nSelected == 0 || bsSelected == null || bsSelected.get(i)) {
				float rA = atom.getVanderwaalsRadiusFloat() + qm_marginAngstroms;
				if (pt.x - rA < xyzMin.x)
					xyzMin.x = pt.x - rA;
				if (pt.x + rA > xyzMax.x)
					xyzMax.x = pt.x + rA;
				if (pt.y - rA < xyzMin.y)
					xyzMin.y = pt.y - rA;
				if (pt.y + rA > xyzMax.y)
					xyzMax.y = pt.y + rA;
				if (pt.z - rA < xyzMin.z)
					xyzMin.z = pt.z - rA;
				if (pt.z + rA > xyzMax.z)
					xyzMax.z = pt.z + rA;
				qm_atoms[iAtom++] = atom;
			}
			else {
				++iAtom;
			}
		}
		if (!Float.isNaN(scale)) {
			Vector3f v = new Vector3f(xyzMax);
			v.sub(xyzMin);
			v.scale(0.5f);
			xyzMin.add(v);
			v.scale(scale);
			xyzMax.set(xyzMin);
			xyzMax.add(v);
			xyzMin.sub(v);
		}

		Logger.info("MO range bohr " + xyzMin + " to " + xyzMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("MO range bohr " + xyzMin + " to " + xyzMax + "\ncalculation type: "
				+ moData.get("calculationType") + "\n");

		int maxGrid = qm_gridMax;

		setVoxelRange(0, xyzMin.x, xyzMax.x, qm_ptsPerAngstrom, maxGrid);
		setVoxelRange(1, xyzMin.y, xyzMax.y, qm_ptsPerAngstrom, maxGrid);
		setVoxelRange(2, xyzMin.z, xyzMax.z, qm_ptsPerAngstrom, maxGrid);
		jvxlFileHeader.append(jvxlGetVolumeHeader(iAtom));
		Point3f pt = new Point3f();
		for (int i = 0; i < nAtoms; i++) {
			Atom atom = atoms[i];
			if (atom.modelIndex != modelIndex)
				continue;
			pt.set(atom);
			pt.scale(1 / ANGSTROMS_PER_BOHR);
			jvxlFileHeader.append(atom.getAtomicAndIsotopeNumber() + " " + atom.getAtomicAndIsotopeNumber() + ".0 "
					+ pt.x + " " + pt.y + " " + pt.z + "\n");
		}
		atomCount = -Integer.MAX_VALUE;
		negativeAtomCount = false;
		precalculateVoxelData = true;
	}

	void generateQuantumCube() {
		QuantumCalculation q;
		float[] origin = { volumetricOrigin.x, volumetricOrigin.y, volumetricOrigin.z };
		switch (qmOrbitalType) {
		case QM_TYPE_GAUSSIAN:
			q = new QuantumCalculation((String) moData.get("calculationType"), qm_atoms, (Vector) moData.get("shells"),
					(float[][]) moData.get("gaussians"), (Hashtable) moData.get("atomicOrbitalOrder"), null, null,
					moCoefficients);
			q.createGaussianCube(voxelData, voxelCounts, origin, volumetricVectorLengths);
			break;
		case QM_TYPE_SLATER:
			q = new QuantumCalculation((String) moData.get("calculationType"), qm_atoms, (Vector) moData.get("shells"),
					null, null, (int[][]) moData.get("slaterInfo"), (float[][]) moData.get("slaterData"),
					moCoefficients);
			q.createSlaterCube(voxelData, voxelCounts, origin, volumetricVectorLengths);
			break;
		default:
		}
	}

	// ///// molecular electrostatic potential ///////

	int mep_gridMax = MepCalculation.MAX_GRID;
	float mep_ptsPerAngstrom = 3f;
	float mep_marginAngstroms = 1f; // may have to adjust this
	int mep_nAtoms;

	Atom[] mep_atoms;

	void setupMep() {
		Atom[] atoms = frame.atoms;
		Point3f xyzMin = new Point3f(Float.MAX_VALUE, Float.MAX_VALUE, Float.MAX_VALUE);
		Point3f xyzMax = new Point3f(-Float.MAX_VALUE, -Float.MAX_VALUE, -Float.MAX_VALUE);
		int modelIndex = viewer.getDisplayModelIndex();
		int iAtom = 0;
		int nSelected = 0;
		iUseBitSets = true;
		int nAtoms = viewer.getAtomCount();
		for (int i = 0; i < nAtoms; i++)
			if (atoms[i].modelIndex == modelIndex) {
				++iAtom;
				if (bsSelected.get(i))
					++nSelected;
			}
		mep_nAtoms = iAtom;
		if (nSelected > 0)
			Logger.info(nSelected + " of " + mep_nAtoms + " atoms will be used in the mep calculation");
		if (mep_nAtoms > 0)
			mep_atoms = new Atom[mep_nAtoms];
		iAtom = 0;
		for (int i = 0; i < nAtoms; i++) {
			Atom atom = atoms[i];
			if (atom.modelIndex != modelIndex)
				continue;
			Point3f pt = new Point3f(atom);
			if (bsSelected.get(i)) {
				float rA = atom.getVanderwaalsRadiusFloat() + mep_marginAngstroms;
				if (pt.x - rA < xyzMin.x)
					xyzMin.x = pt.x - rA;
				if (pt.x + rA > xyzMax.x)
					xyzMax.x = pt.x + rA;
				if (pt.y - rA < xyzMin.y)
					xyzMin.y = pt.y - rA;
				if (pt.y + rA > xyzMax.y)
					xyzMax.y = pt.y + rA;
				if (pt.z - rA < xyzMin.z)
					xyzMin.z = pt.z - rA;
				if (pt.z + rA > xyzMax.z)
					xyzMax.z = pt.z + rA;
				mep_atoms[iAtom++] = atom;
			}
			else {
				++iAtom;
			}
		}
		if (!Float.isNaN(scale)) {
			Vector3f v = new Vector3f(xyzMax);
			v.sub(xyzMin);
			v.scale(0.5f);
			xyzMin.add(v);
			v.scale(scale);
			xyzMax.set(xyzMin);
			xyzMax.add(v);
			xyzMin.sub(v);
		}

		Logger.info("MEP range bohr " + xyzMin + " to " + xyzMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("MEP range bohr " + xyzMin + " to " + xyzMax + "\n");

		int maxGrid = mep_gridMax;

		setVoxelRange(0, xyzMin.x, xyzMax.x, mep_ptsPerAngstrom, maxGrid);
		setVoxelRange(1, xyzMin.y, xyzMax.y, mep_ptsPerAngstrom, maxGrid);
		setVoxelRange(2, xyzMin.z, xyzMax.z, mep_ptsPerAngstrom, maxGrid);
		jvxlFileHeader.append(jvxlGetVolumeHeader(iAtom));
		Point3f pt = new Point3f();
		for (int i = 0; i < nAtoms; i++) {
			Atom atom = atoms[i];
			if (atom.modelIndex != modelIndex)
				continue;
			pt.set(atom);
			pt.scale(1 / ANGSTROMS_PER_BOHR);
			jvxlFileHeader.append(atom.getAtomicAndIsotopeNumber() + " " + atom.getAtomicAndIsotopeNumber() + ".0 "
					+ pt.x + " " + pt.y + " " + pt.z + "\n");
		}
		atomCount = -Integer.MAX_VALUE;
		negativeAtomCount = false;
		precalculateVoxelData = true;
	}

	void generateMepCube() {
		float[] origin = { volumetricOrigin.x, volumetricOrigin.y, volumetricOrigin.z };
		MepCalculation m = new MepCalculation(mep_atoms, mepCharges);
		m.createMepCube(voxelData, voxelCounts, origin, volumetricVectorLengths);
	}

	// /// solvent-accessible, solvent-excluded surface //////

	float solvent_ptsPerAngstrom = 4f;
	int solvent_gridMax = 60;
	int solvent_modelIndex;
	float[] solvent_atomRadius;
	Point3f[] solvent_ptAtom;
	int solvent_nAtoms;
	int solvent_firstNearbyAtom;
	boolean solvent_quickPlane;
	BitSet atomSet = new BitSet();
	BitSet bsSolventSelected;
	Voxel solvent_voxel = new Voxel();

	void setupSolvent() {
		/*
		 * The surface fragment idea:
		 * 
		 * isosurface solvent|sasurface both work on the SELECTED atoms, thus allowing for a subset of the molecule to
		 * be involved. But in that case we don't want to be creating a surface that goes right through another atom.
		 * Rather, what we want (probably) is just the portion of the OVERALL surface that involves these atoms.
		 * 
		 * The addition of Mesh.voxelValue[] means that we can specify any voxel we want to NOT be excluded (NaN). Here
		 * we first exclude any voxel that would have been INSIDE a nearby atom. This will take care of any portion of
		 * the vanderwaals surface that would be there. Then we exclude any special-case voxel that is between two
		 * nearby atoms.
		 * 
		 * Bob Hanson 13 Jul 2006
		 * 
		 */
		bsSolventSelected = new BitSet();
		if (thePlane != null)
			setPlaneParameters(thePlane);
		solvent_quickPlane = true;// viewer.getTestFlag1();
		Atom[] atoms = frame.atoms;
		Point3f xyzMin = new Point3f(Float.MAX_VALUE, Float.MAX_VALUE, Float.MAX_VALUE);
		Point3f xyzMax = new Point3f(-Float.MAX_VALUE, -Float.MAX_VALUE, -Float.MAX_VALUE);
		solvent_modelIndex = -1;
		int iAtom = 0;
		int nAtoms = viewer.getAtomCount();
		iUseBitSets = true;
		if (bsIgnore == null) {
			bsIgnore = new BitSet();
		}
		for (int i = 0; i < nAtoms; i++) {
			if (bsSelected.get(i) && (!bsIgnore.get(i))) {
				if (solvent_quickPlane && thePlane != null
						&& Math.abs(distancePointToPlane(atoms[i], thePlane)) > 2 * solventWorkingRadius(atoms[i])) {
					continue;
				}
				bsSolventSelected.set(i);
			}
		}
		atomSet = new BitSet();
		int firstSet = -1;
		int lastSet = 0;
		for (int i = 0; i < nAtoms; i++)
			if (bsSolventSelected.get(i)) {
				if (solvent_modelIndex < 0)
					solvent_modelIndex = atoms[i].modelIndex;
				if (solvent_modelIndex != atoms[i].modelIndex) {
					bsIgnore.set(i);
					continue;
				}
				++iAtom;
				atomSet.set(i);
				if (firstSet == -1)
					firstSet = i;
				lastSet = i;
			}
		int nH = 0;
		int[] atomNo = null;
		if (iAtom > 0) {
			Point3f[] hAtoms = null;
			if (addHydrogens) {
				hAtoms = viewer.getAdditionalHydrogens(atomSet);
				nH = hAtoms.length;
			}
			solvent_atomRadius = new float[iAtom + nH];
			solvent_ptAtom = new Point3f[iAtom + nH];
			atomNo = new int[iAtom + nH];

			float r = solventWorkingRadius(null);
			for (int i = 0; i < nH; i++) {
				atomNo[i] = 1;
				solvent_atomRadius[i] = r;
				solvent_ptAtom[i] = hAtoms != null ? hAtoms[i] : null;
			}
			iAtom = nH;
			for (int i = firstSet; i <= lastSet; i++) {
				if (!atomSet.get(i))
					continue;
				atomNo[iAtom] = atoms[i].getElementNumber();
				solvent_ptAtom[iAtom] = atoms[i];
				solvent_atomRadius[iAtom++] = solventWorkingRadius(atoms[i]);
			}
		}
		solvent_nAtoms = solvent_firstNearbyAtom = iAtom;
		Logger.info(iAtom + " atoms will be used in the solvent-accessible surface calculation");

		for (int i = 0; i < solvent_nAtoms; i++) {
			Point3f pt = solvent_ptAtom[i];
			float rA = solvent_atomRadius[i];
			if (pt.x - rA < xyzMin.x)
				xyzMin.x = pt.x - rA;
			if (pt.x + rA > xyzMax.x)
				xyzMax.x = pt.x + rA;
			if (pt.y - rA < xyzMin.y)
				xyzMin.y = pt.y - rA;
			if (pt.y + rA > xyzMax.y)
				xyzMax.y = pt.y + rA;
			if (pt.z - rA < xyzMin.z)
				xyzMin.z = pt.z - rA;
			if (pt.z + rA > xyzMax.z)
				xyzMax.z = pt.z + rA;
		}
		Logger.info("surface range " + xyzMin + " to " + xyzMax);
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("solvent-" + (dataType == SURFACE_SASURFACE ? "accesible" : "excluded")
				+ " surface\nrange " + xyzMin + " to " + xyzMax + "\n");

		// fragment idea

		Point3f pt = new Point3f();

		BitSet bsNearby = new BitSet();
		int nNearby = 0;
		firstSet = -1;
		lastSet = 0;
		for (int i = 0; i < nAtoms; i++) {
			if (atomSet.get(i) || bsIgnore.get(i))
				continue;
			float rA = solventWorkingRadius(atoms[i]);
			if (solvent_quickPlane && thePlane != null && Math.abs(distancePointToPlane(atoms[i], thePlane)) > 2 * rA)
				continue;
			pt = atoms[i];
			if (pt.x + rA > xyzMin.x && pt.x - rA < xyzMax.x && pt.y + rA > xyzMin.y && pt.y - rA < xyzMax.y
					&& pt.z + rA > xyzMin.z && pt.z - rA < xyzMax.z) {
				if (firstSet == -1)
					firstSet = i;
				lastSet = i;
				bsNearby.set(i);
				nNearby++;
			}
		}

		if (nNearby != 0) {
			solvent_nAtoms += nNearby;
			solvent_atomRadius = ArrayUtil.setLength(solvent_atomRadius, solvent_nAtoms);
			solvent_ptAtom = (Point3f[]) ArrayUtil.setLength(solvent_ptAtom, solvent_nAtoms);

			iAtom = solvent_firstNearbyAtom;
			for (int i = firstSet; i <= lastSet; i++) {
				if (!bsNearby.get(i))
					continue;
				solvent_ptAtom[iAtom] = atoms[i];
				solvent_atomRadius[iAtom++] = solventWorkingRadius(atoms[i]);
			}
		}

		int maxGrid;
		maxGrid = solvent_gridMax;
		setVoxelRange(0, xyzMin.x, xyzMax.x, solvent_ptsPerAngstrom, maxGrid);
		setVoxelRange(1, xyzMin.y, xyzMax.y, solvent_ptsPerAngstrom, maxGrid);
		setVoxelRange(2, xyzMin.z, xyzMax.z, solvent_ptsPerAngstrom, maxGrid);
		precalculateVoxelData = newSolventMethod;
		int nAtomsWritten = Math.min(solvent_firstNearbyAtom, 100);
		jvxlFileHeader.append(jvxlGetVolumeHeader(nAtomsWritten));

		pt = new Point3f();
		for (int i = 0; i < nAtomsWritten; i++) {
			pt.set(solvent_ptAtom[i]);
			pt.scale(1 / ANGSTROMS_PER_BOHR);
			if (atomNo != null)
				jvxlFileHeader.append(atomNo[i] + " " + atomNo[i] + ".0 " + pt.x + " " + pt.y + " " + pt.z + "\n");
		}
		atomCount = -Integer.MAX_VALUE;
		negativeAtomCount = false;
	}

	float solventWorkingRadius(Atom atom) {
		float r = (solventAtomRadiusAbsolute > 0 ? solventAtomRadiusAbsolute
				: atom == null ? JmolConstants.vanderwaalsMars[1] / 1000f : useIonic ? atom.getBondingRadiusFloat()
						: atom.getVanderwaalsRadiusFloat());
		r *= solventAtomRadiusFactor;
		r += solventExtendedAtomRadius + solventAtomRadiusOffset;
		if (r < 0.1)
			r = 0.1f;
		return r;
	}

	void generateSolventCube() {
		long time = System.currentTimeMillis();
		float rA, rB;
		Point3f ptA;
		Point3f ptY0 = new Point3f(), ptZ0 = new Point3f();
		Point3i pt0 = new Point3i(), pt1 = new Point3i();
		float maxValue = (dataType == SURFACE_NOMAP ? Float.MAX_VALUE : Float.MAX_VALUE);
		for (int x = 0; x < nPointsX; ++x)
			for (int y = 0; y < nPointsY; ++y)
				for (int z = 0; z < nPointsZ; ++z)
					voxelData[x][y][z] = maxValue;
		if (dataType == SURFACE_NOMAP)
			return;
		float maxRadius = 0;
		for (int iAtom = 0; iAtom < solvent_nAtoms; iAtom++) {
			ptA = solvent_ptAtom[iAtom];
			rA = solvent_atomRadius[iAtom];
			if (rA > maxRadius)
				maxRadius = rA;
			boolean isNearby = (iAtom >= solvent_firstNearbyAtom);
			setGridLimitsForAtom(ptA, rA, pt0, pt1);
			voxelPtToXYZ(pt0.x, pt0.y, pt0.z, ptXyzTemp);
			for (int i = pt0.x; i < pt1.x; i++) {
				ptY0.set(ptXyzTemp);
				for (int j = pt0.y; j < pt1.y; j++) {
					ptZ0.set(ptXyzTemp);
					for (int k = pt0.z; k < pt1.z; k++) {
						float v = ptXyzTemp.distance(ptA) - rA;
						if (v < voxelData[i][j][k])
							voxelData[i][j][k] = (isNearby ? Float.NaN : v);
						ptXyzTemp.add(volumetricVectors[2]);
					}
					ptXyzTemp.set(ptZ0);
					ptXyzTemp.add(volumetricVectors[1]);
				}
				ptXyzTemp.set(ptY0);
				ptXyzTemp.add(volumetricVectors[0]);
			}
		}
		if ((dataType == SURFACE_SOLVENT || dataType == SURFACE_MOLECULAR) && solventRadius > 0) {
			Point3i ptA0 = new Point3i();
			Point3i ptB0 = new Point3i();
			Point3i ptA1 = new Point3i();
			Point3i ptB1 = new Point3i();
			for (int iAtom = 0; iAtom < solvent_firstNearbyAtom - 1; iAtom++)
				if (solvent_ptAtom[iAtom] instanceof Atom) {
					ptA = solvent_ptAtom[iAtom];
					rA = solvent_atomRadius[iAtom] + solventRadius;
					setGridLimitsForAtom(ptA, rA - solventRadius, ptA0, ptA1);
					AtomIterator iter = frame.getWithinModelIterator((Atom) ptA, rA + solventRadius + maxRadius);
					while (iter.hasNext()) {
						Atom ptB = iter.next();
						if (ptB.atomIndex <= ((Atom) ptA).atomIndex)
							continue;
						// selected
						// only consider selected neighbors
						if (!bsSolventSelected.get(ptB.atomIndex))
							continue;
						rB = solventWorkingRadius(ptB) + solventRadius;
						if (solvent_quickPlane && thePlane != null
								&& Math.abs(distancePointToPlane(ptB, thePlane)) > 2 * rB)
							continue;

						float dAB = ptA.distance(ptB);
						if (dAB >= rA + rB)
							continue;
						// defining pt0 and pt1 very crudely -- this could be refined
						setGridLimitsForAtom(ptB, rB - solventRadius, ptB0, ptB1);
						pt0.x = Math.min(ptA0.x, ptB0.x);
						pt0.y = Math.min(ptA0.y, ptB0.y);
						pt0.z = Math.min(ptA0.z, ptB0.z);
						pt1.x = Math.max(ptA1.x, ptB1.x);
						pt1.y = Math.max(ptA1.y, ptB1.y);
						pt1.z = Math.max(ptA1.z, ptB1.z);
						voxelPtToXYZ(pt0.x, pt0.y, pt0.z, ptXyzTemp);
						for (int i = pt0.x; i < pt1.x; i++) {
							ptY0.set(ptXyzTemp);
							for (int j = pt0.y; j < pt1.y; j++) {
								ptZ0.set(ptXyzTemp);
								for (int k = pt0.z; k < pt1.z; k++) {
									float dVS = checkSpecialVoxel(ptA, rA, ptB, rB, dAB, ptXyzTemp);
									if (!Float.isNaN(dVS)) {
										float v = solventRadius - dVS;
										if (v < voxelData[i][j][k])
											voxelData[i][j][k] = v;
									}
									ptXyzTemp.add(volumetricVectors[2]);
								}
								ptXyzTemp.set(ptZ0);
								ptXyzTemp.add(volumetricVectors[1]);
							}
							ptXyzTemp.set(ptY0);
							ptXyzTemp.add(volumetricVectors[0]);
						}
					}
				}
		}
		if (thePlane != null) {
			maxValue = 0.001f;
			for (int x = 0; x < nPointsX; ++x)
				for (int y = 0; y < nPointsY; ++y)
					for (int z = 0; z < nPointsZ; ++z)
						if (voxelData[x][y][z] < maxValue) {
							// Float.NaN will also match ">=" this way
						}
						else {
							voxelData[x][y][z] = maxValue;
						}
		}
		Logger.debug("solvent surface time:" + (System.currentTimeMillis() - time));
	}

	void setGridLimitsForAtom(Point3f ptA, float rA, Point3i pt0, Point3i pt1) {
		xyzToVoxelPt(ptA.x - rA, ptA.y - rA, ptA.z - rA, pt0);
		pt0.x -= 1;
		pt0.y -= 1;
		pt0.z -= 1;
		if (pt0.x < 0)
			pt0.x = 0;
		if (pt0.y < 0)
			pt0.y = 0;
		if (pt0.z < 0)
			pt0.z = 0;
		xyzToVoxelPt(ptA.x + rA, ptA.y + rA, ptA.z + rA, pt1);
		pt1.x += 2;
		pt1.y += 2;
		pt1.z += 2;
		if (pt1.x >= nPointsX)
			pt1.x = nPointsX;
		if (pt1.y >= nPointsY)
			pt1.y = nPointsY;
		if (pt1.z >= nPointsZ)
			pt1.z = nPointsZ;
	}

	float getSolventValue(int x, int y, int z) {

		// old method -- not used

		solvent_voxel.setValue(x, y, z, Float.MAX_VALUE);
		float rA, rB;
		Point3f ptA, ptB;
		for (int i = 0; i < solvent_nAtoms && solvent_voxel.value >= -0.5; i++) {
			ptA = solvent_ptAtom[i];
			rA = solvent_atomRadius[i];
			float v = solvent_voxel.ptXyz.distance(ptA) - rA;
			if (v < solvent_voxel.value)
				solvent_voxel.setValue(i >= solvent_firstNearbyAtom ? Float.NaN : v);
		}
		if (solventRadius == 0)
			return solvent_voxel.value;
		Point3f ptV = solvent_voxel.ptXyz;
		for (int i = 0; i < solvent_nAtoms - 1 && solvent_voxel.value >= -0.5; i++) {
			ptA = solvent_ptAtom[i];
			rA = solvent_atomRadius[i] + solventRadius;
			for (int j = i + 1; j < solvent_nAtoms && solvent_voxel.value >= -0.5; j++) {
				if (i >= solvent_firstNearbyAtom && j >= solvent_firstNearbyAtom)
					continue;
				ptB = solvent_ptAtom[j];
				rB = solvent_atomRadius[j] + solventRadius;
				float dAB = ptA.distance(ptB);
				if (dAB >= rA + rB)
					continue;
				float dVS = checkSpecialVoxel(ptA, rA, ptB, rB, dAB, ptV);
				if (!Float.isNaN(dVS))
					solvent_voxel.setValue(solventRadius - dVS);
			}
		}
		return solvent_voxel.value;
	}

	final Point3f ptS = new Point3f();

	float checkSpecialVoxel(Point3f ptA, float rAS, Point3f ptB, float rBS, float dAB, Point3f ptV) {
		/*
		 * Checking here for voxels that are in the situation:
		 * 
		 * A------)-- V ---((--))-- S --(------B |----d--------| or
		 * 
		 * B------)-- V ---((--))-- S --(------A |----d--------|
		 * 
		 * A and B are the two atom centers; V is the voxel; S is a hypothetical PROJECTED solvent center based on the
		 * position of V in relation to first A, then B; ( and ) are atom radii and (( )) are the overlapping
		 * atom+solvent radii.
		 * 
		 * That is, where the projected solvent location for one voxel is within the solvent radius sphere of another,
		 * this voxel should be checked in relation to solvent distance, not atom distance.
		 * 
		 * 
		 * S ++ / \ ++ ++ / | \ ++ + V + x want V such that angle ASV < angle ASB / ****** \ A --+--+----B b
		 * 
		 * A, B are atoms; S is solvent center; V is voxel point objective is to calculate dSV. ++ Here represents the
		 * van der Waals radius for each atom. ***** is the key "trough" location.
		 * 
		 * Getting dVS:
		 * 
		 * Known: rAB, rAS, rBS, giving angle BAS (theta) Known: rAB, rAV, rBV, giving angle VAB (alpha) Determined:
		 * angle VAS (theta - alpha), and from that, dSV, using the cosine law:
		 * 
		 * a^2 + b^2 - 2ab Cos(theta) = c^2.
		 * 
		 * The trough issue:
		 * 
		 * Since the voxel might be at point x (above), outside the triangle, we have to test for that. What we will be
		 * looking for in the "trough" will be that angle ASV < angle ASB that is, cosASB < cosASV
		 * 
		 * If we find the voxel in the "trough", then we set its value to (solvent radius - dVS).
		 * 
		 */
		float dAV = ptA.distance(ptV);
		float dBV = ptB.distance(ptV);
		float dVS = Float.NaN;
		float f = rAS / dAV;
		if (f > 1) {
			ptS.set(ptA.x + (ptV.x - ptA.x) * f, ptA.y + (ptV.y - ptA.y) * f, ptA.z + (ptV.z - ptA.z) * f);
			if (ptB.distance(ptS) < rBS) {
				dVS = solventDistance(ptV, ptA, ptB, rAS, rBS, dAB, dAV, dBV);
				if (!voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV, dBV))
					return Float.NaN;
			}
			return dVS;
		}
		f = rBS / dBV;
		if (f <= 1)
			return dVS;
		ptS.set(ptB.x + (ptV.x - ptB.x) * f, ptB.y + (ptV.y - ptB.y) * f, ptB.z + (ptV.z - ptB.z) * f);
		if (ptA.distance(ptS) < rAS) {
			dVS = solventDistance(ptV, ptB, ptA, rBS, rAS, dAB, dBV, dAV);
			if (!voxelIsInTrough(dVS, rAS * rAS, rBS, dAB, dAV, dBV))
				return Float.NaN;
		}
		return dVS;
	}

	boolean voxelIsInTrough(float dVS, float rAS2, float rBS, float dAB, float dAV, float dBV) {
		// only calculate what we need -- a factor proportional to cos
		float cosASBf = (rAS2 + rBS * rBS - dAB * dAB) / rBS; // /2 /rAS);
		float cosASVf = (rAS2 + dVS * dVS - dAV * dAV) / dVS; // /2 /rAS);
		return (cosASBf < cosASVf);
	}

	float solventDistance(Point3f ptV, Point3f ptA, Point3f ptB, float rAS, float rBS, float dAB, float dAV, float dBV) {
		double angleVAB = Math.acos((dAV * dAV + dAB * dAB - dBV * dBV) / (2 * dAV * dAB));
		double angleBAS = Math.acos((dAB * dAB + rAS * rAS - rBS * rBS) / (2 * dAB * rAS));
		float dVS = (float) Math.sqrt(rAS * rAS + dAV * dAV - 2 * rAS * dAV * Math.cos(angleBAS - angleVAB));
		return dVS;
	}

	// ////// file-based data already in Hashtable /////////

	float[][][] tempVoxelData;

	void setupSurfaceInfo() {
		volumetricOrigin.set((Point3f) surfaceInfo.get("volumetricOrigin"));
		Vector3f[] v = (Vector3f[]) surfaceInfo.get("volumetricVectors");
		for (int i = 0; i < 3; i++) {
			volumetricVectors[i].set(v[i]);
			volumetricVectorLengths[i] = volumetricVectors[i].length();
			unitVolumetricVectors[i].normalize(volumetricVectors[i]);
		}
		int[] counts = (int[]) surfaceInfo.get("voxelCounts");
		for (int i = 0; i < 3; i++)
			voxelCounts[i] = counts[i];
		tempVoxelData = voxelData = (float[][][]) surfaceInfo.get("voxelData");
		precalculateVoxelData = true;
	}

	// ////// user function //////////

	String functionName;

	void setupFunctionXY() {
		jvxlFileHeader = new StringBuffer();
		jvxlFileHeader.append("functionXY\n" + functionXYinfo + "\n");
		functionName = (String) functionXYinfo.get(0);
		volumetricOrigin.set((Point3f) functionXYinfo.get(1));
		if (!isAngstroms)
			volumetricOrigin.scale(ANGSTROMS_PER_BOHR);
		for (int i = 0; i < 3; i++) {
			Point4f info = (Point4f) functionXYinfo.get(i + 2);
			voxelCounts[i] = (int) info.x;
			volumetricVectors[i].set(info.y, info.z, info.w);
			if (!isAngstroms)
				volumetricVectors[i].scale(ANGSTROMS_PER_BOHR);
			volumetricVectorLengths[i] = volumetricVectors[i].length();
			unitVolumetricVectors[i].normalize(volumetricVectors[i]);
		}
		jvxlFileHeader.append(jvxlGetVolumeHeader(2));
		atomCount = 0;
		negativeAtomCount = false;
	}

	float getFunctionValue(int x, int y) {
		return viewer.functionXY(functionName, x, y);
	}

	// // LCAO Cartoons ////

	int nLCAO = 0;

	void drawLcaoCartoon(String lcaoCartoon, Vector3f z, Vector3f x) {
		Vector3f y = new Vector3f();
		boolean isReverse = (lcaoCartoon.length() > 0 && lcaoCartoon.charAt(0) == '-');
		if (isReverse)
			lcaoCartoon = lcaoCartoon.substring(1);
		colorPos = colorPosLCAO;
		colorNeg = colorNegLCAO;
		int sense = (isReverse ? -1 : 1);
		y.cross(z, x);
		String id = (currentMesh == null ? "lcao" + (++nLCAO) + "_" + lcaoCartoon : currentMesh.thisID);
		if (currentMesh == null)
			allocMesh(id);
		defaultColix = Graphics3D.getColix(colorPos);

		if (lcaoCartoon.equals("px")) {
			currentMesh.thisID += "a";
			createLcaoLobe(x, sense);
			setProperty("thisID", id + "b", null);
			createLcaoLobe(x, -sense);
			currentMesh.colix = Graphics3D.getColix(colorNeg);
			return;
		}
		if (lcaoCartoon.equals("py")) {
			currentMesh.thisID += "a";
			createLcaoLobe(y, sense);
			setProperty("thisID", id + "b", null);
			createLcaoLobe(y, -sense);
			currentMesh.colix = Graphics3D.getColix(colorNeg);
			return;
		}
		if (lcaoCartoon.equals("pz")) {
			currentMesh.thisID += "a";
			createLcaoLobe(z, sense);
			setProperty("thisID", id + "b", null);
			createLcaoLobe(z, -sense);
			currentMesh.colix = Graphics3D.getColix(colorNeg);
			return;
		}
		if (lcaoCartoon.equals("pxa")) {
			createLcaoLobe(x, sense);
			return;
		}
		if (lcaoCartoon.equals("pxb")) {
			createLcaoLobe(x, -sense);
			return;
		}
		if (lcaoCartoon.equals("pya")) {
			createLcaoLobe(y, sense);
			return;
		}
		if (lcaoCartoon.equals("pyb")) {
			createLcaoLobe(y, -sense);
			return;
		}
		if (lcaoCartoon.equals("pza")) {
			createLcaoLobe(z, sense);
			return;
		}
		if (lcaoCartoon.equals("pzb")) {
			createLcaoLobe(z, -sense);
			return;
		}
		if (lcaoCartoon.indexOf("sp") == 0 || lcaoCartoon.indexOf("lp") == 0) {
			createLcaoLobe(z, sense);
			return;
		}

		// assume s
		createLcaoLobe(null, 1);
		return;
	}

	Point4f lcaoDir = new Point4f();

	void createLcaoLobe(Vector3f lobeAxis, float factor) {
		initState();
		Logger.debug("creating isosurface " + currentMesh.thisID);
		if (lobeAxis == null) {
			setProperty("sphere", new Float(factor / 2f), null);
			return;
		}
		lcaoDir.x = lobeAxis.x * factor;
		lcaoDir.y = lobeAxis.y * factor;
		lcaoDir.z = lobeAxis.z * factor;
		lcaoDir.w = 0.7f;
		setProperty("lobe", lcaoDir, null);
	}

	void setModelIndex() {
		setModelIndex(atomIndex);
		currentMesh.ptCenter.set(center);
		currentMesh.title = title;
		currentMesh.jvxlDefinitionLine = jvxlGetDefinitionLine(currentMesh);
		if (script != null)
			currentMesh.scriptCommand = fixScript();
	}

	@SuppressWarnings("unchecked")
	Vector getShapeDetail() {
		Vector V = new Vector();
		for (int i = 0; i < meshCount; i++) {
			Hashtable info = new Hashtable();
			Mesh mesh = meshes[i];
			if (mesh == null)
				continue;
			info.put("ID", (mesh.thisID == null ? "<noid>" : mesh.thisID));
			info.put("vertexCount", new Integer(mesh.vertexCount));
			if (mesh.ptCenter.x != Float.MAX_VALUE)
				info.put("center", mesh.ptCenter);
			if (mesh.jvxlDefinitionLine != null)
				info.put("jvxlDefinitionLine", mesh.jvxlDefinitionLine);
			info.put("modelIndex", new Integer(mesh.modelIndex));
			if (mesh.title != null)
				info.put("title", mesh.title);
			V.add(info);
		}
		return V;
	}

	BitSet[] surfaceSet;
	int nSets = 0;
	boolean setsSuccessful;

	BitSet[] getSurfaceSet(int level) {
		if (currentMesh == null)
			return null;
		if (level == 0) {
			surfaceSet = new BitSet[100];
			nSets = 0;
		}
		setsSuccessful = true;
		for (int i = 0; i < currentMesh.polygonCount; i++) {
			int[] p = currentMesh.polygonIndexes[i];
			int pt0 = findSet(p[0]);
			int pt1 = findSet(p[1]);
			int pt2 = findSet(p[2]);
			if (pt0 < 0 && pt1 < 0 && pt2 < 0) {
				createSet(p[0], p[1], p[2]);
				continue;
			}
			if (pt0 == pt1)
				continue;
			if (pt0 >= 0) {
				surfaceSet[pt0].set(p[1]);
				surfaceSet[pt0].set(p[2]);
				if (pt1 >= 0 && pt1 != pt0)
					mergeSets(pt0, pt1);
				if (pt2 >= 0 && pt2 != pt0 && pt2 != pt1)
					mergeSets(pt0, pt2);
				continue;
			}
			if (pt1 >= 0) {
				surfaceSet[pt1].set(p[0]);
				surfaceSet[pt1].set(p[2]);
				if (pt2 >= 0 && pt2 != pt1)
					mergeSets(pt1, pt2);
				continue;
			}
			surfaceSet[pt2].set(p[0]);
			surfaceSet[pt2].set(p[1]);
		}
		int n = 0;
		for (int i = 0; i < nSets; i++)
			if (surfaceSet[i] != null)
				n++;
		BitSet[] temp = new BitSet[n];
		n = 0;
		for (int i = 0; i < nSets; i++)
			if (surfaceSet[i] != null)
				temp[n++] = surfaceSet[i];
		nSets = n;
		surfaceSet = temp;
		if (!setsSuccessful && level < 2)
			getSurfaceSet(++level);
		return surfaceSet;
	}

	int findSet(int vertex) {
		for (int i = 0; i < nSets; i++)
			if (surfaceSet[i] != null && surfaceSet[i].get(vertex))
				return i;
		return -1;
	}

	void createSet(int v1, int v2, int v3) {
		int i;
		for (i = 0; i < nSets; i++)
			if (surfaceSet[i] == null)
				break;
		if (i >= 100) {
			setsSuccessful = false;
			return;
		}
		if (i == nSets)
			nSets = i + 1;
		surfaceSet[i] = new BitSet();
		surfaceSet[i].set(v1);
		surfaceSet[i].set(v2);
		surfaceSet[i].set(v3);
	}