/**
 *  A class of utilities that are useful for multiple programs and so 
 *  get put in their own class.
 */
import java.util.Vector;
import org.apache.oro.text.regex.*;

public final class NCNRUtils extends Object  {

    /**
     *  Returns a number really rounded to the correct number of 
     *  significant figures.  Otherwise if double precision arithmetic
     *  is used on numbers, a number in the least precision bit can
     *  completely throw off the formatting of that number as a string.
     *  This is only necessary because the Java 1.2 formatting classes
     *  cannot be used at this time with all browsers.
     */
    public static double sigfigs(double in, int inSigFigs) {
        if (in == 0.0) return in;
        inSigFigs = inSigFigs-1;
        if (inSigFigs < 1) inSigFigs = 1;
        double tempValue;
        double count;
        long tempSLD;
        tempValue = in;
        count = 1.0;
        double maxValue=1.0;
        for(int i=0;i<inSigFigs;i++)
            maxValue *=10.0;
        if (in > maxValue) {
            while(Math.abs(tempValue) > maxValue) {
                tempValue /= 10.0;
                count /= 10.0;
            }
        } else {
            while(Math.abs(tempValue) < maxValue) {
                tempValue *= 10.0;
                count *= 10.0;
            }
        }
        tempSLD = Math.round(tempValue);
        return (double)tempSLD/count;
    }

    /**
     *  Takes a string which specifies a chemical compound complete with
     *  isotopes in parenthesis and constructs a vector of elements and
     *  numbers of each element in each unit cell.
     */
    public static Vector parseCompound(String inputString) throws 
        BadCompoundException {

        //parse chemical field
        Element thisElement;
        Vector elementVector = new Vector();
        PatternMatcher matcher, matcher1, matcher2, matcher3;
        PatternCompiler compiler;
        Pattern pattern, pattern1, pattern2, pattern3;
        PatternMatcherInput input, input1, input2, input3;
        MatchResult result, result1, result2, result3;
        String elementString, nameString, isotopeString, repeatString;
        compiler = new Perl5Compiler();
        matcher = new Perl5Matcher();
        matcher1 = new Perl5Matcher();
        matcher2 = new Perl5Matcher();
        matcher3 = new Perl5Matcher();
        try{
            pattern = compiler.compile("[A-Z][^A-Z]*");
            pattern1 = compiler.compile("[^\\d|\\(]+");
            pattern2 = compiler.compile("\\((\\d+)\\)");
            pattern3 = compiler.compile("\\d+$");
        } catch(MalformedPatternException exc) {
            System.out.println("Bad pattern.");
            System.out.println(exc.getMessage());
            throw new BadCompoundException("Bad pattern in NCNRUtils.");
        }
        input = new PatternMatcherInput(inputString);
        int repeat = 0;
        int isotope = 0;
        int element = 0;
        int i = 0;
        nameString = "";
        while(matcher.contains(input, pattern)) {
            repeat = 0;
            isotope = 0;
            result = matcher.getMatch();
            elementString = result.toString();
            //System.out.println("ele "+elementString);
            input1 = new PatternMatcherInput(elementString);
            if(matcher1.contains(input1, pattern1)) {
                result1 = matcher1.getMatch();
                nameString = result1.toString();
                //System.out.println("name "+nameString);
                element = getElement(nameString);
                if(element == -1) { 
                    throw new BadCompoundException("Bad Element: "
                                  +nameString);
                } 
            }
            input2 = new PatternMatcherInput(elementString);
            if(matcher2.contains(input2, pattern2)) {
                result2 = matcher2.getMatch();
                isotopeString = result2.group(1);
                //System.out.println("iso "+isotopeString);
                isotope = Integer.parseInt(isotopeString);
            }
            input3 = new PatternMatcherInput(elementString);
            if(matcher3.contains(input3, pattern3)) {
                result3 = matcher3.getMatch();
                repeatString = result3.toString();
                //System.out.println("rep "+repeatString);
                repeat = Integer.parseInt(repeatString);
            }
            if(repeat == 0) repeat = 1;
            thisElement = new Element(element, isotope);
            if(thisElement.isBad()) {
                throw new BadCompoundException("No data for isotope "
                    +nameString+"("+isotope+")");
            }
            elementVector.addElement(thisElement);
            elementVector.addElement(new Integer(repeat));
            i++;
        }

        return elementVector;
    }

