/*******************************************************************************
 * Copyright (c) 2010 Attila Kovacs <attila_kovacs@post.harvard.edu>.
 * All rights reserved. 
 * 
 * This file is part of crush.
 * 
 *     crush is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by
 *     the Free Software Foundation, either version 3 of the License, or
 *     (at your option) any later version.
 * 
 *     crush is distributed in the hope that it will be useful,
 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *     GNU General Public License for more details.
 * 
 *     You should have received a copy of the GNU General Public License
 *     along with crush.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Contributors:
 *     Attila Kovacs <attila_kovacs@post.harvard.edu> - initial API and implementation
 ******************************************************************************/
// Copyright (c) 2009 Attila Kovacs 

package crush.sharc2;

import crush.*;
import nom.tam.fits.*;

import java.io.*;
import java.net.*;
import java.util.Arrays;

import util.ArrayUtil;
import util.HorizontalCoordinates;
import util.Unit;
import util.Util;
import util.Vector2D;

// TODO Split nod-phases into integrations...
public class Sharc2Integration extends Integration<Sharc2, Sharc2Frame> implements GroundBased {
	/**
	 * 
	 */
	private static final long serialVersionUID = 2387643745464766162L;
	
	public boolean hasExtraTimingInfo = false;
	private boolean isDetectorStage = false;

	
	public Sharc2Integration(Sharc2Scan parent) {
		super(parent);
	}

	@Override
	public void validate() {	
		reindex();
		
		// Incorporate the relative instrument gain (under loading) in the scan gain...
		gain *= instrument.averagePixelGain;
		
		// Tau is set here...
		super.validate();	
		
		// Add the residual offsets to the DAC values...
		removeDrifts(instrument, size(), true, true);
		
		boolean directTau = false;
		if(hasOption("tau")) directTau = option("tau").equals("direct"); 
		
		if(!directTau) {
			double measuredLoad = instrument.getLoadTemperature(); 
			if(!hasOption("excessload")) instrument.excessLoad = measuredLoad - getSkyLoadTemperature();
			System.err.println("   Excess optical load on bolometers is " + Util.f1.format(instrument.excessLoad) + " K. (not used)");
			double tauLOS = -Math.log(1.0 - measuredLoad / ((Sharc2Scan) scan).ambientT);
			System.err.println("   Tau from bolometers (not used):");
			printEquivalentTaus(tauLOS * scan.horizontal.sinLat);
		}
		
	}
	
	@Override
	public void fillGaps() {
		// trim gaps...
		if(hasOption("nogaps")) {
			String argument = option("nogaps").getValue();
			if(argument.length() == 0) trimToGap();
			else if(argument.equalsIgnoreCase("JSharc")) if(scan.creator.equals("JSharc")) trimToGap();
		}
		else super.fillGaps();		
	}
	
	@Override
	public void writeProducts() {
		super.writeProducts();
		
		String scanID = getScanID();
		// Needs 'excessload'
		if(hasOption("response.calc")) {
			instrument.calcGainCoefficients(getSkyLoadTemperature());
			try { instrument.writeGainCoefficients(CRUSH.workPath + File.separator + "response-" + scanID + ".dat"); }
			catch(IOException e) { System.err.println("   WARNING! Error writing nonlinearity coefficients: " + e.getMessage()); }
		}
	}
	
	public void printEquivalentTaus(double value) {
		System.err.println("   --->"
				+ " tau(225GHz):" + Util.f3.format(getTau("225GHz", value))
				+ ", tau(350um):" + Util.f3.format(getTau("350um", value))
				+ ", tau(LOS):" + Util.f3.format(value / scan.horizontal.sinLat)
				+ ", PWV:" + Util.f2.format(getTau("pwv", value)) + "mm"
		);		
	}
	
	// TODO ???
	@Override
	public void removeDrifts(ChannelGroup<?> channels, int targetFrameResolution, boolean robust, boolean quick) {
		super.removeDrifts(channels, targetFrameResolution, robust, quick);
		
		// Recalculate the scan gains...
		/*
		if(hasOption("nonlinearity")) {
			comments += "G0";
			gain /= instrument.averagePixelGain;
			instrument.calcAveragePixelGain();
			gain *= instrument.averagePixelGain;
		}
		*/
		
	}
	
