package seed.minerva.imse.old;


import jafama.FastMath;


import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.PrintWriter;
import java.net.InetAddress;
import java.net.UnknownHostException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map.Entry;

import javax.management.RuntimeErrorException;

import otherSupport.RandomManager;

import binaryMatrixFile.BinaryMatrixFile;
import binaryMatrixFile.BinaryMatrixWriter;

import algorithmrepository.Algorithms;
import algorithmrepository.LinearInterpolation1D;
import oneLiners.OneLiners;
import seed.minerva.MinervaOpticsSettings;
import seed.minerva.aug.mse.TubeOptics1;
import seed.minerva.optics.Util;
import seed.minerva.optics.tracer.Tracer;
import seed.minerva.optics.types.*;
import seed.minerva.optics.surfaces.*;
import seed.minerva.optics.drawing.AsciiOutForWendel;
import seed.minerva.optics.drawing.SVGRayDrawing;
import seed.minerva.optics.drawing.VRMLDrawer;
import seed.minerva.optics.interfaces.Absorber;
import seed.minerva.optics.interfaces.IsoIsoInterface;
import seed.minerva.optics.interfaces.IsoIsoStdFresnel;
import seed.minerva.optics.interfaces.Reflector;
import seed.minerva.optics.materials.IsotropicFixedIndexGlass;
import seed.minerva.optics.materials.SchottSFL6;
import seed.minerva.optics.materials.Vacuum;
import seed.minerva.optics.optics.*;
import svg.SVGSplitView3D;

public class AugImseOld {
	final static double reHitTolerence = 1e-6;
	
	String outPath = MinervaOpticsSettings.getAppsOutputPath() + "/rayTracing/augImse";
	
	double nbiStart[] = {  -1.8792922,  -1.4524737,    0  }; //Tangential
	double nbiUnit[] = {  0.9639529,  0.2519894, 0.0854169  }; //Tangential
	//double nbiStart[] = { -1.9022743 , -1.3847706  , 0 }; //Normal
	//double nbiUnit[] = {  0.9181567, 0.3869007, 0.0854169  }; //Normal
	
	double mirrorPos[] = new double[]{ -0.302726,  -2.22957,  -0.4 };
	double mirrorNorm[] = new double[]{ -0.876326,  -0.286449, 0.43};
	
	TubeOptics1 tubeOptics;
	Iris mirrorBoxIris;
	Disc mirrorBoxHole;
	Square mainMirror;	
	Box fusionInABox;	
	Optic all;
	
	FileOutputStream asciiOut;
	PrintWriter asciiWriter;
	int nAsciiLines;
	
