package seed.minerva.polarimetry;

import oneLiners.OneLiners;
import oneLiners.AsciiMatrixFile;
import seed.minerva.ConnectionPoint;
import seed.minerva.MinervaSettings;
import seed.minerva.NodeImpl;
import seed.minerva.StateFullNodeImpl;
import seed.minerva.diagnostics.Diagnostic;
import seed.minerva.diagnostics.LineOfSightSystem;
import seed.minerva.diagnostics.ServiceManager;
import seed.minerva.nodetypes.*;
import seed.minerva.physics.ElectronDensity;
import seed.minerva.physics.MagneticField;
import sun.reflect.generics.reflectiveObjects.NotImplementedException;
import algorithmrepository.Algorithms;



/** <pre><!-- for HTML'ified javadoc -->
 *  Based on:
 * "Review of plasma polarimetry" S. E. Segre
 *  
 * Stokes representation:
 * s_0 = cos 2chi . cos 2psi 
 * s_1 = cos 2chi . sin 2psi
 * s_2 = sin 2chi
 * 
 * Final stokes from initial:
 * sF = M . sI   
 * 
 * Equations 6.15, 6.6, 6.5, 6.8:
 * M =  / 	0		-W3		W2	\
 * 		|	W3		0		-W1	|
 * 		\	-W2		W1		0	/
 * 
 * W1 = c1l3 int[ ne . (Bx^2 - By^2)  dz ]
 * W2 = c1l3 int[ ne . Bx . By  dz ]
 * W3 = c3l2 int[ ne . Bz  dz ]
 * 
 * The 'analytical' prediction is based on these:
 * 
 * sF_0 = - W3  sI_1 + W2 sI_2
 * sF_1 = W3 sI_0 - W1 sI_2
 * sF_2 = -W2 sI_0 + W1 sI_1
 * 
 * It assumes that chi0 = 0:
 * (1)	cos 2psiF	=  - W3 . sin 2psi0
 * (2)	sin 2psiF	= W3 . cos 2psi0
 * (3)	sin 2chiF	= -W2 . cos 2psi0 + W1 sin 2psi0
 * 
 * For psi0 ~ 0:
 * DPsi = 1/2 W3		= 1/2 c3l2 int[ ne . Bz  dz ]
 * 2 chiF = - 1/2 W2	= - 1/2 c1l3 int[ ne . Bx . By  dz ]
 *  
 * For psi0 ~ 45':
 * cos (90' + 2Dpsi) 	= -W3, so:
 * DPsi = - 1/2 W3		= - 1/2 c3l2 int[ ne . Bz  dz ]
 * 2 chiF = 1/2 W1		= 1/2 c1l3 int[ ne . (Bx^2 - By^2)  dz ]
 * (same as before, but -ve?)
 * 
 * The coordinates assumed here have z // with los, y perp to both los and toroidal dir and x ^ y = z.
 * So for any LOS in the poloidal plane, Bx = Bphi and By ~ 0, so chiF = 0 if psi0 = 0.
 * 
 * Final analytic is therefor:
 * DPsi =  1/2 c3l2 int[ ne . Bz  dz ]
 * chiF = 0 for psi0 ~ 0
 * chiF = 1/2 c1l3 int[ ne . Btor^2 dz ] for psi0 ~ 45'
 * 
 * Wasn't that fun.
 *  
 * </pre>
 */
public class PolarimeterPredictionNode extends StateFullNodeImpl {
	private static final double d2r = Math.PI / 180;
	private static final double r2d = 180 / Math.PI;
	private static final double c = 299792458;
	private static final double e = 1.60217646e-19;       //electron charge, not exponential!
	private static final double me = 9.10938188e-31;
	private static final double ep0 = 8.85418782e-12;
	private static final double c1 = e*e*e*e / (16 * Math.PI*Math.PI*Math.PI * ep0 * me*me*me * c*c*c*c);
	private static final double c3 = e*e*e / ( 4 * Math.PI*Math.PI * ep0 * me*me * c*c*c);
	