	public void filter13Hz(ChannelGroup<Sharc2Pixel> channels) {
		System.err.println("   Filtering 13.4 Hz resonances.");
		
		double f = 13.38 * Unit.Hz;
		double nyquist = 0.5 / instrument.samplingInterval;
		
		// Make it s.t. notch filter has ~0.3 Hz width...
		if(nyquist < f) throw new IllegalStateException("   WARNING! Data no longer contains unaliased 13.38 Hz signals.");
		
		int n = ArrayUtil.getTruncatedFFTSize(2 * (int)Math.ceil(nyquist / (0.3 * Unit.Hz)));
		
		double[] response = new double[n];
		Arrays.fill(response, 1.0);
		
		int F = (int)Math.floor(f / nyquist * response.length);
		response[F] = 0.0;
		response[F+1] = 0.0;
		
		filter(channels, response);
		
		comments += "n";
	}
	
	@Override
	public void setTau() throws Exception {
		String source = option("tau").getValue().toLowerCase();
		
		try { setTau(source); }
		catch(Exception e) {
			if(source.equals("maitau")) if(hasOption("maitau.fallback")) {
				System.err.print("   WARNING! MaiTau lookup failed. ");
				source = option("maitau.fallback").getValue().toLowerCase();
				if(source.equals("maitau")) System.err.println("Deadlocked fallback option.");
				System.err.println("Falling back to '" + source + "'.");
				option("tau").value = source;
				setTau(source);
			}
			else throw e;
		}
	}
	
	public void setTau(String source) throws Exception {
		
		if(source.equals("direct")) setZenithTau(getDirectTau());
		else if(source.equals("maitau")) {
			try { setTau("350um", getMaiTau("350um")); }
			catch(NumberFormatException no350) { setTau("225GHz", getMaiTau("225GHz")); }
		}
		else if(source.equalsIgnoreCase("225GHz")) {
			if(hasOption("tau.225GHz")) setTau("225GHz", option("tau.225GHz").getDouble());
			else throw new IllegalArgumentException("tau.225GHz is undefined.");
		}
		else if(source.equalsIgnoreCase("350um")) {
			if(hasOption("tau.350um")) setTau("350um", option("tau.350um").getDouble());
			else throw new IllegalArgumentException("tau.350um is undefined.");
		}
		else setZenithTau(Double.parseDouble(source));
			
		printEquivalentTaus(zenithTau);
		
		// TODO move to obslog...
		double tauLOS = zenithTau / scan.horizontal.sinLat;
		System.err.println("   Optical load is " + Util.f1.format(((Sharc2Scan) scan).ambientT * (1.0 - Math.exp(-tauLOS))) + " K.");
	
	}
	
	
	public void trimToGap() {
		Sharc2Frame first = getFirstFrame();
		Sharc2Frame last = getLastFrame();
		
		int t = first.index + 1;
		
		for(; t<=last.index; t++) {
			Sharc2Frame exposure = get(t);
			if(exposure != null) {
				int dN = (t - first.index);
				if(exposure.MJD > first.MJD + (dN + 1) * instrument.samplingInterval) break;
				if(hasExtraTimingInfo) {
					if(exposure.frameNumber > first.frameNumber + (dN + 1)) break;
					if(exposure.dspTime > first.dspTime + (dN + 1) * instrument.samplingInterval) break;
				}
			}
		}
		
		if(t < last.index) {
			System.err.println("   WARNING! Gap detected. Discarding data after gap.");
			this.trim(0, t);
		}
	}
	
	
	@Override
	public Sharc2Frame getFrameInstance() {
		return new Sharc2Frame((Sharc2Scan) scan);
	}

	