	public AugImseOld() {
		
		//as built, the tube runs along the x axis from [0,0,0] at the centre of the main viewing mirror
		tubeOptics = new TubeOptics1();
			
		//mirror box, also with X along east central
		mirrorBoxIris = new Iris("mirrorBoxIris", new double[]{ 2.175, 0.058, -0.400 }, new double[]{ 0, -1, 0 }, 0.30, 0.065, Absorber.ideal());
		mirrorBoxHole = new Disc("mirrorBoxHole", new double[]{ 2.175, 0.060, -0.400 }, new double[]{ 0, -1, 0 }, 0.065, null, null, IsoIsoInterface.ideal());
		Util.rotateOnX(mirrorBoxIris, mirrorBoxHole.getCentre(), -23*Math.PI/180);
		Util.rotateOnX(mirrorBoxHole, mirrorBoxHole.getCentre(), -23*Math.PI/180);
		Util.rotateOnZ(mirrorBoxIris, new double[]{ 0, 0, 0 }, -9*2*Math.PI/32); //spin to East and 1/32th past (so that X is now pu the centre of the MSE sector)
		Util.rotateOnZ(mirrorBoxHole, new double[]{ 0, 0, 0 }, -9*2*Math.PI/32); //spin to East and 1/32th past (so that X is now pu the centre of the MSE sector)
		//mirrorBoxHoleG.rotateOnZ(new double[]{ 0, 0, 0 }, -9*2*Math.PI/32); //spin to East and 1/32th past (so that X is now pu the centre of the MSE sector)
		//rotate the tube up to the correct elevation angle
		Util.rotateOnY(tubeOptics, new double[]{ 0, 0, 0}, 9.8 * Math.PI / 180);
		
		//no wrotate the tube in the horizontal plane so that it's at the correct angle
		//our coordinate system always has x as North
		double ang = 
			-Math.PI/2 + //from North, come back around so that tube is out to the east
			-2*Math.PI/32 + //East is on the MSE sector edge, so the center of the sector is 1/32th further around (16 sectors)
			-3.11 * Math.PI / 180; //The actual angle of the tube from east - Measured from sheet ü3 (Tokamak with crude pipe)
			//-3.8 * Math.PI / 180; //The actual angle of the tube from east - Measured from sheet ü2 (Detailed pipe, no Tokamak)
			
		Util.rotateOnZ(tubeOptics, new double[]{ 0, 0, 0 }, ang);
		
		//now move it's 0,0,0 to the real mirror position
		tubeOptics.shift(mirrorPos);
		
		double mirrorHoriz[] = Util.reNorm(Util.cross(mirrorNorm, new double[]{ 0, 0, 1 }));
		
		mainMirror = new Square("mainMirror", mirrorPos, mirrorNorm, mirrorHoriz, 0.200, 0.100, Reflector.ideal());
		
		fusionInABox = new Box("fusionInABox", new double[]{-5,-5,-1}, new double[]{5,5,1}, null, Absorber.ideal());
		
		all = new Optic("all", new Element[]{ mirrorBoxHole, mirrorBoxIris, mainMirror, tubeOptics, /* fusionInABox*/ });
		
		OneLiners.makePath(outPath + "/dummy");
		
	}
		
	public static void main(String[] args) { 
		
		(new AugImseOld()).gridRays();
		//(new AugImseOptics()).focalScan();
		//(new AugImseOptics()).generateRandomDistribs(args);
		//(new AugImseOptics()).sumDistribs();
	} 