	public final static int PLASMA_MODEL_COLD = 1;
	public final static int PLASMA_MODEL_WARM_SEGRE = 2;
	public final static int PLASMA_MODEL_WARM_MIRNOV = 3;
	public final static int PLASMA_MODEL_WARM = PLASMA_MODEL_WARM_MIRNOV;

	
	// Parents
	ScalarND electronDensity, electronTemperature;
	LineOfSightSystem los;
	IntegerArray dPsiStatus, chiStatus;
	MagneticField magneticField;
	/** Angle between initial linear polarisation's E vector and the direction perp to LOS which is 
	 * most parallel to +ve toroidal dir */
	DoubleArray psi0, chi0;
	DoubleValue laserWavelength;    				//laser angular freq
	BooleanArray returnPass;
	
	// Properties 
	private int predEnableRequest[]; /* Diagnotic.ON, Diagnostic.OFF or Diagnostic.AUTO to take from channelOK */
	private int plasmaModel = PLASMA_MODEL_COLD;
	
	// State
	private int nChannels;
	private double dPsi[];
	private double chi[];
	private boolean doublePass[];
	//laser angular freq
	double w = Double.NaN;
	
	
	/** whether or not we should calc DPsi and chi - dependant on prediction enable and psiStatus/chiStatus */
	private boolean calcDPsi[], calcChi[];
	
	/** Perform the analytical integrals instead of the full evolution evaluation.
	 * The analytical version expects the initial polarisation to have chi0 ~ 0
	 * and psi0 ~ 0 or psi0 ~ 45 degrees.  
	 */
	private boolean approxAnalytic = false;	
	private boolean psiIs45Degrees[];
	
	/** points along LOS [channel][{x,y,z}][point] */
	private double losPoints[][][];	
	/** dists along LOS */
	private double[][] losDists;	
	/** unit vector along LOS [channel][{x,y,z}] */
	private double uPara[][];
	/** unit vector perp to both LOS and toroidal dir [channel][{x,y,z}] */
	private double uPerp[][];
	/** magnetic field on los in [channel][x,y,z][point] */
	private double[][][] Bxyz;
	/** magnetic field on los in [channel][mag,theta,alpha][point] */
	private double[][][] Bmta;
	/** electron density and temperature on los in [channel][point] */
	private double[][] ne, Te;
	/** cold Omega with ne = 1 [channel][point][para,perp,perpperp] */
	private double[][][] coldOmegaFact;
	/** addition to cold Omega to get warm omega, with ne=1 and tau=1 [channel][point][para,perp,perpperp] */
	private double[][][] warmOmegaFact;
	
	final double uPhi[] = new double[]{ 0, 1, 0 };
	
	private double neUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronDensity");
	private double TeUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronTemperature");
	
	
	public PolarimeterPredictionNode() {
		this("Polarimetry predictions");
	}
	
	public PolarimeterPredictionNode(String name) {
		super(name);
		
		addConnectionPoint(new ConnectionPoint("magneticField", MagneticField.class, false, getField("magneticField")));
		addConnectionPoint(new ConnectionPoint("electronDensity", ScalarND.class, false, getField("electronDensity")));
		addConnectionPoint(new ConnectionPoint("electronTemperature", ScalarND.class, false, getField("electronTemperature")));
		addConnectionPoint(new ConnectionPoint("los", LineOfSightSystem.class, false, getField("los")));
		addConnectionPoint(new ConnectionPoint("dPsiStatus", IntegerArray.class, false, getField("dPsiStatus")));
		addConnectionPoint(new ConnectionPoint("chiStatus", IntegerArray.class, false, getField("chiStatus")));
		addConnectionPoint(new ConnectionPoint("psi0", DoubleArray.class, false, getField("psi0")));
		addConnectionPoint(new ConnectionPoint("chi0", DoubleArray.class, false, getField("chi0")));
		addConnectionPoint(new ConnectionPoint("laserWavelength", DoubleValue.class, false, getField("laserWavelength")));
		addConnectionPoint(new ConnectionPoint("returnPass", BooleanArray.class, false, getField("returnPass")));
		
	}
	
