package gen;
import mmc.*;

public class IntFlux extends PhysicsModel{

    final private static double X0=1.148;

    private double z0=0;
    private double R=Atmosphere.R+z0;

    private int P=0;
    private int M=0;
    private int S=0;

    private double xpar[][][]={
	{  /* Standard US Atmosphere average parameters */
	    { 0.00851164, 0.124534, 0.059761, 2.32876, 19.279 },
	    { 0.102573, -0.068287, 0.958633, 0.0407253, 0.817285 },
	    { -0.017326, 0.114236, 1.15043, 0.0200854, 1.16714 },
	    { 1.3144, 50.2813, 1.33545, 0.252313, 41.0344 }
	}
    };

    private double fpar[][]={
	{ 0.701459, 2.71461, 1, 1, 1 },
	{ 0.590968, 2.69506, 1, 0, 0 },
	{ 0.594662, 2.69671, 0, 0, 0 }
    };

    private double xx=0, xs=0, xo=0, xl=0;
    private Integral I;

    public IntFlux(int model, int flag){
	I = new Integral(iromb, imaxs, iprec);
	M=model;
	S=flag;
    }

    public static void main(String[] args){
	double x;
	IntFlux F = new IntFlux(0, 0);

	switch(1){
	case 1:
	    for(int i=0; i<=1000; i++){
		x=i/1000.;
		Output.out.println(Output.f(x)+" "+Output.f(F.getIntFlux(600, x)));
	    }
	    break;
	case 2:
	    x=1;
	    Output.err.println("At cos="+Output.f(x)+", total flux is "+Output.f(F.getIntFlux(600, x)));
	    for(int i=0; i<=1000000; i++){
		Output.out.println(Output.f(F.getE(600, x, Math.random()))+" "+Output.f((F.o(x)-F.X0)));
	    }
	    break;
	default:
	}
    }

    public double getIntFlux(double E, double x){
	return getE(E, x, -1);
    }

    public double getE(double E0, double x, double rnd){
	xx=x;
	switch(S){
	case 1:
	    xs=Math.sqrt(1-x*x)/(1+h(x)/R);
	    xs=Math.sqrt(1-xs*xs);
	    break;
	case 0:
	default:
	    xs=c(x);
	}
	if(P==0){
	    xo=fpar[M][3]*o(x)-X0; if(xo<0) xo=0;
	    xl=fpar[M][4]*l(x)*1.e2*Mmu/(Lmu*C);
	}
	if(rnd<0) return I.integrateWithSubstitution(1, -E0, this, -2);
	else{
	    I.integrateWithSubstitution(1, -E0, this, -2, rnd);
	    return -I.getUpperLimit();
	}
    }

    public double h(double x){
	double s=Math.sqrt(1-x*x);
	double y=s/(1+xpar[M][0][0]);
	return xpar[M][0][4]*(1+xpar[M][0][2]*Math.pow(s, xpar[M][0][3]))/
	    Math.pow(1-(y*y), xpar[M][0][1]);
    }

    public double c(double x){
	return Math.sqrt((x*x+fpar[M][2]*
			  (xpar[M][1][0]*xpar[M][1][0]+
			   xpar[M][1][1]*Math.pow(x, xpar[M][1][2])+
			   xpar[M][1][3]*Math.pow(x, xpar[M][1][4])))/
			 (1+fpar[M][2]*(xpar[M][1][0]*xpar[M][1][0]+
					xpar[M][1][1]+xpar[M][1][3])));
    }

    public double o(double x){
	return 1/(xpar[M][2][0]+xpar[M][2][1]*Math.pow(x, xpar[M][2][2])+
		  xpar[M][2][3]*Math.pow(1-x*x, xpar[M][2][4]));
    }

    public double l(double x){
	return 1e3/(xpar[M][3][0]+xpar[M][3][1]*Math.pow(x, xpar[M][3][2])+
				  xpar[M][3][3]*Math.pow(1-x*x, xpar[M][3][4]));
    }

    public double function(double x){
	double ei, ef;
	ef=-x;
	switch(P){
	case 1:
	    {
		final double E1=121;
		final double E2=897;
		final double pK=0.213;
		final double A=2.85e-2;
		ei=ef;
		return A*fpar[M][0]*Math.pow(ei,-fpar[M][1])*(1/(1+6*ei*xs/E1)+pK/(1+1.44*ei*xs/E2));
	    }
	case 2:
	    {
		final double E1=121;
		final double E2=897;
		final double E3=194;
		final double A=2.4e-3;

		double xi, a, b;
		ei=ef;

		if(xx<0.3){
		    a=0.11-2.4*xx;
		    b=-0.22+0.69*xx;
		}
		else{
		    a=-0.46-0.54*xx;
		    b=-0.01+0.01*xx;
		}
		xi=a+b*Math.log(ei)/Log10;
		return A*fpar[M][0]*Math.pow(ei,-fpar[M][1])*(0.05/(1+1.5*ei*xs/E2)+0.185/(1+1.5*ei*xs/E3)+
							    11.4*Math.pow(ei, xi)/(1+1.21*ei*xs/E1));
	    }
	case 0:
	default:
	    {
		final double E1=115;
		final double E2=850;
		final double pK=0.054;
		final double a=0.262;
		final double b=0.350e-3;
		final double as=-0.487;
		final double bs=8.766e-3;
		final double A=0.14;

		double zfi, zff, p, aux, aux1, aux2, auf1, auf2, f1, f2, f3, f4, de;
		ei=((a+b*ef)*Math.exp(b*xo)-a)/b;

		aux1=(bs*a-as*b);
		aux1*=2*aux1;
		aux2=2*b*b*b;
		zfi=(bs*ei*b*(-2*bs*a+4*as*b+bs*ei*b)+aux1*Math.log(a+ei*b))/aux2;
		zff=(bs*ef*b*(-2*bs*a+4*as*b+bs*ef*b)+aux1*Math.log(a+ef*b))/aux2;

		aux=1.1*xs;
		auf1=aux/E1;
		auf2=aux/E2;
		aux1=1/(1+auf1*ei);
		aux2=1/(1+auf2*ei);
		aux=aux1+pK*aux2;

		f1=Math.pow(ei, -fpar[M][1]-1)*aux;
		f2=Math.pow(ei, -fpar[M][1]-1)*(-(fpar[M][1]+1)*aux/ei-auf1*aux1*aux1-pK*auf2*aux2*aux2);
		f3=(-fpar[M][1]-1)*f2/ei+Math.pow(ei, -fpar[M][1]-1)
		    *((fpar[M][1]+1)*aux/(ei*ei)
		      +((fpar[M][1]+1)/ei)*(auf1*aux1*aux1+pK*auf2*aux2*aux2)
		      +2*auf1*auf1*aux1*aux1*aux1+2*pK*auf2*auf2*aux2*aux2*aux2);
		aux=f2/f1;
		f4=f3/f1-aux*aux;

		aux1=as+bs*ef;
		aux2=as+bs*ei;
		de=Math.exp(b*xo)-(f2/(2*f1))*(aux1*aux1-aux2*aux2)/(a+b*ef)-(f4/2)*(zff-zfi);

		ei=ei-(f2/(2*f1))*(zff-zfi);

		p=Math.exp(-xl/ei);
		return A*fpar[M][0]*p*de*Math.pow(ei,-fpar[M][1])*(1/(1+1.1*ei*xs/E1)+pK/(1+1.1*ei*xs/E2));
	    }
	}
    }
}