package seed.minerva.apps.imse;


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.LinkedList;
import java.util.List;
import java.util.Map.Entry;

import javax.management.RuntimeErrorException;

import otherSupport.ColorMaps;
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.AugMSESystem;
import seed.minerva.imse.IMSEOptics105_35;
import seed.minerva.imse.IMSEOptics135_50;
import seed.minerva.imse.IMSEOptics135_50_rescaled;
import seed.minerva.imse.IMSEOptics180_50;
import seed.minerva.imse.IMSEOptics75_25;
import seed.minerva.imse.IMSEOptics75_25_imgFielded;
import seed.minerva.imse.IMSEOptics75_25_imgFielded_fibreFielded;
import seed.minerva.imse.IMSEOptics75_25_rescaled;
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.collection.HitsCollector;
import seed.minerva.optics.collection.ImageCollector;
import seed.minerva.optics.collection.IntensityInfo;
import seed.minerva.optics.collection.IntersectionProcessor;
import seed.minerva.optics.collection.PlaneAngleInfo;
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 seed.minerva.optics.pointSpread.DualGaussianPSF;
import seed.minerva.optics.pointSpread.GaussianPSF;
import seed.minerva.optics.pointSpread.MiniImagePSF;
import seed.minerva.optics.pointSpread.PointSpreadBuilder;
import seed.minerva.optics.pointSpread.PointSpreadFunction;
import seed.minerva.optics.pointSpread.PointsPSF;
import svg.SVGSplitView3D;



/** Collect average and spread of ray directions through protection
 * glass, for faraday rotation calculation.
 * 
 * @author oliford
 */
public class CollectProtectionGlassAngles implements Runnable {
	final static double reHitTolerence = 1e-6;
	
	public static String outPath = MinervaOpticsSettings.getAppsOutputPath() + "/rayTracing/augImse/protectionGlass";
	
	/** n Rays per start point */ 
	public static int nRays = 50000;
	
	public static int nThreads = 4;
	
	/** Start points along beam length. roughly 0.25 to 0.80 is in image */
	/*public static double L0 =  0.200, L1 = 0.850, dL1 = 0.010; 
	public static int lGridStep = 5;
	public static double A0 = -0.300, A1 = 0.300, dA1 = 0.010, dA2 = 0.050; //*/	
	
	public static double L0 =  0.200, L1 = 0.850, dL1 = 0.045; 
	public static int lGridStep = 1;
	public static double U0 = -0.300, U1 = 0.300, dU1 = 0.050, dU2 = 0.050; //*/
		
	public static double L[] = OneLiners.linSpace(L0, L1, dL1);
	
	/** Start points in vertical-ish direction, roughly -0.3 to +0.27 is in image */
	public static double Uhr[] = OneLiners.linSpace(U0, U1, dU1); 
	public static double Ulr[] = OneLiners.linSpace(U0, U1, dU2);
	//public static double Ahr[] = new double[]{ 0 }, Alr[] = new double[]{ 0 };
	
	/** Start points in other beam perp direction */
	//public static double B[] = OneLiners.linSpace(-0.3, 0.3, 16);
	//public static double B[] = new double[]{ -0.1, 0.0, 0.1 };
	public static double B[] = new double[]{ -0.00 };
	
	public final static double wavelength = 653e-9;
	
	/** What rays to trace */
	public static double minIntensity = 0.01;
	public static boolean traceReflections = false;
	
	/** Optical system setup */
	public static AugMSESystem sys = new AugMSESystem(null); //new IMSEOptics75_25_imgFielded_fibreFielded());
	/**we require three optics;
	/* 1) The imaging plane, which will collect the light INTENSITY.
	/* 2) The polarisation sensitive plane, which collects the light polarisation.
	/* 3) The inital target to fire rays at.
	/* e.g. for the normal MSE system, 1 is thee fibre plane, 2 is the PEMs and 3 is the main mirror:
	*/
	//public static Plane imagePlane = sys.imseOptics.ccd;
	public static Plane imagePlane = sys.tubeOptics.fibreEnds;
	public static Plane polarisationPlane = sys.tubeOptics.PEMsFront;
	public static Element initRaysTarget = sys.mirrorBox.protectionCoverFront;
	
