package sansmodels;

public class PolyHardSphere extends SANSModel{
 
    private double polyDisp, radius, volumeFraction;
    private double deltaRho, background;

    //Default constructor
    public PolyHardSphere() {
        radius = 100.0;
        polyDisp = 0.12;
        volumeFraction = 0.1;
        deltaRho = 2.e-6;
        background = 0.0;
        setNumberOfParameters(5);
        String[] parameters = {"Radius (A)",
            "Polydispersity (0-1)", "Volume Fraction (0-1)",
            "Contrast (A-2)", "Background (cm-1)"};
        setParametersText(parameters);
    }

    public PolyHardSphere(double inRadius, 
                          double inPolyDisp, double inVolFrac,
                          double inDeltaRho, double inBackground) {
        radius = inRadius; 
        polyDisp = inPolyDisp;
        volumeFraction = inVolFrac;
        if (volumeFraction >= 1.0) volumeFraction = 0.99999;
        deltaRho = inDeltaRho;
        background = inBackground;
        setNumberOfParameters(5);
        String[] parameters = {"Radius (A)",
            "Polydispersity (0-1)", "Volume Fraction (0-1)",
            "Contrast (A-2)", "Background (cm-1)"};
        setParametersText(parameters);
    }

    public double getFormFactor(double inX) {
        double mu,mu1,d1,d2,d3,d4,d5,d6,capd,rho,beta;
        double l,l1,b,c,chi,chi1,chi2,e,t,t1,t2,t3,p;
        double ka,zz,v1,v2,p1,p2;
        double h1,h2,h3,h4,e1,y,y1,s,s1,s2,s3,hint1,hint2;
        double capl,capl1,capmu,capmu1,r3,pq,happ;
        double ka2,r1,heff,k;
        double h, sigma, cont, rad, z1, bkg, phi, z2;
 
        if (polyDisp <= 0.0)polyDisp = 0.0001;
        if (polyDisp >= 1.0)polyDisp = 1.0;
        k = inX;
        rad = radius;
        sigma = 2.0*radius;
        z2 = polyDisp;
        phi = volumeFraction;
        cont = deltaRho*1.0e4;
        bkg = background;

        zz=1/(z2*z2)-1.0;
        b = sigma/(zz+1);
        c = zz+1;
 
        r1 = sigma/2.0;
        r3 = r1*r1*r1;
        r3 *= (zz+2)*(zz+3)/((zz+1)*(zz+1));
        rho=phi/(1.3333333333*Math.PI*r3);
        t1 = rho*b*c;
        t2 = rho*b*b*c*(c+1);
        t3 = rho*b*b*b*c*(c+1)*(c+2);
        capd = 1-Math.PI*t3/6;
 
        v1=1/(1+b*b*k*k);
        v2=1/(4+b*b*k*k);
        p=Math.pow(v1,(c/2))*Math.sin(c*Math.atan(b*k));
        p1=b*c*Math.pow(v1,((c+1)/2))*Math.sin((c+1)*Math.atan(b*k));
        p2=c*(c+1)*b*b*Math.pow(v1,((c+2)/2))*Math.sin((c+2)*Math.atan(b*k));
        mu=Math.pow(2,c)*Math.pow(v2,(c/2))*Math.sin(c*Math.atan(b*k/2));
        mu1=Math.pow(2,(c+1))*b*c*Math.pow(v2,((c+1)/2))*Math.sin((c+1)*Math.atan(k*b/2));
        s1=b*c;
        s2=c*(c+1)*b*b;
        s3=c*(c+1)*(c+2)*b*b*b;
        chi=Math.pow(v1,(c/2))*Math.cos(c*Math.atan(b*k));
        chi1=b*c*Math.pow(v1,((c+1)/2))*Math.cos((c+1)*Math.atan(b*k));
        chi2=c*(c+1)*b*b*Math.pow(v1,((c+2)/2))*Math.cos((c+2)*Math.atan(b*k));
        l=Math.pow(2,c)*Math.pow(v2,(c/2))*Math.cos(c*Math.atan(b*k/2));
        l1=Math.pow(2,(c+1))*b*c*Math.pow(v2,((c+1)/2))*Math.cos((c+1)*Math.atan(k*b/2));
        d1=(Math.PI/capd)*(2+(Math.PI/capd)*(t3-(rho/k)*(k*s3-p2)));
        d2=Math.pow((Math.PI/capd),2.0)*(rho/k)*(k*s2-p1);
        d3=(-1.0)*Math.pow((Math.PI/capd),2.0)*(rho/k)*(k*s1-p);
        d4=(Math.PI/capd)*(k-(Math.PI/capd)*(rho/k)*(chi1-s1));
        d5=Math.pow((Math.PI/capd),2.0)*((rho/k)*(chi-1)+0.5*k*t2);
        d6=Math.pow((Math.PI/capd),2.0)*(rho/k)*(chi2-s2);
        e1=Math.pow((Math.PI/capd),2.0)*Math.pow((rho/k/k),2.0)
           *((chi-1)*(chi2-s2)
           -(chi1-s1)*(chi1-s1)-(k*s1-p)*(k*s3-p2)+Math.pow((k*s2-p1),2.0));
        e=1-(2*Math.PI/capd)*(1+0.5*Math.PI*t3/capd)*(rho/k/k/k)*(k*s1-p)
           -(2*Math.PI/capd)*rho/k/k*((chi1-s1)
           +(0.25*Math.PI*t2/capd)*(chi2-s2))-e1;
        y1=Math.pow((Math.PI/capd),2.0)*Math.pow((rho/k/k),2.0)
           *((k*s1-p)*(chi2-s2)
           -2*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1));
        y = (2*Math.PI/capd)*(1+0.5*Math.PI*t3/capd)*(rho/k/k/k)
           *(chi+0.5*k*k*s2-1)-(2*Math.PI*rho/capd/k/k)*(k*s2-p1
           +(0.25*Math.PI*t2/capd)*(k*s3-p2))-y1;
        capl=2.0*Math.PI*cont*rho/k/k/k*(p-0.5*k*(s1+chi1));
        capl1=2.0*Math.PI*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
        capmu=2.0*Math.PI*cont*rho/k/k/k*(1-chi-0.5*k*p1);
        capmu1=2.0*Math.PI*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
 