	// TODO Move higher up in the hierarchy...
	/*
	public void detectChopper() {		
		if(isAutoChopRecognition()) {
			if(scan.serialNo < 8156) { scan.isChopped = false; return; }
			System.out.print("   Detecting Chopper. ");
		}
		else if(!scan.isChopped) return;
		else System.out.print("   Chopped Reduction. ");

		// First calculate the mean
		double mean = 0.0;
		for(Frame exposure : this) if(exposure != null) mean += exposure.chopperPosition.x;
		mean /= size();

		// Now calculate the mean hi and low windows
		double hi=0.0, lo=0.0;
		int nhi=0, nlo=0;

		for(Frame exposure : this) if(exposure != null) {
			double chopOffset = exposure.chopperPosition.x;
			if(chopOffset > mean) { hi += chopOffset; nhi++; }
			else if(chopOffset < mean) { lo += chopOffset; nlo++; }
		}

		// nhi and nlo will be 0 if chopOffset is identical
		// in that case, the hi and lo windows are simply the identical chopper values.
		if(nhi > 0) hi /= nhi; else hi=mean;
		if(nlo > 0) lo /= nlo; else lo=mean;	

		// Measure the chop throw
		scan.chopperThrow = hi - lo;

		if(isAutoChopRecognition()) {
			if(scan.chopperThrow > 5.0 * Unit.arcsec) {
				scan.isChopped = true;
				System.out.print("OK. ");
			}
			else {
				scan.isChopped = false;
				System.out.println("No Chopper.");
				scan.chopperThrow = 0.0;
				return;
			}
		}


		// Remove the center from the chopper offsets
		System.out.print("Centering. "); 
		final float center = (float)(hi + lo) / 2.0F;
		for(Frame exposure : this) if(exposure != null) exposure.chopperPosition.x -= center;

		hi -= center;
		lo -= center;

		double firstMJD = Double.NaN, lastMJD = Double.NaN;
		int transitions = 0;
		boolean isHi = get(0).chopperPosition.x > center;

		for(int t=1; t < size(); t++) {
			Frame exposure = get(t);

			if(exposure != null) {
				double chopOffset = exposure.chopperPosition.x;

				// Check for transition
				if((isHi && chopOffset < lo / 2.0) || (!isHi && chopOffset > hi / 2.0)) {
					if(transitions == 0) firstMJD = exposure.MJD;
					else lastMJD = exposure.MJD;

					transitions++;
					isHi = !isHi;
				}
			}
		}

		scan.chopperFrequency =  (transitions-1) / (2 * (lastMJD - firstMJD) * Unit.day);

		System.out.print("Throw " + Util.f1.format(scan.chopperThrow / Unit.arcsec) + " \". ");
		System.out.print("Frequency " + Util.f1.format(scan.chopperFrequency * Unit.sec) + " Hz. "); 

		System.out.println();
	}
	 */
	