	/** Don't draw rays that don't hit these elements */
	public static Surface mustHitToDraw = sys.tubeOptics.fibreEnds;
	//public static Surface mustHitToDraw = imagePlane;
	
	/** Outputting, set null to ignore */
	public static VRMLDrawer vrmlOut; //set in main because outPath get changed
	public static SVGRayDrawing svg;
	private static int nDrawingSkipRays = 77;
	
	/** Colour table for rays */
	private static double col[][];
	
	/** Direct image builder */
	public static ImageCollector imgCollect = new ImageCollector( 
													-0.015, 0.015, 600, 
													-0.015, 0.015, 600);
	
	private static String vecInfoNames[] = new String[]{ 
		"meanVecX", "meanVecY", "meanVecZ", "meanAng", 
		"covarXX", "covarXY", "covarXZ", "covarYY", "covarYZ", "covarZZ",
		"meanGlassIncid"
	};
	
	//multithreading globals
	private static double avgGlassVectorInfo[][][];	
	private static double R[][], Z[][];
	private int threadID;
	
	public static void main(String[] args) { 
		
		outPath += "-" + (sys.imseOptics == null ? "null" : sys.imseOptics.getName());
		System.out.println("outPath = " + outPath);
		
		//vrmlOut = new VRMLDrawer(outPath + "/tube.vrml", 0.005);
		//public static VRMLDrawer vrmlOut = null;
		//svg = new SVGRayDrawing(outPath + "/tube-", new double[]{ -5, -5, -1, 5, 5, 1}, true);
		//public static SVGRayDrawing svg = null;
	
		System.out.println("Grid: L = " + L[0] + ":" + (L[0]-L[1]) + ":" + L[L.length-1] + " = " + L.length + " hr x" + lGridStep);
		System.out.println("Grid: hrU = " + Uhr[0] + ":" + (Uhr[0]-Uhr[1]) + ":" + Uhr[Uhr.length-1] + " = " + Uhr.length + "");
		System.out.println("Grid: lrU = " + Ulr[0] + ":" + (Ulr[0]-Ulr[1]) + ":" + Ulr[Ulr.length-1] + " = " + Ulr.length + "");
		if(B.length > 1)
			System.out.println("Grid: B = " + B[0] + ":" + (B[0]-B[1]) + ":" + B[B.length-1] + " = " + B.length + "");
		
		initDrawing();
				
		avgGlassVectorInfo = new double[vecInfoNames.length][Ulr.length][L.length];
		
		R = new double [Ulr.length][L.length];
		Z = new double [Ulr.length][L.length];
		
		runThreads();
		
		endDrawing();
		if(imgCollect != null)
			imgCollect.writeImage(outPath + "/image.bin");
		
		BinaryMatrixFile.mustWrite(outPath + "/R.bin", R, false);
		BinaryMatrixFile.mustWrite(outPath + "/Z.bin", Z, false);
		
		for(int k=0; k < avgGlassVectorInfo.length; k++){
			BinaryMatrixFile.mustWrite(outPath + "/glassVecInfo-"+vecInfoNames[k]+".bin", avgGlassVectorInfo[k], false);			
		}		
	}
	
	public CollectProtectionGlassAngles(int threadID) {
		this.threadID = threadID;
	}
	