	/** Calcs LOS vectors uPara and uPerp */
	private void updateLOS(){
		losPoints = los.getPoints();
		losDists = los.getDistances();
		nChannels = los.getNumLOS();
		doublePass = returnPass.getBooleanArray();
		
		uPara = new double[nChannels][];		
		uPerp = new double[nChannels][];
		
		for(int i=0; i < nChannels; i++){
			/* one vector parallel to los */
			uPara[i] = new double[] {	losPoints[i][0][1] - losPoints[i][0][0], 
										losPoints[i][1][1] - losPoints[i][1][0],
										losPoints[i][2][1] - losPoints[i][2][0]		};
										
			double dl = Math.sqrt((uPara[i][0] * uPara[i][0])
								 + (uPara[i][1] * uPara[i][1])
								 + (uPara[i][2] * uPara[i][2]));
			uPara[i][0] /= dl; 
			uPara[i][1] /= dl; 
			uPara[i][2] /= dl;
			
			/* and one perp to both los and toroidal dir: uPerp = (u x phi) / |u x phi| */
			uPerp[i] = new double[] { -uPara[i][2], 0, uPara[i][0] };
			dl = Math.sqrt((uPerp[i][0] * uPerp[i][0]) + (uPerp[i][2] * uPerp[i][2]));
			if(dl <= 0)
				throw new RuntimeException("LOS is parallel to toroidal dir (y) os polarisation is ill defined.");
				
			uPerp[i][0] /= dl; 
			uPerp[i][2] /= dl;
			
		}
	}
	
	/** Resets the channel prediction enable flags to all automatic */
	private void initChannelEnableRequests() {
		/* If it already exists and is the right size, leave as is */
		if(predEnableRequest != null && predEnableRequest.length == nChannels)
			return;
			
		predEnableRequest = new int[nChannels];		
		for(int i=0; i < nChannels; i++)
			predEnableRequest[i] = Diagnostic.AUTOMATIC;
	}
	

	/** Recalculates calcDPsi/calcChi - whether or not to predict a channel based on the current
	 * state of the predictionEnables and the data statuses dPsiStatus/chiStatus.
	 */
	private void recalcCalculationEnable() {
		int nChannels = los.getNumLOS();
		calcDPsi = new boolean[nChannels];		
		calcChi = new boolean[nChannels];
		int statusDPsi[] = dPsiStatus.getIntegerArray();
		int statusChi[] = chiStatus.getIntegerArray();

		for(int i=0; i < nChannels; i++){
		
			calcDPsi[i] = 
				((predEnableRequest[i] == Diagnostic.AUTOMATIC) && (statusDPsi[i] == Diagnostic.DATA_VALID ))
					|| ((predEnableRequest[i] == Diagnostic.ON));
		
			calcChi[i] = 
				((predEnableRequest[i] == Diagnostic.AUTOMATIC) && (statusChi[i] == Diagnostic.DATA_VALID ))
					|| ((predEnableRequest[i] == Diagnostic.ON));
		}
	
	}
	
	/** Checks that psi0 and chi0 are valid for the analytic approximated calcuation
	 * and sets psiIs45Degrees[] accordingly */
	private void checkInitialPolarisation(){
		double psi0[] = this.psi0.getDoubleArray();
		double chi0[] = this.chi0.getDoubleArray();
		psiIs45Degrees = new boolean[nChannels];
		
		for(int i=0; i < nChannels; i++){			
			if( (chi0[i]*chi0[i]) > (5*5*d2r*d2r)  )
				throw new IllegalArgumentException("Analytic approx not valid for chi0 != 0");
			
			/* Get psi0 in degrees in 0 < psi0 < 180 */
			double psiDegrees = psi0[i]*r2d;			
			while(psiDegrees < 0) psiDegrees += 180;
			while(psiDegrees > 180) psiDegrees -= 180;
									
			if( (psiDegrees > 40 && psiDegrees < 50) || (psiDegrees > 130 && psiDegrees < 140))
				psiIs45Degrees[i] = true;
			else if( psiDegrees < 5 || psiDegrees > 175 || (psiDegrees > 85 && psiDegrees < 95)) /* psi0 ~= 0 */
				psiIs45Degrees[i] = false;						
			else{
				System.err.println("WARNING: Analytic approx not valid for psi0 = "+psiDegrees+" != {0 or +/-45} on channel "+i+". Disabling this channel's prediction. ");
				calcDPsi[i] = false;
				calcChi[i] = false;
			}
		}
	}