	/** Calculates distributions on the image plane, from random points in the emmission volume (the beam).
	 * The Gaussian approximations to those distributions are used by the camera forward model to generate
	 * the proper image.
	 */
	/*public void generateRandomDistribs(String args[]){
		int maxRayAttempts = 20000; //rays per emission point
		int maxHits = 1000;
		int nEmissionPoints = 10000000; //number of random emission points
		double nbiD0 = 0.1; //start of emission volume along beam
		double nbiD1 = 1.1; //end of emission volume  along beam
		double nbiSigmaR = 0.3; // 1 sigma for radius of emission volume.
		
		int nFullOut = 20; //number of emission points from which to output the full hit points list

		//some specific stuff for running this on a cluster
		multiProcSetup(args);

		BinaryMatrixWriter gaussOut = new BinaryMatrixWriter(outPath + "/genGauss.bin", 31);
		BinaryMatrixWriter fullHits = new BinaryMatrixWriter(outPath + "/genGauss-fullHits.bin", 19);
		
		double beamUp[] = Util.reNorm(Util.cross(nbiUnit, new double[]{ 0, 1, 0 }));  
		double beamRight[] = Util.reNorm(Util.cross(nbiUnit, beamUp));  
		
		ArrayList<double[]> hitList = new ArrayList<double[]>();
		ArrayList<double[]> startVecs = new ArrayList<double[]>();
		
		Square focalPlate = (Square) Element.findElement(tubeOptics.elems, "fibreEnds");
		int lineNo = 0;

		int i=0;
		while(i < nEmissionPoints){
		
			//generate a random point in the emission volume
			double d = RandomManager.instance().nextUniform(nbiD0, nbiD1);
			//double r = FastMath.sqrt(RandomManager.instance().nextUniform(0, 1)) * nbiR;
			double r = Math.sqrt(Math.abs(RandomManager.instance().nextNormal(0, nbiSigmaR)));
			double phi = RandomManager.instance().nextUniform(0, 2*Math.PI);
			double cosPhi = FastMath.cos(phi), sinPhi = FastMath.sin(phi);				
			double start[] = new double[]{
					nbiStart[0] + d * nbiUnit[0] + r*beamUp[0] * cosPhi + r*beamRight[0] * sinPhi,
					nbiStart[1] + d * nbiUnit[1] + r*beamUp[1] * cosPhi + r*beamRight[1] * sinPhi,
					nbiStart[2] + d * nbiUnit[2] + r*beamUp[2] * cosPhi + r*beamRight[2] * sinPhi
				};
			
			double startRPZ[] = new double[]{
					FastMath.sqrt(start[0]*start[0] + start[1]*start[1]),
					FastMath.atan2(start[1],  start[0]),
					start[2],
			};
			
			//calculate distance from this point, to the centre of the mirror's bounding sphere that
			//we're going to fire rays at
			double avgFireDir[] = new double[]{
					start[0] - mainMirror.getCentre()[0],
					start[1] - mainMirror.getCentre()[1],
					start[2] - mainMirror.getCentre()[2]
				};
		
			//the intensity that we put this point in at, has to compensate for both the generation fraction
			//and the fact that we didn't draw evenly over the distribution volume
			
			double emissionDistSq = avgFireDir[0]*avgFireDir[0] + avgFireDir[1]*avgFireDir[1] + avgFireDir[2]*avgFireDir[2];
			double targetSphereRadius = mainMirror.boundingSphereRadius();
		
			//calculate the fraction of the total light emission from this point, that would have
			//actually gone in that direction
			double genFrac = (targetSphereRadius * targetSphereRadius) / (4 * emissionDistSq * emissionDistSq);
			
			double drawProb = FastMath.exp(-r*r / (2*nbiSigmaR*nbiSigmaR)) / (nbiSigmaR * FastMath.sqrt(2*Math.PI));
			
		
			hitList.clear();
			startVecs.clear();
			
			int nRaysAttempted = 0;
			
			double clusterWorseCosAng = 1;
			double clusterMeanOfRayWorstCosAng = 1;
			for(nRaysAttempted=0; nRaysAttempted < maxRayAttempts; nRaysAttempted++) {
				
				double s[] = start.clone();
				
				//generate random rays toward the main mirror's bounding sphere
				double unit[] = Tracer.generateRandomRayTowardSurface(s, mainMirror);
				
				//trace it through until it hits an absorber or after 30 hits 
				RaySegment ray;
				Tracer.trace(s, unit.clone(), all, 30, null);
				
				Element hitElems[] = ray.getElements();
				
				if(hitElems.length < 4)
					continue;
				
				int fpIdx;				
				for(fpIdx=0; fpIdx < hitElems.length; fpIdx++){
					if(hitElems[fpIdx] == focalPlate)
						break;
				}
				
				if(fpIdx >= hitElems.length)
					continue;
				
				double hitPos[] = ray.getPositonsStore()[fpIdx];
				double hitProj[] = focalPlate.planeProjection(hitPos);
				
				hitList.add(hitProj);
				startVecs.add(unit);
				
				//find largest intersection angle
				double rayWorstCosAng = 1;
				double norms[][] = ray.getNormalsStore();
				double pos[][] = ray.getPositonsStore();
				for(int j=1; j < ray.getNSegments(); j++){			
					if(hitElems[j].getMatType() != Element.MATTYPE_TRANSPARENT || hitElems[j].getRIndexRatio() == 1.0)
						continue; //we're only interested in refracting interfaces
					
					//not interested in lenses 3A/B
					if(hitElems[j].name.startsWith("lens3"))
						continue;
					
					double segUnit[] = Tracer.reNorm(new double[]{
							pos[j][0] - pos[j-1][0], 
							pos[j][1] - pos[j-1][1], 
							pos[j][2] - pos[j-1][2], 
					   });
					double cosAng = - (norms[j][0]*segUnit[0] + norms[j][1]*segUnit[1] + norms[j][2]*segUnit[2]); 
					System.out.println(hitElems[j].name + "\t" + FastMath.acos(cosAng)*180/Math.PI);
					if(cosAng < rayWorstCosAng)
						rayWorstCosAng = cosAng;
				}
				
				if(nFullOut > 0){
					fullHits.writeRow(d, r, phi, start, startRPZ, hitList.size(), nRaysAttempted, genFrac, drawProb, hitProj, hitPos, FastMath.acos(rayWorstCosAng));
				}
				
				if(rayWorstCosAng < clusterWorseCosAng)
					clusterWorseCosAng = rayWorstCosAng;
				
				clusterMeanOfRayWorstCosAng += rayWorstCosAng;
				
				//got enough hits already
				if(hitList.size() >= maxHits)
					break;
			}
			
			int nHits = hitList.size();
			System.out.println(i + "/" + nEmissionPoints + ": d="+d+", r=" + r + ", nHit="+ nHits +" / " + nRaysAttempted + "\tnFullOut = " + nFullOut + "\tclusterMeanOfRayWorstAng = " + FastMath.acos(clusterMeanOfRayWorstCosAng/nHits)*180/Math.PI);
			
			if(startVecs.size() != nHits)
				throw new RuntimeException("Sanity Failure: I can't count.");
			
			if(nHits <= 20)
				continue;
			
			double mc[] = new double[5];
			for(double p[] : hitList){
				mc[0] += p[0];
				mc[1] += p[1];
			}
			mc[0] /= nHits;
			mc[1] /= nHits;
			for(double p[] : hitList){
				mc[2] += (p[0] - mc[0])*(p[0] - mc[0]);
				mc[3] += (p[0] - mc[0])*(p[1] - mc[1]);
				mc[4] += (p[1] - mc[1])*(p[1] - mc[1]);
			}
			mc[2] /= nHits;
			mc[3] /= nHits;
			mc[4] /= nHits;
			
			//now collected LOS-like direction info

			double losVec[] = new double[3]; 
			for(double u[] : startVecs){
				losVec[0] += u[0];
				losVec[1] += u[1];
				losVec[2] += u[2];
			}
			losVec[0] /= nHits;
			losVec[1] /= nHits;
			losVec[2] /= nHits;
			
			double losRight[] = Tracer.reNorm(Tracer.cross(losVec, new double[]{ 0, 0, 1})); //LOS right
			if(FastMath.sqrt(losRight[0]*losRight[0] + losRight[1]*losRight[1] + losRight[2]*losRight[2]) == 0)
				throw new RuntimeException("Sanity: Not expecting a vertical average LOS");
			
			double losUp[] = Tracer.reNorm(Tracer.cross(losRight, losVec));
			double stdevRayDotR = 0;
			double stdevRayDotU = 0;
			for(double u[] : startVecs){
				double uDotR = u[0]*losRight[0] + u[1]*losRight[1] + u[2]*losRight[2];
				double uDotU = u[0]*losUp[0] + u[1]*losUp[1] + u[2]*losUp[2];
				stdevRayDotR += uDotR*uDotR;
				stdevRayDotU += uDotU*uDotU;
			}
			stdevRayDotR /= nHits;
			stdevRayDotU /= nHits;
			
			//now, for information, the stdev of the angles away from that
			
			gaussOut.writeRow(/* 1 / d, r, phi, start, startRPZ, 
							  /* 10 / nHits, nRaysAttempted, genFrac, drawProb, 
							  /* 14 / mc, 
							  /* 19 / losVec, losRight, losUp, 
							  /* 28 / stdevRayDotU, stdevRayDotR,
							  /* 30 / clusterWorseCosAng, clusterMeanOfRayWorstCosAng/nHits);
			
			if(nFullOut > 0)
				nFullOut--;
			
			if( ((i+1) % 20) == 0) 
				gaussOut.flush(); //gentle reminder to commit JavaOneLineUtils on NervousEnergy
			
			i++;
		}

		
	}
	*/
	