	//threaded main program
	public void run() {
		RandomManager randMan = new RandomManager();
		long sumDT = 0;
		int nDT = 0;
		
		double uB[] = Util.reNorm(Util.cross(AugMSESystem.nbiUnit, new double[]{0,0,1}));
		double uA[] = Util.cross(uB, AugMSESystem.nbiUnit);
		
		for(int iL=0; iL < L.length; iL++){ //for each start pos along beam
			boolean isHighResColumn = ((iL % lGridStep) == 0); 
			double U[] = isHighResColumn ? Uhr : Ulr;
			for(int iU=0; iU < U.length; iU++){ //for each vertical pos
				if((iU % nThreads) != threadID)
					continue;
				for(int iB=0; iB < B.length; iB++){ //for each horiz. beam perp pos
					
					if(svg!=null && threadID == 0) 
						svg.getSVG3D().startGroup("start "+iL+", "+iU+", "+iB);
						
				
					double start[] = new double[]{
							AugMSESystem.nbiStart[0] + L[iL]*AugMSESystem.nbiUnit[0] + U[iU]*uA[0] + B[iB]*uB[0],
							AugMSESystem.nbiStart[1] + L[iL]*AugMSESystem.nbiUnit[1] + U[iU]*uA[1] + B[iB]*uB[1],
							AugMSESystem.nbiStart[2] + L[iL]*AugMSESystem.nbiUnit[2] + U[iU]*uA[2] + B[iB]*uB[2],
						};
					
					synchronized (R) {
						R[iU][iL] = FastMath.sqrt(start[0]*start[0] + start[1]*start[1]);
						Z[iU][iL] = start[2];
					}
					
					BinaryMatrixWriter vecsOut = new BinaryMatrixWriter(outPath + "/vecs-iL_"+iL+"-iU_"+iU+".bin", 3);
					
					int nHits = 0;
					double vec[][] = new double[nRays][3];
					
					double avgIncidenceAng = 0;
					
					long t0 = System.currentTimeMillis();
					for(int j=0;j<nRays;j++){
						
						RaySegment ray = new RaySegment();
						ray.startPos = start.clone();
						ray.dir = Tracer.generateRandomRayTowardSurface(randMan, ray.startPos, initRaysTarget);
						//ray.dir = Util.reNorm(Util.minus(initRaysTarget.getBoundarySphereCentre(), ray.startPos));
	
						double right[] = Util.cross(ray.dir, new double[]{ 0, 0, 1 });
						ray.up = Util.reNorm(Util.cross(right, ray.dir));
						
						//For now, we're going to ignore the problem of what the MSE emission
						//makes as an initial polarisation and just do things in terms of 'up' and 'right'
						ray.E0 = PointSpreadFunction.getInputStatesForMuellerCalc();
						ray.length = Double.POSITIVE_INFINITY;
						ray.wavelength = wavelength;
						ray.startHit = null;
										
						Tracer.trace(sys, ray, 1000, minIntensity, traceReflections);
							
						if(threadID == 0)
							drawRay((iL*Ulr.length)+iU*Ulr.length/U.length, ray);
						
						List<Intersection> hits = ray.getIntersections(imagePlane);
						
						for(Intersection hit : hits){
							Intersection glassHit = hit.incidentRay.findFirstEarlierIntersection(sys.mirrorBox.protectionCoverFront);
							
							vec[nHits][0] = glassHit.incidentRay.dir[0];
							vec[nHits][1] = glassHit.incidentRay.dir[1];
							vec[nHits][2] = glassHit.incidentRay.dir[2];
							
							avgIncidenceAng += Util.dot(glassHit.incidentRay.dir, glassHit.normal);
							
							vecsOut.writeRow(vec[nHits]);
							
							nHits++;
						}
								
						Pol.recoverAll();
					}
					
					vecsOut.close();
					
					long dt = System.currentTimeMillis() - t0;
					sumDT += dt;
					nDT++;
					
					avgIncidenceAng = FastMath.acos(avgIncidenceAng / nHits);
					
					System.out.println("iL=" + iL+", iU="+iU+", iB="+iB+": R="+R[iU][iL]+", Z="+Z[iU][iL] +
							", hits="+ nHits+"/"+nRays + "\t" + nHits + "\t<θ> = " + avgIncidenceAng);
					//System.out.println("time: this=" + dt + "\tavg=" + (sumDT/nDT));
					
					if(svg!=null && threadID == 0) svg.getSVG3D().endGroup();
					
					double vecStats[] = Util.vecStats(vec, nHits);
					for(int k=0; k < vecStats.length; k++) {
						avgGlassVectorInfo[k][iU][iL] = vecStats[k];
					}
					avgGlassVectorInfo[10][iU][iL] = avgIncidenceAng;

					
				}
			}
			if(isHighResColumn && imgCollect != null)
				synchronized (imgCollect) { //only sync on write, .nextIntersection() might be called
					imgCollect.writeImage(outPath + "/image.bin");
				}
			
		}
		
	}
	
