package sansmodels;

public class SANSResolution extends Object {
    private double lambda, lambdaWidth; 
    private double L1, L2, S1, S2, BS, DDet, apOff;

    // default constructor
    public SANSResolution() {
        lambda = 6.00;
        lambdaWidth = 0.10;
        L1 = 16.32;
        L2 = 11.0;
        S1 = 50.0;
        S2 = 12.7;
        BS = 50.8;
        DDet = 1.0;
        apOff = 0.0;
    }

    public SANSResolution( double inLam, double inLamW,
                           double inL1, double inL2,
                           double inS1, double inS2,
                           double inBS, double inDDet,
                           double inApOff) {
        lambda = inLam;
        lambdaWidth = inLamW;
        L1 = inL1;
        L2 = inL2;
        S1 = inS1;
        S2 = inS2;
        BS = inBS;
        DDet = inDDet;
        apOff = inApOff;
    }

    //returns Sigma-Q, Q-bar, and shadow corrections
    public double[] getResolution(double inQ) {

        //lots of calculation variables
        double a2, q_small, lp, v_lambda, v_b, v_d, vz, yg, v_g;
        double r0, delta, inc_gamma, fr, fv, rmd, v_r1, rm, v_r;

        //Output Quantities
        double fSubS, QBar, SigmaQ;
        double outResolution[] = new double[3];

        //Constants
        double del_r = .1;
        double vz_1 = 3.956e5;
        double g = 981.0;

        
        //rename for working variables
        double wLam, wLW, wL1, wL2, wS1, wS2;
        double wBS, wDDet, wApOff;
        wLam = lambda; 
        wLW = lambdaWidth;
        wDDet = DDet;
        wApOff = apOff;
        wS1 = 0.5*0.1*S1;
        wS2 = 0.5*0.1*S2;

        wL1 = 100.0*L1;
        wL1 -= apOff;

        wL2 = 100.0*L2;
        wL2 += apOff;

        wBS = 0.5*0.1*BS;
        
        //Start resolution calculation
        a2 = wS1*wL2/wL1 + wS2*(wL1+wL2)/wL1;
        q_small = 2.0*Math.PI*(wBS-a2)*(1.0-wLW)/(wLam*wL2);
        lp = 1.0/( 1.0/wL1 + 1.0/wL2);

        v_lambda = wLW*wLW/6.0;
        v_b = 0.25*(wS1*wS1*wL2*wL2/wL1/wL1) 
              +0.25*(wS2*wS2*wL2*wL2/lp/lp);
        v_d = (wDDet*wDDet/2.3548/2.3548) + del_r*del_r/12.0;
        vz = vz_1 / wLam;
        yg = 0.5*g*wL2*(wL1+wL2)/vz/vz;
        v_g = 2.0*yg*yg*v_lambda;
        
        r0 = wL2*Math.tan(2.0*arcsin(wLam*inQ/4.0/Math.PI));
        delta = 0.5*(wBS - r0)*(wBS - r0)/v_d;
        if (r0 < wBS) {
            inc_gamma=Math.exp(gammln(1.5))*(1-gammp(1.5,delta));
        } else {
            inc_gamma=Math.exp(gammln(1.5))*(1+gammp(1.5,delta));
        }

        fSubS = 0.5*(1.0+erf( (r0-wBS)/Math.sqrt(2.0*v_d) ) );
        if (fSubS <= 0.0) fSubS = 1.e-10;
        fr = 1.0 + Math.sqrt(v_d)*Math.exp(-1.0*delta)
             /(r0*fSubS*Math.sqrt(2.0*Math.PI));
        fv = inc_gamma/(fSubS*Math.sqrt(Math.PI)) 
             - r0*r0*(fr-1.0)*(fr-1.0)/v_d;

        rmd = fr*r0;
        v_r1 = v_b + fv*v_d +v_g;

        rm = rmd + 0.5*v_r1/rmd;
        v_r = v_r1 - 0.5*(v_r1*v_r1/rmd/rmd);
        if (v_r < 0.0) v_r = 0.0;
        QBar = (4.0*Math.PI/wLam)*Math.sin(0.5*Math.atan(rm/wL2));
        SigmaQ = QBar*Math.sqrt(v_r/rmd/rmd +v_lambda);

        outResolution[0] = SigmaQ;
        outResolution[1] = QBar;
        outResolution[2] = fSubS;
        return outResolution;
    }