	@Override
	public int getNodPhase() {
		// TODO Auto-generated method stub
		return 0;
	}
	
	
	protected void read(BasicHDU[] HDU, int firstDataHDU) throws IllegalStateException, HeaderCardException, FitsException {
		Sharc2Scan sharcscan = (Sharc2Scan) scan;
		
		int nDataHDUs = HDU.length - firstDataHDU, records = 0;
		for(int datahdu=0; datahdu<nDataHDUs; datahdu++) records += HDU[firstDataHDU + datahdu].getAxes()[0];

		System.err.println(" Processing scan data:");
		
		System.err.println("   " + nDataHDUs + " HDUs,  " + records + " x " +
				(int)(instrument.integrationTime/Unit.ms) + "ms frames" + " -> " + 
				Util.f1.format(records*instrument.integrationTime/Unit.min) + " minutes total."); 
		
		clear();
		ensureCapacity(records);
		
		final Vector2D equatorialOffset = new Vector2D();
		double chopCenter = 0.0;
			
		for(int n=0; n<nDataHDUs; n++) {
			BinaryTableHDU hdu = (BinaryTableHDU) HDU[firstDataHDU+n]; 
			records = hdu.getAxes()[0];
			
			//final int iRA = hdu.findColumn("RA");
			//final int iDEC = hdu.findColumn("DEC");
			final int iAZ = hdu.findColumn("AZ");
			final int iEL = hdu.findColumn("EL");
			final int iAZE = hdu.findColumn("AZ_ERROR");
			final int iELE = hdu.findColumn("EL_ERROR");
			final int iPA = hdu.findColumn("Parallactic Angle");
			final int iUT = hdu.findColumn("UT");
			final int iLST = hdu.findColumn("LST");
			final int iDT = hdu.findColumn("Detector Time");
			final int iSN = hdu.findColumn("Sequence Number");
			final int iRAO = hdu.findColumn("RAO");
			final int iDECO = hdu.findColumn("DECO");
			final int iAZO = hdu.findColumn("AZO");
			final int iELO = hdu.findColumn("ELO");
			//final int iFlag = hdu.findColumn("Celestial");
			final int iCHO = hdu.findColumn("CHOP_OFFSET");
			//final int iTRK = hdu.findColumn("Tracking");
			
			hasExtraTimingInfo = iDT >= 0 && iSN >= 0; 

			final boolean isDoubleUT = hdu.getHeader().getStringValue("TFORM" + (iUT+1)).equalsIgnoreCase("1D");
			
			for(int i=0; i<records; i++) {		
				Sharc2Frame frame = new Sharc2Frame(sharcscan);
				Object[] row = hdu.getRow(i);
				
				final double UT = isDoubleUT ? ((double[]) row[iUT])[0] * Unit.hour : ((float[])row[iUT])[0] * Unit.hour;
				frame.MJD = sharcscan.iMJD + UT / Unit.day;
				
				frame.parseData((float[][]) row[0]);
				
				// The Equatorial coordinates will be calculated at based on the scan's tracking center
				// and the horizontal offsets
				/*
				frame.equatorial = new EquatorialCoordinates(
						((float[])row[iRA])[0] * Unit.hourAngle, 
						((float[])row[iDEC])[0] * Unit.deg, scan.apparentEpoch);
				*/
				
				// Enforce the calculcation of the equatorial coordinates
				frame.equatorial = null;
				
				final double AZE = ((float[])row[iAZE])[0] * Unit.arcsec;
				final double ELE = ((float[])row[iELE])[0] * Unit.arcsec;
				
				frame.horizontal = new HorizontalCoordinates(
						((float[])row[iAZ])[0] * Unit.deg + AZE,
						((float[])row[iEL])[0] * Unit.deg + ELE);
				
				final double PA = ((float[])row[iPA])[0] * Unit.deg;
				frame.sinPA = Math.sin(PA);
				frame.cosPA = Math.cos(PA);
				
				frame.LST = ((float[])row[iLST])[0] * Unit.hour;

				if(hasExtraTimingInfo) {
					frame.dspTime = ((double[])row[iDT])[0] * Unit.sec;
					frame.frameNumber = ((int[])row[iSN])[0];
				}		
				
				frame.horizontalOffset = new Vector2D(
						((float[])row[iAZO])[0] * Unit.arcsec + AZE * frame.horizontal.cosLat,
						((float[])row[iELO])[0] * Unit.arcsec + ELE);
				
				frame.chopperPosition.x = ((float[])row[iCHO])[0] * Unit.arcsec;
				chopCenter += frame.chopperPosition.x;	
				
				// Add in the scanning offsets...
				frame.horizontalOffset.add(sharcscan.horizontalOffset);

				// Add in the equatorial sweeping offsets
				// Watch out for the sign of the RA offset, which is counter to the native coordinate direction
				equatorialOffset.setNative(((float[])row[iRAO])[0] * Unit.arcsec, ((float[])row[iDECO])[0] * Unit.arcsec);			
				frame.toHorizontal(equatorialOffset);
				frame.horizontalOffset.add(equatorialOffset);
					
				//Skip reading the relescope flags...
				//They aren't really used for anything...
				//for(int bit=0; bit<32; bit++) frame.telescopeFlags |= ((boolean[])row[iFlag+bit])[0] ? 1<<bit : 0;
				
				add(frame);
			}
		}
		
		if(size() > 0) chopCenter /= size();
			
		// Add to the offsets the centered chopper signal, then set the instrument rotation...
		for(Sharc2Frame frame : this) {
			frame.chopperPosition.x -= chopCenter;

			// The chopper position is uncentered so better leave it out unless needed...
			frame.horizontalOffset.x += frame.chopperPosition.x;
			
			// Add the chopper offet to the actual coordinates as well...
			frame.horizontal.x += frame.chopperPosition.x / frame.horizontal.cosLat;
		}

	}
	
	
	@Override
	public synchronized void detectorStage() {
		if(isDetectorStage) return;
		
		for(Frame exposure : this) if(exposure != null) for(Sharc2Pixel pixel : instrument) 
			exposure.data[pixel.index] /= instrument.rowGain[pixel.row];
		isDetectorStage = true;
	}