	/** Predicts dPsi and chi from full integration. 
	 * ds.dz = W ^ s
	 * 
	 * W = 
	 * 
	 */  
	private void doFullPrediction() {
		double psi0[] = this.psi0.getDoubleArray();
		double chi0[] = this.chi0.getDoubleArray();
		
	
		dPsi = new double[nChannels];
		chi = new double[nChannels];
				
		for(int i=0;i<nChannels;++i) {
			if(!calcDPsi[i] && !calcChi[i]){
				dPsi[i] = Double.NaN;
				chi[i] = Double.NaN;
				continue;
			}
			
			double s[] = dPsiChiToStokes(psi0[i], chi0[i]);
			
			singlePass(i, s, 1);
			
			if(doublePass[i])
				singlePass(i, s, -1);
			
			double dPsiChi[] = stokesToDPsiChi(s);
			
			dPsi[i] = dPsiChi[0] - psi0[i];
			chi[i] = -dPsiChi[1];
			
		}
	}
	
	private void singlePass(int i, double s[], double dir){
		for(int jj=1; jj < losPoints[i][0].length; jj++){
			int j = (dir > 0) ? jj : (losPoints[i][0].length - jj);
			
			double dl = losDists[i][j] - losDists[i][j-1];
			
			double Omega[] = new double[]{
					ne[i][j] * neUnit * coldOmegaFact[i][j][0],
					ne[i][j] * neUnit * coldOmegaFact[i][j][1],
					ne[i][j] * neUnit * coldOmegaFact[i][j][2]	};
			
			if(plasmaModel == PLASMA_MODEL_WARM_SEGRE || plasmaModel == PLASMA_MODEL_WARM_MIRNOV){
				Omega[0] += ne[i][j] * neUnit * Te[i][j] * TeUnit * warmOmegaFact[i][j][0];
				Omega[1] += ne[i][j] * neUnit * Te[i][j] * TeUnit * warmOmegaFact[i][j][1];
				Omega[2] += ne[i][j] * neUnit * Te[i][j] * TeUnit * warmOmegaFact[i][j][2];
			}
			
			

			/* This next line fixes things, it is essentially a hack. */
			//Omega[2] = -Omega[2];   
			/* It relates in the end to the differences in definitions of rotations and ellipticities. Whether you call ellipticity
			 * clockwise or anticlockwise from view direction of travel or looking at it oncoming. Whether the rotations away from the 
			 * toroidal direction are in toward the core or outward. It is defined here in some way for the plasma as according to B[] 
			 * from the equilibrium provider and also in the calibration code ScrewLoose depending on how the half-wave plate position relates
			 * to the vessel.
			 * It's complicated and this makes it work. */
			
			doRotation(Omega, s, dl);
			
		}
	}
	
	/**
	 * s_0 = cos 2chi . cos 2psi 
	 * s_1 = cos 2chi . sin 2psi
	 * s_2 = sin 2chi
	 */ 
	private static final double[] dPsiChiToStokes(double psi, double chi) {
		return new double[]{
				Math.cos(2*chi) * Math.cos(2*psi),
				Math.cos(2*chi) * Math.sin(2*psi),
				Math.sin(2*chi)
			};
	}
	
	private static final double[] stokesToDPsiChi(double s[]) {
		//get the rotation and ellipticity as we leave the plasma
		double psi = Math.atan2( s[1], s[0]) / 2;  //rotation is relative to initial direction
		double chi = Math.asin( s[2] ) / 2; //ellipcticity is absolute
		return new double[]{ psi, chi };
	}