    private double gammln(double inX) {
        double stp, half, one, fpf, x, tmp, ser;
        double[] cof = { 76.18009173,
                        -86.50532033,
                         24.01409822,
                         -1.231739516,
                           .120858003e-2,
                          -.536382e-5};
        stp = 2.50662827465; 
        half = 0.5; one = 1.0; fpf = 5.5;

        x = inX-one;
        tmp = x+fpf;
        tmp = (x+half)*Math.log(tmp)-tmp;
        ser=one;
        for(int i=0;i<6;i++) {
           x += one;
           ser += cof[i]/x;
        }
        return tmp + Math.log(stp*ser);
    }

    private double gser( double a, double x) {
        if(x <= 0.0) return 0.0;
        
        double ap, sum, del;
        int itmax = 100;
        double eps = 3.0e-7;
        int i=0;

        ap = a;
        sum = 1.0/a;
        del = sum;
       
        while(i<itmax) {
            ap++;
            del *= x/ap;
            sum += del;
            i++;
            if(Math.abs(del) < Math.abs(sum)*eps) i=itmax;
        }
        return sum*Math.exp(-1.0*x + a*Math.log(x) - gammln(a));
    }

    private double gcf( double a, double x) {

        double gold, a0, a1, b0, b1, fac, an, ana, anf, g;
        int itmax = 100;
        double eps = 3.0e-7;
        int i = 1;

        gold = 0.0;
        a0 = 1.0;
        a1 = x;
        b0 = 0.0;
        b1 = 1.0;
        fac = 1.0;
        g = 1.0;
 
        while(i<itmax) {
            an = (double)i;
            ana = an - a;
            a0 = (a1+ a0*ana)*fac;
            b0 = (b1+ b0*ana)*fac;
            anf = an*fac;
            a1 = x*a0 + anf*a1;
            b1 = x*b0 + anf*b1;
            if(a1 != 0.0) {
                fac = 1.0/a1;
                g=b1*fac;
                if(Math.abs((g-gold)/g) < eps) i=itmax;
                gold = g;
            }
        }
        return g*Math.exp(-1.0*x +a*Math.log(x) -gammln(a));
    }

    private double gammp( double a, double x) {
        if ((x < 0.0) || (a <= 0.0)) return 0.0;
        if (x < a+1.0) {
           return gser(a,x);
        } else {
           return 1.0-gcf(a,x);
        }
    }

    private double erf( double inX) {
        if (inX < 0.0) {
           return -1.0*gammp(0.5,inX*inX);
        } else {
           return gammp(0.5,inX*inX);
        }
    }
    
    private double arcsin( double inX) {
        if((1.0-inX) < 1.0e-6) {
            return Math.PI/2.0;
        }
        return Math.asin(inX);
    }
    
    private double gauss(double q, double qexp, double sigma) {
        return Math.exp(-0.5*(q-qexp)*(q-qexp)/sigma/sigma)
               /2.4755/sigma;
    }
    
    public double[] getParameters() {
        double outParameters[] = new double[9];
        
        outParameters[0]=lambda;
        outParameters[1]=lambdaWidth;
        outParameters[2]=L1;
        outParameters[3]=L2;
        outParameters[4]=S1;
        outParameters[5]=S2;
        outParameters[6]=BS;
        outParameters[7]=DDet;
        outParameters[8]=apOff;
        return outParameters;
    }

    public void setParameters(double[] inParameters) {
        if(inParameters.length != 9) return;
        lambda = inParameters[0];
        lambdaWidth = inParameters[1];
        L1 = inParameters[2];
        L2 = inParameters[3];
        S1 = inParameters[4];
        S2 = inParameters[5];
        BS = inParameters[6];
        DDet = inParameters[7];
        apOff = inParameters[8];
    }

}