class CPSADescriptor extends AbstractMolecularDescriptor implements IMolecularDescriptor {

    @Override
    public DescriptorValue calculate(IAtomContainer atomContainer) {
        DoubleArrayResult retval = new DoubleArrayResult();

        if (!GeometryUtil.has3DCoordinates(atomContainer)) {
            for (int i = 0; i < 29; i++)
                retval.add(Double.NaN);
            return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                    getDescriptorNames(), new CDKException("Molecule must have 3D coordinates"));
        }

        IAtomContainer container;
        try {
            container = (IAtomContainer) atomContainer.clone();
        } catch (CloneNotSupportedException e) {
            logger.debug("Error during clone");
            for (int i = 0; i < 29; i++)
                retval.add(Double.NaN);
            return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                    getDescriptorNames(), new CDKException("Error during clone" + e.getMessage()));
        }

        //        IsotopeFactory factory = null;
        //        try {
        //            factory = IsotopeFactory.getInstance(container.getNewBuilder());
        //        } catch (Exception e) {
        //            logger.debug(e);
        //        }

        GasteigerMarsiliPartialCharges peoe;
        try {
            peoe = new GasteigerMarsiliPartialCharges();
            peoe.assignGasteigerMarsiliSigmaPartialCharges(container, true);
        } catch (Exception e) {
            logger.debug("Error in assigning Gasteiger-Marsilli charges");
            for (int i = 0; i < 29; i++)
                retval.add(Double.NaN);
            return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                    getDescriptorNames(), new CDKException("Error in getting G-M charges"));
        }

        NumericalSurface surface;
        try {
            surface = new NumericalSurface(container);
            surface.calculateSurface();
        } catch (NullPointerException npe) {
            logger.debug("Error in surface area calculation");
            for (int i = 0; i < 29; i++)
                retval.add(Double.NaN);
            return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                    getDescriptorNames(), new CDKException("Error in surface area calculation"));
        }

        //double molecularWeight = mfa.getMass();
        double[] atomSurfaces = surface.getAllSurfaceAreas();
        double totalSA = surface.getTotalSurfaceArea();

        double ppsa1 = 0.0;
        double ppsa3 = 0.0;
        double pnsa1 = 0.0;
        double pnsa3 = 0.0;
        double totpcharge = 0.0;
        double totncharge = 0.0;
        for (int i = 0; i < container.getAtomCount(); i++) {
            if (container.getAtom(i).getCharge() > 0) {
                ppsa1 += atomSurfaces[i];
                ppsa3 += container.getAtom(i).getCharge() * atomSurfaces[i];
                totpcharge += container.getAtom(i).getCharge();
            } else {
                pnsa1 += atomSurfaces[i];
                pnsa3 += container.getAtom(i).getCharge() * atomSurfaces[i];
                totncharge += container.getAtom(i).getCharge();
            }
        }

        double ppsa2 = ppsa1 * totpcharge;
        double pnsa2 = pnsa1 * totncharge;

        // fractional +ve & -ve SA
        double fpsa1 = ppsa1 / totalSA;
        double fpsa2 = ppsa2 / totalSA;
        double fpsa3 = ppsa3 / totalSA;
        double fnsa1 = pnsa1 / totalSA;
        double fnsa2 = pnsa2 / totalSA;
        double fnsa3 = pnsa3 / totalSA;

        // surface wtd +ve & -ve SA
        double wpsa1 = ppsa1 * totalSA / 1000;
        double wpsa2 = ppsa2 * totalSA / 1000;
        double wpsa3 = ppsa3 * totalSA / 1000;
        double wnsa1 = pnsa1 * totalSA / 1000;
        double wnsa2 = pnsa2 * totalSA / 1000;
        double wnsa3 = pnsa3 * totalSA / 1000;

        // hydrophobic and poalr surface area
        double phobic = 0.0;
        double polar = 0.0;
        for (int i = 0; i < container.getAtomCount(); i++) {
            if (Math.abs(container.getAtom(i).getCharge()) < 0.2) {
                phobic += atomSurfaces[i];
            } else {
                polar += atomSurfaces[i];
            }
        }
        double thsa = phobic;
        double tpsa = polar;
        double rhsa = phobic / totalSA;
        double rpsa = polar / totalSA;

        // differential +ve & -ve SA
        double dpsa1 = ppsa1 - pnsa1;
        double dpsa2 = ppsa2 - pnsa2;
        double dpsa3 = ppsa3 - pnsa3;

        double maxpcharge = 0.0;
        double maxncharge = 0.0;
        int pidx = 0;
        int nidx = 0;
        for (int i = 0; i < container.getAtomCount(); i++) {
            double charge = container.getAtom(i).getCharge();
            if (charge > maxpcharge) {
                maxpcharge = charge;
                pidx = i;
            }
            if (charge < maxncharge) {
                maxncharge = charge;
                nidx = i;
            }
        }

        // relative descriptors
        double rpcg = maxpcharge / totpcharge;
        double rncg = maxncharge / totncharge;
        double rpcs = atomSurfaces[pidx] * rpcg;
        double rncs = atomSurfaces[nidx] * rncg;

        // fill in the values
        retval.add(ppsa1);
        retval.add(ppsa2);
        retval.add(ppsa3);
        retval.add(pnsa1);
        retval.add(pnsa2);
        retval.add(pnsa3);

        retval.add(dpsa1);
        retval.add(dpsa2);
        retval.add(dpsa3);

        retval.add(fpsa1);
        retval.add(fpsa2);
        retval.add(fpsa3);
        retval.add(fnsa1);
        retval.add(fnsa2);
        retval.add(fnsa3);

        retval.add(wpsa1);
        retval.add(wpsa2);
        retval.add(wpsa3);
        retval.add(wnsa1);
        retval.add(wnsa2);
        retval.add(wnsa3);

        retval.add(rpcg);
        retval.add(rncg);
        retval.add(rpcs);
        retval.add(rncs);

        retval.add(thsa);
        retval.add(tpsa);
        retval.add(rhsa);
        retval.add(rpsa);

        return new DescriptorValue(getSpecification(), getParameterNames(), getParameters(), retval,
                getDescriptorNames());
    }

    public  CPSADescriptor();

}