	private void prepMagField(){
		Bxyz = new double[nChannels][][];
		Bmta = new double[nChannels][][];
		
		for(int i=0; i < nChannels; i++){
			if(!calcDPsi[i] && !calcChi[i]){
				Bxyz[i] = null;
				Bmta[i] = null;
				continue;
			}
			
			Bxyz[i] = magneticField.magneticField(losPoints[i]);
			Bmta[i] = new double[3][losPoints[i][0].length];

			for(int j=1; j < losPoints[i][0].length; j++){				

				Bmta[i][0][j] = Math.sqrt( Bxyz[i][0][j]*Bxyz[i][0][j] +
						Bxyz[i][1][j]*Bxyz[i][1][j] +
						Bxyz[i][2][j]*Bxyz[i][2][j]);

				//		get magnitude of B and angles relative to LOS and toroidal dir
				//cos(theta) = u . B / |B|
				//tan(alpha) = Bperp / Bphi   (with Bperp = B . uPerp)

				double uBB = (uPara[i][0]*Bxyz[i][0][j] + uPara[i][1]*Bxyz[i][1][j] + uPara[i][2]*Bxyz[i][2][j]) / Bmta[i][0][j];
				if(uBB > 1)uBB = 1;
				else if(uBB < -1)uBB = -1;

				Bmta[i][1][j] = Math.acos( uBB);
				Bmta[i][2][j] = Math.atan2(Bxyz[i][0][j]*uPerp[i][0] + Bxyz[i][1][j]*uPerp[i][1] + Bxyz[i][2][j]*uPerp[i][2], Bxyz[i][1][j]);
			}
		}
	}


	private void prepColdOmegaFacts(){
		coldOmegaFact = new double[nChannels][][];
		
		for(int i=0;i<nChannels;++i) {
			if(!calcDPsi[i] && !calcChi[i])
				continue;
			
			coldOmegaFact[i] = new double[losPoints[i][0].length][3];

			for(int j=1; j < losPoints[i][0].length; j++)
				coldOmegaFact[i][j] = coldOmega(1.0, Bmta[i][0][j], Bmta[i][2][j], Bmta[i][1][j], w);
		}
	}

	private void prepWarmOmegaFacts(){
		warmOmegaFact = new double[nChannels][][];
		
		for(int i=0;i<nChannels;++i) {
			if(!calcDPsi[i] && !calcChi[i])
				continue;
			
			warmOmegaFact[i] = new double[losPoints[i][0].length][3];

			for(int j=1; j < losPoints[i][0].length; j++)
				warmOmegaFact[i][j] = warmOmegaTerm(1.0, 1.0, Bmta[i][0][j], Bmta[i][2][j], Bmta[i][1][j], w, plasmaModel);
			
		}
	}

	/** Rotate s around W by angle |W| dl */
	private static final void doRotation(double W[], double s[], double dl){
		double Wmag, Wmag2, a[],b[],c[],Ws;
		Wmag2 = W[0]*W[0] + W[1]*W[1] + W[2]*W[2];
		
		if(Wmag2 == 0)return; //no rotation
		
		Wmag = Math.sqrt(Wmag2);
		
		//a = (W ^ s) / Wmag
		a = new double[]{
				(   W[1]*s[2] - W[2]*s[1]) / Wmag,
			    ( - W[0]*s[2] + W[2]*s[0]) / Wmag,
				(   W[0]*s[1] - W[1]*s[0]) / Wmag
		};
		
		//b = s - (W . s)W / Wmag2
		Ws = W[0]*s[0] + W[1]*s[1] + W[2]*s[2]; //W . s
		b = new double[]{
				s[0] - (Ws*W[0]) / Wmag2,
				s[1] - (Ws*W[1]) / Wmag2,
				s[2] - (Ws*W[2]) / Wmag2
		};
		
		//c = (W . s)W / Wmag2
		c = new double[]{
				Ws * W[0] / Wmag2,
				Ws * W[1] / Wmag2,
				Ws * W[2] / Wmag2
		};
		
		int j;
		for(j=0;j<3;j++)
			s[j] = a[j]*Math.sin(Wmag*dl) + b[j]*Math.cos(Wmag*dl) + c[j];
				
	}
	