	@Override
	public synchronized void readoutStage() {
		if(!isDetectorStage) return;
		for(Frame exposure : this) if(exposure != null) for(Sharc2Pixel pixel : instrument) 
			exposure.data[pixel.index] *= instrument.rowGain[pixel.row];
		isDetectorStage = false;
	}
	

	public double getMaiTau(String id) throws IOException {
		// Return immediately if ID does not match 225GHz or 350um, which are the only values in
		// the Mai-Tau lookup at present
		if(!id.equalsIgnoreCase("225GHz") && !id.equalsIgnoreCase("350um")) 
			throw new IllegalArgumentException("No MaiTau lookup for '" + id + "'.");
		
		if(!hasOption("maitau.server")) 
			throw new IllegalArgumentException(" WARNING! MaiTau server not set. Use 'maitau.server' configuration key.");
		
		Socket tauServer = new Socket();
		tauServer.setSoTimeout(3000);
		tauServer.setTcpNoDelay(true); 
		tauServer.setPerformancePreferences(0, 1, 2); // connection time, latency, throughput
		tauServer.setTrafficClass(0x10); // low latency
		tauServer.connect(new InetSocketAddress(option("maitau.server").getValue(), 63225));
		
		PrintWriter out = new PrintWriter(tauServer.getOutputStream(), true);
		BufferedReader in = new BufferedReader(new InputStreamReader(tauServer.getInputStream()));

		while(in.read() != '>'); // Seek for prompt

		out.println("set noexit"); // Enter into interactive mode (do not disconnect after first command).

		while(in.read() != '>'); // Seek for prompt

		out.println("set " + id); // Select which tau value to query...

		while(in.read() != '>'); // Seek for prompt

		out.println("get tau " + scan.timeStamp); // Request tau for the specified date	
		
		double value = Double.NaN;
		
		try { 
			value = Double.parseDouble(in.readLine().trim());
			System.err.println("   Got MaiTau! tau(" + id + ") = " + Util.f3.format(value));
		}
		catch(NumberFormatException e) {}
		
		out.println("exit"); // Disconnect from Mai-Tau server 
		in.close();
		out.close();
		
		if(Double.isNaN(value)) throw new NumberFormatException("No " + id + " value for date in MaiTau database.");
		
		return value;
	}

	
	public void setTau(String id, double value) {
		Vector2D t = (Vector2D) getTauCoefficients(id);
		Vector2D sharct = (Vector2D) getTauCoefficients("sharc2");
		setZenithTau(sharct.x / t.x * (value - t.y) + sharct.y);
	}

	public double getTau(String id) {
		return getTau(id, zenithTau);
	}
	
	public double getTau(String id, double value) {
		Vector2D t = (Vector2D) getTauCoefficients(id);
		Vector2D sharct = (Vector2D) getTauCoefficients("sharc2");
		return t.x / sharct.x * (value - sharct.y) + t.y;
	}
	
	protected Vector2D getTauCoefficients(String id) {
		String key = "tau." + id.toLowerCase();
		
		if(!hasOption(key + ".a")) throw new IllegalStateException("   WARNING! " + key + " has undefined relationship.");
		
		Vector2D coeff = new Vector2D();
		coeff.x = option(key + ".a").getDouble();
		if(hasOption(key + ".b")) coeff.y = option(key + ".b").getDouble();
	
		return coeff;
	}
	
	// Reminder: eps_obs = eps * areaFactor
	public double getSkyLoadTemperature() {
		double transmission = 0.5 * (getFirstFrame().transmission + getLastFrame().transmission);
		return (1.0 - transmission) * ((Sharc2Scan) scan).ambientT;
	}


	public double getDirectTau() { 
		double eps = (instrument.getLoadTemperature() - instrument.excessLoad) / ((Sharc2Scan) scan).ambientT; 	
		return -Math.log(1.0-eps) * scan.horizontal.sinLat;
	}

	
}