        h1=capl*(capl*(y*d1-e*d6)+capl1*(y*d2-e*d4)
           +capmu*(e*d1+y*d6)+capmu1*(e*d2+y*d4));
        h2=capl1*(capl*(y*d2-e*d4)+capl1*(y*d3-e*d5)
           +capmu*(e*d2+y*d4)+capmu1*(e*d3+y*d5));
        h3=capmu*(capl*(e*d1+y*d6)+capl1*(e*d2+y*d4)
           +capmu*(e*d6-y*d1)+capmu1*(e*d4-y*d2));
        h4=capmu1*(capl*(e*d2+y*d4)+capl1*(e*d3+y*d5)
           +capmu*(e*d4-y*d2)+capmu1*(e*d5-y*d3));

        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(e*e+y*y));
 
        pq=4*Math.PI*cont*(Math.sin(k*sigma/2)-0.5*k*sigma*Math.cos(k*sigma/2));
        hint2=8*Math.PI*Math.PI*rho*cont*cont/(k*k*k*k*k*k)
              *(1-chi-k*p1+0.25*k*k*(s2+chi2));
 
        ka=k*(sigma/2);
        h=hint1+hint2;
        heff=1.0+hint1/hint2;
        ka2=ka*ka;
 
        return (h+bkg);
    }


    public void setParameters(double[] inParameters) {
        radius = inParameters[0];
        polyDisp = inParameters[1];
        volumeFraction = inParameters[2];
        if (volumeFraction >= 1.0) volumeFraction = 0.99999;
        deltaRho = inParameters[3];
        background = inParameters[4];
    }
 
    public double[] getParameters() {
        double[] outParameters = new double[5];
        outParameters[0] = radius;
        outParameters[1] = polyDisp;
        outParameters[2] = volumeFraction;
        outParameters[3] = deltaRho;
        outParameters[4] = background;
        return outParameters;
    }
}