	/**
	 * Calculate Omega(less ne factor) for a wave travelling in the z dir through the field B
	 * @returns The vector Omega (the Stokes vector about which a polarisation will rotate) 
	 */ 
	/*private static final double[] coldOmega(double ne, double BMag, double alpha, double theta, double w){
		final double wpFact = 3182.60729652443; //factor for plasma freq = e^2 / (epsilon0 me)
		final double em = 175882017145.163; //factor for e- cyc freq = (e/m)
		final double c = 299792458; // speed of light
		
		double Bc[] = new double[]{ BMag*Math.cos(alpha)*Math.sin(theta),
									BMag*Math.sin(alpha)*Math.sin(theta),
									BMag*Math.cos(theta) };
		
		double wc2,wp2; //electron cyclotron and plasma frequencies	(squared)	
		//double BMag; // B magnitude
		double s,v[] = new double[3]; //scalar and vector parts of Omega respectivly
		
		wp2 = wpFact * ne; // plasma freq^2 
		wc2 = em * em * BMag * BMag ; //ion cyc freq^2
		
		//the scalar part of Omega
		s = wp2 / (2 * c * (w*w*w) * ( 1 - (wc2 / (w*w) )));
		
		//the vector part of Omega		
		v[0] = em * em * (Bc[0]*Bc[0] - Bc[1]*Bc[1]); // vx = (e/m)^2 * (Bx^2 - By^2)
		v[1] = em * em * 2 * Bc[0] * Bc[1]; //vy = (e/m)^2 * 2BxBy
		v[2] = 2 * w * em * Bc[2];
		
		return new double[]{ s*v[0], s*v[1], s*v[2] };		
	}*/
	
	private static final double[] coldOmega(double ne, double BMag, double alpha, double theta, double w){
		final double wpFact = 3182.60729652443; //factor for plasma freq = e^2 / (epsilon0 me)
		final double em = 175882017145.163; //factor for e- cyc freq = (e/m)
		final double c = 299792458; // speed of light
		
		double cosThe = Math.cos(theta);
		double sinThe = Math.sin(theta);
		double cos2Al = Math.cos(2 * alpha);
		double sin2Al = Math.sin(2 * alpha);
		
		
		double wc,wc2,wp2; //electron cyclotron and plasma frequencies	(squared)	
		//double BMag; // B magnitude
		double s, v[] = new double[3]; //scalar and vector parts of Omega respectivly
		
		wp2 = wpFact * ne; // plasma freq^2 
		wc = em * BMag ; //ion cyc freq
		wc2 = wc*wc ; //ion cyc freq^2
		
		double X = wp2 / (w*w);
		double Y = wc / w;
		double Y2 = Y*Y;
		
		//the scalar part of Omega
		s = X / (2 * c * w * ( 1 - Y2));
		
		//the vector part of Omega		
		v[0] = w * X / (2 * c * (1 - Y2)) * Y2 * sinThe * sinThe * cos2Al;
		v[1] = w * X / (2 * c * (1 - Y2)) * Y2 * sinThe * sinThe * sin2Al;
		v[2] = w * X / (2 * c * (1 - Y2)) * 2 * Y * cosThe;
		
		return v;		
	}
	
	
	