	/** Returns the ray path of the axial ray, used for sequential tracing */ 
	private static List<Surface> getCentralRayPath() {
		
		// First fire a ray backwards down the optic axis of the IMSE system
		RaySegment backRay = new RaySegment(
				imagePlane.getCentre(), 
				Util.minus(new double[]{0,0,0}, imagePlane.getNormal()),
				wavelength);
		
		Tracer.trace(sys, backRay, 500, minIntensity, traceReflections);
		
		//find where it hits the 'beam plane'
		Intersection hit = backRay.getIntersections(sys.beamPlane).get(0);
				
		//create a new ray in the opposite dir
		RaySegment fwdRay = new RaySegment(
				hit.pos, Util.minus(new double[]{0,0,0}, hit.incidentRay.dir), 
				wavelength); 
		
		Tracer.trace(sys, fwdRay, 500, minIntensity, traceReflections);
		
		System.out.println("Ray Path: ");
		fwdRay.dumpPath();
		
		VRMLDrawer.dumpRay(outPath + "/initRay.vrml", sys, fwdRay);
		
		List<Surface> primaryRayPath = fwdRay.getPrimarySurfaceSequence();
		primaryRayPath.remove(sys.beamPlane);//some hit this and some don't, so don't require it for the sequential trace
		
		return primaryRayPath;
	}

	/* **************************** Drawing Stuff **************************** */

	public static void initDrawing(){
		if(vrmlOut != null){
			vrmlOut.setDrawPolarisationFrames(false);
			vrmlOut.setDrawOnlyStrongest(false);
			vrmlOut.setSkipRays(nDrawingSkipRays);
		}
		
		col = ColorMaps.alternating2D3x3(L.length, Ulr.length);
		
		if(svg != null){
			svg.generateLineStyles(col, 0.0005);
			svg.setSkipRays(nDrawingSkipRays);
		}
		
	}
	
	private static void drawRay(int startIdx, RaySegment ray) {
		if(svg == null && vrmlOut == null)
			return;
				
		if(ray.getIntersections(mustHitToDraw).size() <= 0)
			return;
		
		if(svg!=null) svg.drawRay(ray, startIdx);
		if(vrmlOut != null) vrmlOut.drawRay(ray, col[startIdx]);
	}
	
	private static void endDrawing(){
		
		if(svg != null){
			svg.drawElement(sys);
			svg.destroy();
		}
		
		if(vrmlOut != null){
			vrmlOut.drawOptic(sys);
			vrmlOut.destroy();
		}
	}
	
	private static void runThreads() {
		CollectProtectionGlassAngles rnr[] = new CollectProtectionGlassAngles[nThreads];
		Thread thrd[] = new Thread[nThreads];
		for(int iT=0; iT < nThreads; iT++){
			rnr[iT] = new CollectProtectionGlassAngles(iT);
			thrd[iT] = new Thread(rnr[iT]);
			thrd[iT].start();
		}
		
		int nAlive;
		do{
			try {
				Thread.sleep(100);
			} catch (InterruptedException e) { }
		
			nAlive = 0;
			for(int iT=0; iT < nThreads; iT++){
				if(thrd[iT].isAlive())
					nAlive++;
			}
			
		}while(nAlive > 0);
		
		System.out.println("All threads done");
	}

	
}