import com.imsl.math.*;
import sansmodels.Gaussian;

public class Refl5a extends ReflectivityModel {

    public Refl5a() {
    }

    public Complex getReflectivity( double inQ ) {
        Complex gj, paj, rij, arg, alpj, betj, gamj;
        Complex delj, alp, bet, gam, del, palp, pbet;
        Complex pgam, pdel, s1, s2, s3, s4, srr;
        double are, bre, cre, dre, aim, bim, cim, dim;
        double tj;
        
        double[] deltaZ = getDeltaZ();
        Complex[] rho = getRho();
        Complex fronting = getFronting();
        Complex backing = getBacking();
        int numlayers = getNumLayers();

        Complex ci = new Complex(0.0,1.0);
        Complex cr = new Complex(1.0,0.0);
        Complex kz = Complex.sqrt(Complex.plus(inQ*inQ/4.0,
            Complex.times(4.0*Math.PI,fronting)));
        Complex ribf = Complex.sqrt(Complex.minus(1.0,
            Complex.times(4.0*Math.PI,Complex.over(fronting,
            Complex.times(kz, kz)))));
        Complex riaf = Complex.sqrt(Complex.minus(1.0,
            Complex.times(4.0*Math.PI,Complex.over(backing,
            Complex.times(kz, kz)))));
        palp = new Complex(1.0);
        pbet = new Complex(0.0);
        pgam = new Complex(0.0);
        pdel = new Complex(1.0);

        alp = new Complex(0.0);
        bet = new Complex(0.0);
        gam = new Complex(0.0);
        del = new Complex(0.0);

        for (int j=0;j<numlayers;j++) {
            tj = deltaZ[j];
            gj = new Complex(rho[j]);
            paj = Complex.times(kz,tj);
            rij = Complex.sqrt(Complex.minus(1.0,
                Complex.times(4.0*Math.PI,Complex.over(gj,
                Complex.times(kz, kz)))));
            arg = Complex.times(rij,paj);
            alpj = Complex.cos(arg);
            betj = Complex.over(Complex.sin(arg),rij);
            gamj = Complex.times(Complex.negative(rij),Complex.sin(arg));
            delj = new Complex(alpj);
            alp = Complex.plus(Complex.times(alpj,palp),
                               Complex.times(betj,pgam));
            bet = Complex.plus(Complex.times(alpj,pbet),
                               Complex.times(betj,pdel));
            gam = Complex.plus(Complex.times(gamj,palp),
                               Complex.times(delj,pgam));
            del = Complex.plus(Complex.times(gamj,pbet),
                               Complex.times(delj,pdel));
            palp = new Complex(alp);
            pbet = new Complex(bet);
            pgam = new Complex(gam);
            pdel = new Complex(del);
        }
        alp = Complex.times(riaf,alp);
        bet = Complex.times(Complex.times(ribf,riaf),bet);
        gam = new Complex(gam);
        del = Complex.times(ribf,del);
        
        are = alp.real();
        bre = bet.real();
        cre = gam.real();
        dre = del.real();
        
        aim = alp.imag();
        bim = bet.imag();
        cim = gam.imag();
        dim = del.imag();
        
        s1 = Complex.plus(bet, gam);
        s2 = Complex.minus(del, alp);
        s3 = Complex.minus(bet, gam);
        s4 = Complex.plus(del, alp);
        
        srr = Complex.over(Complex.plus(s1,Complex.times(ci,s2)),
                           Complex.plus(s3,Complex.times(ci,s4)));
        return srr;
    }

    public double getSmearedReflectivity( double inQ, int inResol ) {
        int nord = 20;
        double va, vb, zi, yyy, summ, resol, intensity;
        double resolution = getResolution();
        Complex reflectivity;
        double sigma = 0.0;
        if(inResol == DELTAQ) {
            sigma = resolution;
        } else if(inResol == DELTAQOVERQ) {
            sigma = resolution*inQ;
        }
        
        va = -3.0*sigma + inQ;
        if (va < 0) va=0.0;
        vb = 3*sigma + inQ;
 
        reflectivity = new Complex(); 
        summ= 0.0;
        yyy = 0.0;

        for(int i=0;i<nord;i++) {
            zi = (Gaussian.z20[i]*(vb-va) + vb + va)/2.0;
            resol = 1.0/Math.sqrt(2.0*Math.PI*sigma*sigma);
            resol *= Math.exp((-1.0*(zi-inQ)*(zi-inQ))
                                   /(2.0*sigma*sigma));
            reflectivity = getReflectivity(zi);
            intensity = reflectivity.real()*reflectivity.real() +
                        reflectivity.imag()*reflectivity.imag();
            yyy = Gaussian.w20[i]*resol*intensity;
            summ += yyy;
        }
 
        return (vb-va)/2.0*summ;
    }

}