	/**
	 * Calculate Omega(less ne factor) for a wave travelling in the z dir through the field B
	 * @returns The vector Omega (the Stokes vector about which a polarisation will rotate) 
	 */ 
	private static final double[] warmOmegaTerm(double Ne, double Te, double B, double alpha, double theta, double w, int model){
		final double wpFact = 3182.60729652443; //factor for plasma freq = e^2 / (epsilon0 me)
		final double em = 175882017145.163; //factor for e- cyc freq = (e/m)
		final double c = 299792458; // speed of light
		final double tauFact = 1.95695136e-6; //to get from eV to tau ( 1eV / mc^2 )
		
		double wc2,wp2; //electron cyclotron and plasma frequencies	(squared)	
		double X,Y,Y2; //normalised
		double sinThe, cosThe, cos2Al, sin2Al;
		double tau;
		
		
		tau = tauFact * Te;
		wp2 = wpFact * Ne; // plasma freq^2 
		wc2 = em * em * B * B ; //ion cyc freq^2
		
		X = wp2 / (w*w);
		Y2 = wc2 / (w*w);
		Y = Math.sqrt(Y2);
		
		sinThe = Math.sin(theta);		
		cosThe = Math.cos(theta);
		sin2Al = Math.sin(2 * alpha);
		cos2Al = Math.cos(2 * alpha);
		
		
		double v[] = new double[3]; //Omega
		
		//Omega0
		v[0] = 0; //(w * X / (2 * c))  * ( Y2 * sinThe * sinThe * cos2Al );
		v[1] = 0; //(w * X / (2 * c))  * ( Y2 * sinThe * sinThe * sin2Al );
		v[2] = 0; // (w * X / (2 * c))  * ( 2 * Y * cosThe );
		
		//+ deltaOmega
		if(model == PLASMA_MODEL_WARM_SEGRE){
		//Segre's correction (wrong)
			v[0] += (3 * tau * w * X / c)  * ( 2 * Y2 * sinThe * sinThe * cos2Al );
			v[1] += (3 * tau * w * X / c)  * ( 2 * Y2 * sinThe * sinThe * sin2Al );
			v[2] += (3 * tau * w * X / c)  * ( Y * cosThe );
			
		}else if(model == PLASMA_MODEL_WARM_MIRNOV){
			//Mirnov's correction (more correct)		

			double N = 1; //TODO: Do something here or justify it				
			v[0] += (3 * tau * w * X / c)  * ( (2*N*N - (5.0/4.0)) * Y2 * sinThe * sinThe * cos2Al );
			v[1] += (3 * tau * w * X / c)  * ( (2*N*N - (5.0/4.0)) * Y2 * sinThe * sinThe * sin2Al );
			v[2] += (3 * tau * w * X / c)  * ( (N*N - (5.0/3.0)) * Y * cosThe );
		}else
			throw new IllegalArgumentException("Unknown warm plasma model.");
		return v;
		
	}
	
	
	/** Predicts dPsi and chi from simple analytic 1D integrals
	 * <pre> 
	 * DPsi =  1/2 c3l2 int[ ne . Bz  dz ]
	 * chiF = 0 for psi0 ~ 0
	 * chiF = 1/2 c1l3 int[ ne . Btor^2 dz ] for psi0 ~ 45'
	 * </pre>
	 * */  
	private void doAnalyticPrediction() {

		double lambda = laserWavelength.getDouble();
		double c1l3 = c1 * lambda*lambda*lambda;
		double c3l2 = c3 * lambda*lambda;
		
		double neUnit = ServiceManager.getInstance().getUnitManager().getUnit("ElectronDensity");
				
		dPsi = new double[nChannels];
		chi = new double[nChannels];
							
		for(int i=0;i<nChannels;++i) {
			if(!calcDPsi[i] && !calcChi[i]){
				dPsi[i] = Double.NaN;
				chi[i] = Double.NaN;
				continue;
			}	
			
			double dx = losDists[i][1] - losDists[i][0];
						 
			double neBPara[] = new double[losPoints[i][0].length];			
			double neBTor2[] = new double[losPoints[i][0].length];
			
			// Calc BPara (no ne at this point) - it's needed for both DPsi and chi
			if(calcDPsi[i]){
				for(int j=0; j < losPoints[i][0].length; j++){				
					neBPara[j] = (	(Bxyz[i][0][j] * uPara[i][0]) +
									(Bxyz[i][1][j] * uPara[i][1]) +
									(Bxyz[i][2][j] * uPara[i][2])	);
				
					neBPara[j] *= ne[i][j] * neUnit;
				}
			
				
				dPsi[i] = 0.5 * c3l2 * Algorithms.trapz(neBPara, dx);
				//a return pass would pick up as much again
				//(since it rotates the opposite way but propagates in the opposite direction)
				if(doublePass[i])
					dPsi[i] *= 2;
			}
			
			                      
			if (calcChi[i]){ // neBTor2 is only needed for chi
				// if passing twice, the ellipticity should be entirely removed (TODO: Check that again)
				if(psiIs45Degrees[i] && !doublePass[i]) {
					for(int j=0; j < losPoints[i][0].length; j++)
						neBTor2[j] = (Bxyz[i][1][j]*Bxyz[i][1][j]) * ne[i][j] * neUnit;
						
					chi[i] = 0.5 * c1l3 * Algorithms.trapz(neBTor2, dx);
					
				} else {
					chi[i] = 0;
				}
			}
			
		}	
				
	}	