    /**
     *  Internal method to this class.  Matches an element string with
     *  its number in the set of element data.
     */
    private static int getElement(String in) {
        int outElement = -1;
        int i = 0;
        while((outElement == -1)&&(i < ElementData.elements.length)) {
           if(in.equalsIgnoreCase(ElementData.elements[i])) outElement = i+1;
           i++;
        }
        return outElement;
    }

    /**
     *  Actually calculates more than just the SLD's of a set of a compound
     *  it also calculates some cross sections.
     *  The nine output parameters included in the output array are
     *  output[0] = Real neutron SLD
     *  output[1] = Imaginary neutron SLD
     *  output[2] = Real Cu Kalpha SLD
     *  output[3] = Imaginary Cu Kalpha SLD
     *  output[4] = Real Mo Kalpha SLD
     *  output[5] = Imaginary Mo Kalpha SLD
     *  output[6] = Spin Incoherent neutron macroscopic cross section
     *  output[7] = Neutron macroscopic absorption cross section
     *  output[8] = Neutron 1/e length
     *  output[9] = Alloy Incoherent neutron macroscopic cross section
     */
    public static double[] calculateSLD(Vector elementVector, 
                                        double density,
                                        double wavelength) {
        int numElements = elementVector.size()/2;
        int[] repeats = new int[numElements];
        int repeat = 0;
        for(int i=0;i<numElements;i++) {
            repeat = i + 1;
            repeats[i] = 
                ((Integer)elementVector.elementAt(repeat)).intValue();
            elementVector.removeElementAt(repeat);
        }
        int numAtoms = 0;
        double cellMass = 0.0;
        double cellRealB = 0.0;
        double cellImB = 0.0;
        double cellRealXCu = 0.0;
        double cellImXCu = 0.0;
        double cellRealXMo = 0.0;
        double cellImXMo = 0.0;
        double cellIncohXS = 0.0;
        double cellScatXS = 0.0;
        double cellIncohSpinXS = 0.0;
        double cellIncohNucXS = 0.0;
        double cellAbsorbXS = 0.0;
        double output[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
        for(int i=0;i<numElements;i++) {
            numAtoms +=repeats[i];
        }
        for(int i=0;i<numElements;i++) {
            cellMass += ((Element)elementVector.elementAt(i)).mass
                        * (double)repeats[i];
            cellRealB += ((Element)elementVector.elementAt(i)).realCohB
                        * (double)repeats[i];
            cellImB += ((Element)elementVector.elementAt(i)).imCohB
                        * (double)repeats[i];
            cellRealXCu += ( (double)((Element)elementVector.elementAt(i)).Z
                          +((Element)elementVector.elementAt(i)).CuKafp)
                         * (double)repeats[i];
            cellImXCu += ((Element)elementVector.elementAt(i)).CuKafpp
                         * (double)repeats[i];
            cellRealXMo += ( (double)((Element)elementVector.elementAt(i)).Z
                          +((Element)elementVector.elementAt(i)).MoKafp)
                         * (double)repeats[i];
            cellImXMo += ((Element)elementVector.elementAt(i)).MoKafpp
                         * (double)repeats[i];
            cellScatXS += ((Element)elementVector.elementAt(i)).scatXS
                           * (double)repeats[i];
            cellAbsorbXS += ((Element)elementVector.elementAt(i)).absXS
                           * (double)repeats[i];
        }

     
        // extra numAtoms added in denominator below because of squared 
        //      term cellRealB, volume already has numAtoms in it...
        //      this is the macroscopic incoherent cross section
        // also factor of 100 in denominator because cellRealB is in 10^-15m
        //     units and cellScatXS is in 10^-24 cm units
        cellIncohNucXS = cellScatXS - 4.*Math.PI*cellRealB*cellRealB/numAtoms/100.;

        double volume = (cellMass/density)/ElementData.AVOGADRO *1.0e24;
//        System.out.println("Inc Nuc XS: "+   Double.toString(cellIncohNucXS/volume)+
//                           "\nInc Spin XS: "+Double.toString(cellIncohSpinXS/volume)+
//                           "\nCell Mass: "+Double.toString(cellMass)+
//                           "\nNumber Density: "+Double.toString(1.0/volume));

        
        System.out.println("Inc XS: "+   Double.toString(cellIncohNucXS/volume)+
                          "\nScat XS: "+ Double.toString(cellScatXS/volume)+
                          "\nAbs XS: "+ Double.toString((cellAbsorbXS/volume)*(wavelength/1.8))+
                           "\nCell Mass: "+Double.toString(cellMass)+
                           "\nNumber Density: "+Double.toString(1.0/volume));
        output[0] = cellRealB/volume*1.e-5;
        output[1] = cellImB/volume*1.e-5;
        output[2] = cellRealXCu*ElementData.RE/volume*1.0e8;
        output[3] = cellImXCu*ElementData.RE/volume*1.0e8;
        output[4] = cellRealXMo*ElementData.RE/volume*1.0e8;
        output[5] = cellImXMo*ElementData.RE/volume*1.0e8;
        output[6] = cellIncohNucXS/volume;
        output[7] = cellAbsorbXS/volume*(wavelength/1.8);
// cellScatXS = Sigma_coherent + Sigma_incoherent so dont need to add in incoh again
        output[8] = 1.0/(output[7] + cellScatXS/volume); 
        output[9] = cellIncohSpinXS/volume;
 
        return output;
    }

    /**
     *  Calculates mean free path and thickness for 90% transmission
     *  for samples for Zema.
     *  output[0] = mean free path
     *  output[1] = Thickness 
     */
    public static double[] calculateTrans(Vector elementVector, 
                                        double density,
                                        double percent,
                                        int geometry,
                                        double wavelength) {
        int numElements = elementVector.size()/2;
        int[] repeats = new int[numElements];
        int repeat = 0;
        for(int i=0;i<numElements;i++) {
            repeat = i + 1;
            repeats[i] = 
                ((Integer)elementVector.elementAt(repeat)).intValue();
            elementVector.removeElementAt(repeat);
        }
        int numAtoms = 0;
        double cellMass = 0.0;
        double cellScatXS = 0.0;
        double cellAbsorbXS = 0.0;
        double totalXS = 0.0;
        double noAbsXS = 0.0;
        double output[] = {0.0, 0.0};
        percent = percent/100.0;
        for(int i=0;i<numElements;i++) {
            numAtoms +=repeats[i];
        }
        for(int i=0;i<numElements;i++) {
            cellMass += ((Element)elementVector.elementAt(i)).mass
                        * (double)repeats[i];
            cellScatXS += ((Element)elementVector.elementAt(i)).incohXS
                           * (double)repeats[i];
            cellAbsorbXS += ((Element)elementVector.elementAt(i)).absXS
                           * (double)repeats[i];
        }
        double volume = (cellMass/density)/ElementData.AVOGADRO *1.0e24;
        cellScatXS = cellScatXS/volume;
        cellAbsorbXS = cellAbsorbXS/volume*(wavelength/1.8);
        // macroscopic cross sections now
        totalXS = (cellScatXS + cellAbsorbXS); 
        output[0] = 1.0/totalXS; // mfp in cm
        if (geometry == 0 ) {
           output[1] = -1.0/(totalXS*Math.PI)*Math.log(1.0 - (percent*totalXS/cellScatXS)); // in cm
        } else {
           output[1] = -1.0/totalXS*Math.log(1.0 - (percent*totalXS/cellScatXS)); // in cm
        }
        return output;
    }

    public static double getVolumeFraction( double massDensityA,
                                            double massDensityB,
                                            double weightFractionA ) {
        double answer;
        double c1 = massDensityA;
        double c2 = massDensityB;
        double c3 = weightFractionA;
        
        answer = c2*c3/(c2*c3 + c1*(1.0 - c3));
 
        return answer;
    }

    public static double getWeightFraction( double massDensityA,
                                            double massDensityB,
                                            double volumeFractionA ) {
        double answer;
        double c1 = massDensityA;
        double c2 = massDensityB;
        double c3 = volumeFractionA;
        
        answer = c1*c3/(c1*c3 + c2*(1.0 - c3));
        

        return answer;
    }

}