	private void multiProcSetup(String[] args) {
		int localID = -1;
		
		for(int i=0; i < args.length; i++){
			if (args[i].equalsIgnoreCase("--localid")) {
				localID = Integer.parseInt(args[i+1]);
			}
		}
		
		String hostName;
		
		System.out.print("Getting hostname... ");
		try{
			hostName = InetAddress.getLocalHost().getCanonicalHostName();
			System.out.println("OK: " + hostName);
		} catch (UnknownHostException e2) {
			hostName = "ERROR_GETTING_HOSTNAME";
			System.out.println("FAILED.");
		}
		
		if(localID == -1)
			localID = hostName.hashCode() + (int)RandomManager.instance().nextUniform(0, Integer.MAX_VALUE);
		
		//some vague attempt at a random seed based both on time and on our machine ID, to stop  multiple cluster slaves using the same seed
		//we can use the global one, as it's only global to this process/machine.
		RandomManager.instance().setSeed(hostName.hashCode() + localID + (int)System.currentTimeMillis());
		
		outPath = outPath + "/" + localID;
		System.out.println("outPath = " + outPath);
		
		
	}

	public void sumDistribs(){

		//double d[][] = BinaryMatrixFile.mustLoad(outPath + "/genGauss-nSrc_1e5-nRays_1e3-drawProbCorrected.bin", false);
		double d[][] = BinaryMatrixFile.mustLoad(outPath + "/genGauss.bin", false);
		//double d[][] = 	BinaryMatrixFile.mustLoad("/afs/ipp-garching.mpg.de/u/oliford/minerva/results/rayTracing/augImse/genGauss.bin", false);
		
		int nX = 320;
		int nY = 240;
		double x0 = -0.007, x1 = 0.007;
		double y0 = -0.007, y1 = 0.007;
		
		double x[] = OneLiners.linSpace(x0, x1, nX);
		double y[] = OneLiners.linSpace(y0, y1, nY);
		
		double I[][] = new double[nX+1][nY+1];
		System.arraycopy(y, 0, I[0], 1, nY);
		for(int ipx=0; ipx < nX; ipx++)
			I[ipx+1][0] = x[ipx];
		
		double sum = 0;
		
		for(int i=0; i < d.length; i++){
			double mX = d[i][13];
			double mY = d[i][14];
			double sXX = d[i][15];
			double sXY = d[i][16];
			double sYY = d[i][17];
			double det = sXX * sYY - sXY * sXY;
			double intens = d[i][11] * (d[i][9]/d[i][10]);// / d[i][12];
			
			double cXX = sYY / det;
			double cXY = -sXY / det;
			double cYY = sXX / det;

			
			if(det <= 0)
				continue;
			
			for(int ipy=0; ipy < nY; ipy++){
				double dmy = mY - y[ipy];
				if(dmy*dmy > 7*sXX)
					continue;
				
				for(int ipx=0; ipx < nX; ipx++){
					double dmx = mX - x[ipx];
					if(dmx*dmx > 7*sYY)
						continue;
					
					double v = intens * FastMath.exp(-0.5 *  (cXX*dmx*dmx + 2*cXY*dmx*dmy + cYY*dmy*dmy)) / (2*Math.PI*FastMath.sqrt(det));
					I[ipx+1][ipy+1] += v;
					sum += v;
				}	
			}
			
			System.out.print(".");
			
			if((i+1) % 100 == 0)
				System.out.println();
			 
			if((i+1) % 500 == 0)
				System.out.println(i + " / " + d.length + ", sum = " + sum);
		}
		
		BinaryMatrixFile.mustWrite(outPath + "/genGaussSum.bin", I, false);
		
	}
	
	
	/** Used to focus the image plane.
	 * Output the image plane points from a single source point, with the image plane at a series of points. */
/*	public void focalScan(){
		int nRays = 500000;
		int nShifts = 1;
		double shift0 = -0.003;
		double shift1 = +0.003;
		double nbiD = 0.65;
		
		BinaryMatrixWriter hitsOut = new BinaryMatrixWriter(outPath + "/focalScan-hits.bin", 8);

		double start[] = new double[]{
				nbiStart[0] + nbiD * nbiUnit[0],
				nbiStart[1] + nbiD * nbiUnit[1],
				nbiStart[2] + nbiD * nbiUnit[2]
			};
		
		
		start = new double[]{
				start[0] - 0.0 * (start[0] - mirrorPos[0]), 
				start[1] - 0.0 * (start[1] - mirrorPos[1]), 
				start[2] - 0.0 * (start[2] - mirrorPos[2]), 
		};

		Square focalPlate = (Square) Element.findElement(tubeOptics.elems, "fibreEnds");
		Disc lens3ABack = (Disc) Element.findElement(tubeOptics.elems, "lens3ABack");
		double centre0[] =  focalPlate.getCentre();
		double normal[] =  focalPlate.normal;
		
		//double col[][] = jetColorMap(nShifts);
		double col[][] = jetColorMap(300);

		//startAsciiOut(outPath + "/tube-asciiLines.txt");
		
		for(int iS=0; iS<nShifts; iS++){
			double shift = (nShifts <= 1) ? (shift0+shift1)/2 : (shift0 + (shift1 - shift0) * iS / (nShifts - 1));
			
			double hits[][] = new double[2][nRays];
		
			focalPlate.setCentre(new double[]{
						centre0[0] + shift * normal[0],
						centre0[1] + shift * normal[1],
						centre0[2] + shift * normal[2],
					});
			
			int nHitFibrePlane = 0;
			for(int j=0;j<nRays;j++){
				
				//start the ray off heading randomly toward the first lens primary surface
				double s[] = start.clone();
				double unit[] = Tracer.generateRandomRayTowardSurface(s, mainMirror);
				
				//trace it through until it hits an absorber or after 30 hits 
				RayPath ray = Tracer.traceRay(s, unit, all, 30, null);
				
				Element hitElems[] = ray.getElements();
				
				if(hitElems == null || hitElems.length < 4 || hitElems[3] != tubeOptics.elems[1])
					continue;
				
				int cIdx = 0;
				for(int k=0; k < hitElems.length; k++){
					if(hitElems[k] == lens3ABack){
						double hitPos[] = ray.getPositonsStore()[k];
						double hitProj[] = lens3ABack.planeProjection(hitPos);
						double dx = hitProj[0];
						double dy = hitProj[1];
						cIdx = (int)(col.length * Math.sqrt(dx*dx+dy*dy) / lens3ABack.getRadius());
						if (cIdx >= col.length) cIdx = col.length-1;						
					}
					if(hitElems[k] == focalPlate){
						double hitPos[] = ray.getPositonsStore()[k];
						double hitProj[] = focalPlate.planeProjection(hitPos);
						
						
						hitsOut.writeRow(hitPos, col[cIdx], hitProj[0]+iS*0.002, hitProj[1]);
						hits[0][nHitFibrePlane] = hitProj[0];
						hits[1][nHitFibrePlane] = hitProj[1];

						asciiLineOut(ray.toLineData(), col[cIdx]);
						
						nHitFibrePlane++;
						break;
					}
				}
				
			}
			
			double mX=0, mY=0;
			for(int i=0;i < nHitFibrePlane; i++){
				mX += hits[0][i];
				mY += hits[1][i];
			}
			mX /= nHitFibrePlane;
			mY /= nHitFibrePlane;
			
			double mDist = 0;
			for(int i=0;i < nHitFibrePlane; i++){
				double dX = hits[0][i] - mX;
				double dY = hits[1][i] - mY;
				mDist += dX*dX + dY*dY;
			}
			mDist = FastMath.sqrt(mDist / nHitFibrePlane);
			
			System.out.println(shift +"\t" + nHitFibrePlane + "\t" + mDist);
		}
		
		//and the txt output
		for(double lineData[][] : all.draw()) {
			asciiLineOut(lineData, new double[]{ 0, 0.5, 0 });
		}
		
		closeAsciiOut();
		
	}
*/	
	/** Trace and draw all rays from a rgid of starting positions along the beam and above and below it */
	public void gridRays(){
		int nRays = 500;
		int nStartsL = 9;//90;
		int nStartsZ = 6; //60;
	
		double nbiD0 = 0.3;
		double nbiD1 = 1.0;
		
		double Z0 = -0.3;
		double Z1 = 0.3;
		
		SVGRayDrawing svg = null;//new SVGRayDrawing(outPath + "/tube-", new double[]{ -5, -5, -1, 5, 5, 1}, true);
		VRMLDrawer vrmlOut = new VRMLDrawer(outPath + "/tube.vrml", 0.005);
		if(vrmlOut != null){
			vrmlOut.setDrawPolarisationFrames(true);
			vrmlOut.setDrawOnlyStrongest(true);
		}
		
		BinaryMatrixWriter hitsOut = new BinaryMatrixWriter(outPath + "/gridRays-hits.bin", 12);
		HashMap<Element, Integer> termCounts = new HashMap<Element, Integer>();
		
		double col[][] = alternatingColorMap(nStartsL, nStartsZ);
		if(svg != null) svg.generateLineStyles(col, 0.0005);
		
		OneLiners.dumpArray(mainMirror.getCentre());
		OneLiners.dumpArray(mainMirror.getNormal());
		
		//Square focalPlate = (Square) Util.findElement(tubeOptics, "fibreEnds");
		Disc focalPlate = (Disc) Util.findElement(tubeOptics, "lens2Front");
		for(int iL=0; iL<nStartsL; iL++){
			for(int iZ=0; iZ < nStartsZ; iZ++){
				int i = iL * nStartsZ + iZ;
				if(svg!=null) svg.getSVG3D().startGroup("start "+iL+", "+iZ);
			
				double dL = (nStartsL <= 1) ? (nbiD0+nbiD1)/2 : (nbiD0 + ((double)iL * (nbiD1-nbiD0))/(nStartsL-1));
				double dZ = (nStartsZ <= 1) ? (Z0+Z1)/2 : (Z0 + (Z1 - Z0) * iZ / (nStartsZ-1));
				double start[] = new double[]{
						nbiStart[0] + dL * nbiUnit[0],
						nbiStart[1] + dL * nbiUnit[1],
						nbiStart[2] + dL * nbiUnit[2] + dZ
					};
				
				double startRPZ[] = new double[]{
						FastMath.sqrt(start[0]*start[0] + start[1]*start[1]),
						FastMath.atan2(start[1],  start[0]),
						start[2],
				};
				
				//if(iL == (nStartsL-1)){
				//	start = new double[]{ -2.1, -0.1, 0.6 };
				//}
				
				//termCounts.clear();
				double polAng[][] = new double[nRays][];
				
				
				int nHitFibrePlane = 0;
				for(int j=0;j<nRays;j++){
					
					//start the ray off heading randomly toward the first lens primary surface
					//double start[] = new double[]{ -1.475 - 0.2 + ((double)i * 0.4)/(nStarts-1), 0, 0 };
					
					RaySegment ray = new RaySegment();
					ray.startPos = start.clone();
					ray.dir = Tracer.generateRandomRayTowardSurface(ray.startPos, mainMirror);
					/*ray.dir = Util.reNorm(new double[]{
							mainMirror.getCentre()[0] - ray.startPos[0], 
							mainMirror.getCentre()[1] - ray.startPos[1], 
							mainMirror.getCentre()[2] - ray.startPos[2], 
						});//*/ 
					double right[] = Util.cross(ray.dir, new double[]{ 0, 0, 1 });
					ray.up = Util.reNorm(Util.cross(right, ray.dir));
					
					//We're effectivly interested in the light that would come from dipole 
					// oscillations in each of the x,y and z plane
					
					ray.E0 = new double[][]{ 
							{ray.up[0], 0, right[0], 0 },
							{ray.up[2], 0, right[2], 0 },
							{ray.up[1], 0, right[1], 0 },
					};
					
					ray.length = Double.POSITIVE_INFINITY;
					ray.wavelength = 653e-9;
					ray.startHit = null;
					
					Tracer.trace(all, ray, 30, 0.01, false);
												
					List<Intersection> mainMirrorHits = ray.getIntersections(mainMirror);
					
					if(mainMirrorHits.size() <= 0)
						continue;
					
					List<Intersection> focalPlaneHits = ray.getIntersections(focalPlate);
				
					if(focalPlaneHits.size() <= 0)
						continue;
						
					for(Intersection hit : focalPlaneHits) {
						double hitProj[] = focalPlate.posXYZToPlaneUR(hit.pos);
						hitsOut.writeRow(hit.pos, col[i], hitProj, startRPZ, dZ);
						hit.incidentRay.rotatePolRefFrame(focalPlate.getUp());
						
						polAng[nHitFibrePlane] = new double[]{
								Pol.psi(hit.incidentRay.E1[0]),
								Pol.psi(hit.incidentRay.E1[1]),
								Pol.psi(hit.incidentRay.E1[2])
						};
						
						nHitFibrePlane++;
						
					}
					
					if(svg!=null) svg.drawRay(ray, i);
					if(vrmlOut != null) vrmlOut.drawRay(ray, col[i]);
				}
				
				if(svg!=null) svg.getSVG3D().endGroup();
				
				if(nHitFibrePlane > 0){
					double mAngX=0, mAngY=0, mAngZ=0;
					for(int j=0; j < nHitFibrePlane; j++){
						mAngX += polAng[j][0];
						mAngY += polAng[j][1];
						mAngZ += polAng[j][2];
					}
					mAngX /= nHitFibrePlane;
					mAngY /= nHitFibrePlane;
					mAngZ /= nHitFibrePlane;
					
					double sAngX=0, sAngY=0, sAngZ=0;
					for(int j=0; j < nHitFibrePlane; j++){
						sAngX += FastMath.pow2(polAng[j][0] - mAngX);
						sAngY += FastMath.pow2(polAng[j][1] - mAngY);
						sAngZ += FastMath.pow2(polAng[j][2] - mAngZ);
					}
					sAngX = FastMath.sqrt(sAngX / nHitFibrePlane);
					sAngY = FastMath.sqrt(sAngY / nHitFibrePlane);;
					sAngZ = FastMath.sqrt(sAngZ / nHitFibrePlane);;
				
					System.out.println("  "+mAngX*180/Math.PI + "\t" + sAngX*180/Math.PI);
					System.out.println("  "+mAngY*180/Math.PI + "\t" + sAngY*180/Math.PI);
					System.out.println("  "+mAngZ*180/Math.PI + "\t" + sAngZ*180/Math.PI);
					
				}
				
				
				System.out.println(Math.sqrt(start[0]*start[0]+start[1]*start[1]) +"\t"  +dZ +"\t"+ nHitFibrePlane);
				
				
			}	
			
		}
		
		if(svg != null){
			svg.drawElement(all);
			svg.destroy();
		}
		if(vrmlOut != null){
			vrmlOut.drawOptic(all);
			vrmlOut.destroy();
		}
			
		
	}


	private static double[][] alternatingColorMap(int nStartsL, int nStartsZ) {
		double quad[][][] = {
				{ { 1, 1, 0 },  { 1, 0, 0 } }, // yel, red
				{ { 0, 0, 1 },  { 0, 1, 0 } }  // blu,  grn
			};
			
		                
		double col[][] = new double[nStartsL*nStartsZ][];
		for(int iL=0; iL<nStartsL; iL++){
			for(int iZ=0; iZ < nStartsZ; iZ++){
				int i = iL * nStartsZ + iZ;

				col[i] = quad[iL % 2][iZ % 2];
			}
		}
		return col;
	}
		
}