	@Override
	public void updateState() {
		boolean warm = !approxAnalytic && (plasmaModel == PLASMA_MODEL_WARM_MIRNOV || plasmaModel == PLASMA_MODEL_WARM_SEGRE);
		if(isPropertyChanged("plasmaModel")){
			warmOmegaFact = null;
			dPsi = null;
		}
		
		if(losPoints == null || isAncestorChanged("los") || isAncestorChanged("returnPass")){
			updateLOS();
			initChannelEnableRequests();
			calcDPsi = null;
		}
		
		if(calcDPsi == null || isAncestorChanged("dPsiStatus") || isAncestorChanged("chiStatus")){
			recalcCalculationEnable();
			Bxyz = null;
			ne = null;
			dPsi = null;
			psiIs45Degrees = null;
		}
		
		if( psiIs45Degrees == null || isAncestorChanged("psi0") || isAncestorChanged("chi0") || isPropertyChanged("approxAnalytic")){
			if(approxAnalytic)
				checkInitialPolarisation();
			dPsi = null;
		}
		
		if(Double.isNaN(w) || isAncestorChanged("laserWavelength")){
			w = 2 * Math.PI * c / laserWavelength.getDouble();
			coldOmegaFact = null;	
			warmOmegaFact = null;
		}
		
		if( Bxyz == null || isAncestorChanged("magneticField")){
			prepMagField();
			coldOmegaFact = null;	
			dPsi = null;
		}
				
		if(!approxAnalytic && coldOmegaFact == null){
			prepColdOmegaFacts();
			dPsi = null;
		}
				
		if(warm && warmOmegaFact == null){
			prepWarmOmegaFacts();
			dPsi = null;
		}
		
		if(ne == null || isAncestorChanged("electronDensity")){
			ne = new double[nChannels][];
			for(int i=0; i < nChannels; i++)
				ne[i] = (calcDPsi[i] || calcChi[i]) ? electronDensity.eval(losPoints[i]) : null;			
			dPsi = null;
		}
		
		if(warm && (Te == null || isAncestorChanged("electronTemperature"))){
			Te = new double[nChannels][];
			for(int i=0; i < nChannels; i++)
				Te[i] = (calcDPsi[i] || calcChi[i]) ? electronTemperature.eval(losPoints[i]) : null;				
			dPsi = null;
		}
		
		if(dPsi == null){
			if(approxAnalytic)
				doAnalyticPrediction(); //Just that for now
			else
				doFullPrediction();				
		}
		
		
	}
	
	/** Return an array of Psi for all channels, NaN where the channel prediction disabled */
	public double[] getDPsi() { update(); return dPsi; }	
	/** Return an array of Chi for all channels, NaN where the channel prediction disabled */
	public double[] getChi() { update(); return chi;	}
	
	public void setPredictionEnableRequest(int predEnableRequest){
		for(int i=0; i < nChannels; i++)
			this.predEnableRequest[i] = predEnableRequest;
		setChanged("predEnableRequest");
	}
	
	public void setPredictionEnableRequest(int[] predEnableRequest){
		update();
		
		if(predEnableRequest.length != nChannels)
			throw new IllegalArgumentException("Have "+nChannels+" channels but were passed "+predEnableRequest.length+" enables");
		this.predEnableRequest = predEnableRequest;
		setChanged("predEnableRequest");
	}
	
	/** Enable or disable analytic approximations (Faraday Rotation / Cotton-Mouton) */
	public void setAnalyticApprox(boolean approxAnalytic){ this.approxAnalytic = approxAnalytic; setChanged("approxAnalytic"); }
	
	/** Set model for treatment of plasma electron temperature effects, one of PLASMA_MODEL_* */ 
	public void setPlasmaModel(int plasmaModel){ this.plasmaModel = plasmaModel; setChanged("plasmaModel"); }
		
}