/*******************************************************************************
 * 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) 2007 Attila Kovacs 

package crush.sourcemodel;

import nom.tam.fits.*;
import nom.tam.util.*;



import java.io.*;
import java.util.*;

import crush.CRUSH;
import crush.GenericInstrument;
import crush.Instrument;
import crush.Scan;

import util.*;


public class SkyMap extends Image {
	public double[][] weight, count;
	
	public Vector<Scan<?, ?>> scans = new Vector<Scan<?, ?>>();
	
	public String commandLine;
	
	public int generation = 0;
	public double integrationTime = 0.0;	
	public double weightFactor = 1.0;
 
	public double filterBlanking = Double.NaN;
	public double clippingS2N = Double.NaN;
	
	public Unit janskyPerBeam = new Unit("Jy/beam", Double.NaN);
	public Unit fitsUnit = dataUnit;
	
	private double jansky = Double.NaN;
	
	private Vector2D reuseOffset = new Vector2D();

	public SkyMap() { Locale.setDefault(Locale.US); }
	
	public SkyMap(String fileName, double janskyPerDataUnit) throws FitsException, HeaderCardException, IOException { 
		this();
		setJansky(janskyPerDataUnit);
		read(fileName);		
	}

	public SkyMap(int i, int j) { 
		this();
		setSize(i, j);
	}
	
	@Override
	public Object clone() {
		SkyMap clone = (SkyMap) super.clone();
		clone.reuseOffset = new Vector2D();
		return clone;
	}

	@Override
	public Image copy() {
		SkyMap copy = (SkyMap) super.copy();

		// Make a copy of the fundamental data
		copy.weight = (double[][]) copyOf(weight);
		copy.count = (double[][]) copyOf(count);
		
		// Make a copy of the beam dependent units as these get recalculated....
		copy.janskyPerBeam = (Unit) janskyPerBeam.clone();
		copy.fitsUnit = dataUnit;
		
		return copy;
	}
	
	
	
	public void setJansky(double value) {
		jansky = value;
		janskyPerBeam.value = jansky/dataUnit.value * dataUnitPerBeam.value;
	}
	
	public double getJansky() {
		return jansky;
	}
	
	
	// Todo: move a general version to Unit...
	// SkyMap should just re-define jansky and dataUnit values, and the beam...
	public Unit getUnit(String name) {
		if(name.equals(dataUnitPerBeam.name)) return dataUnitPerBeam;
		else if(name.equals("Jy/beam")) return janskyPerBeam;
		else if(name.equals("jansky/beam")) return janskyPerBeam;
		else if(name.equals("Jy/bm")) return janskyPerBeam;
		else if(name.equals("jansky/bm")) return janskyPerBeam;
		else if(name.contains("/")) {
			int index = name.lastIndexOf('/');
			String baseUnit = name.substring(0, index);
			String area = name.substring(index + 1);
			Unit unit = new Unit(name, Double.NaN);
			if(area.equals("arcsec**2") || area.equals("arcsec2") || area.equals("arcsec^2") || area.equals("sqarcsec"))
				unit.value = getUnit(baseUnit).value / Unit.arcsec2;
			else if(area.equals("arcmin**2") || area.equals("arcmin2") || area.equals("arcmin^2") || area.equals("sqarcmin"))
				unit.value = getUnit(baseUnit).value / Unit.arcmin2;
			else if(area.equals("deg**2") || area.equals("deg2") || area.equals("deg^2") || area.equals("sqdeg"))
				unit.value = getUnit(baseUnit).value / Unit.deg2;
			else if(area.equals("rad**2") || area.equals("rad2") || name.equals("rad^2") || area.equals("sr"))
				unit.value = getUnit(baseUnit).value / Unit.sr;
			else if(area.equals("mas") || area.equals("mas2") || area.equals("mas^2") || area.equals("sqmas"))
				unit.value = getUnit(baseUnit).value / (Unit.mas * Unit.mas);
			else if(name.equals("pixel") || name.equals("pix"))
				unit.value = getUnit(baseUnit).value / (delta * delta);
			return unit;
		}
		// Else assume there is a standard multiplier in front, such as k, M, m, u...
		else {
			Unit unit = new Unit(name, Double.NaN);
			unit.value = Unit.getMultiplier(name.charAt(0)) * getUnit(name.substring(1)).value;
			return unit;
		}
	}
	
	@Override
	protected void calcBeamUnits() {
		super.calcBeamUnits();
		janskyPerBeam.value = jansky /dataUnit.value * dataUnitPerBeam.value;		
	}
	

	public void shortInfo() {
		System.err.println("\n\n  [" + sourceName + "]\n" + super.toString());
	}

	public void setUnit(String name) {
		if(name.equalsIgnoreCase(dataUnitPerBeam.name)) fitsUnit = dataUnitPerBeam;
		else if(name.equalsIgnoreCase(dataUnit.name)) fitsUnit = dataUnit;
		else if(name.equalsIgnoreCase("Jy/beam")) fitsUnit = janskyPerBeam;
		else System.err.println("ERROR Can't understand unit " + name);
	}
	

	@Override
	public void setSize(int i, int j) {
		super.setSize(i, j);
		weight = new double[i][j];
		count = new double[i][j];
	}

	@Override
	protected void crop(int imin, int jmin, int imax, int jmax) {
		double[][] oldweight = weight;
		double[][] oldcount = count;
		
		super.crop(imin, jmin, imax, jmax);
			
		for(int i=imin, i1=0; i<=imax; i++, i1++) for(int j=jmin, j1=0; j<=jmax; j++, j1++) {
			weight[i1][j1] = oldweight[i][j];
			count[i1][j1] = oldcount[i][j];
		}
	}


	@Override
	public void clear() {
		super.clear();
		for(int i=0; i<sizeX(); i++) {
			Arrays.fill(weight[i], 0.0);
			Arrays.fill(count[i], 0.0);
		}
	}

	@Override
	public void reset() {
		super.reset();
		integrationTime = 0.0;
	}
	
	@Override
	public double getWeight(int i, int j) { return weight[i][j]; }
	
	@Override
	public void scale(int i, int j, double factor) {
		data[i][j] *= factor;
		weight[i][j] /= factor * factor;
	}

	public final void addPointAt(Vector2D mapOffset, final double value, final double g, final double w, final double time) {
		addPointAt(indexOfdX(mapOffset.x), indexOfdY(mapOffset.y), value, g, w, time);
	}
	
	public final void addPointAt(SphericalCoordinates coords, final double value, final double g, final double w, final double time) {
		projection.project(coords, reuseOffset);
		addPointAt(indexOfdX(reuseOffset.x), indexOfdY(reuseOffset.y), value, g, w, time);
	}

	public final void addPointAt(final int i, final int j, final double value, final double g, final double w, final double time) {
		final double wg = w * g;
		data[i][j] += wg * value;
		weight[i][j] += wg * g;
		count[i][j] += time;
	}
	
	public synchronized void addDirect(SkyMap map, double w) {
		for(int i=0; i<sizeX(); i++) for(int j=0;j<sizeY(); j++) if(map.flag[i][j] == 0) {
			data[i][j] += w * map.weight[i][j] * map.data[i][j];
			weight[i][j] += w * map.weight[i][j];
			count[i][j] += map.count[i][j];
			flag[i][j] = 0;
		}
		integrationTime += map.integrationTime;
	}

	// Normalize assuming weighting was sigma weighting
	public void normalize() {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
			if(weight[i][j] <= 0.0) {
				data[i][j] = 0.0;
				weight[i][j] = 0.0;
				flag[i][j] = 1;
			} 
			else {
				data[i][j] /= weight[i][j];	
				flag[i][j] = 0;	    
			}	    
		}
	}

	public void sanitize() {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(weight[i][j] <= 0.0 || flag[i][j] != 0) {
			weight[i][j] = 0.0;
			flag[i][j] |= 1;
			data[i][j] = 0.0;
		}

	}
	
	public void applyCorrection(double filtering, double[][] significance) {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(weight[i][j] > 0.0) if(significance[i][j] <= clippingS2N) {
			data[i][j] /= filtering;
			weight[i][j] *= filtering * filtering;
		}	
	}

	// It's important to completely reset clipped points, otherwise they can be used...
	// One possibility for future is to raise flag only, then call sanitize()
	public void rmsClip(double maxRelativeRMS) {
		double[][] rms = getRMSImage().data;
		double maxRMS = maxRelativeRMS * new Image(rms).getMin();
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			if(rms[i][j] > maxRMS) flag[i][j] = 1;
	}

	public void exposureClip(double minRelativeExposure) {
		double minIntTime = minRelativeExposure * new Image(count).getMax();
		directExposureClip(minIntTime);
	}
	
	public void directExposureClip(double minIntTime) {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			if(count[i][j] < minIntTime) flag[i][j] = 1;
	}

	public void s2nClip(double level) {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			if(data[i][j] / getRMS(i,j) < level) flag[i][j] = 1;
	}

	
	public boolean[][] getMask(double minS2N, int minNeighbours) {
		boolean[][] mask = new boolean[sizeX()][sizeY()];
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) if(data[i][j] / getRMS(i,j) > minS2N) mask[i][j] = true;
		
		if(minNeighbours > 0) {
			boolean[][] cleaned = new boolean[sizeX()][sizeY()];
			
			for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(mask[i][j]) {
				int neighbours = -1;
				final int fromi = Math.max(0, i-1);
				final int toi = Math.min(sizeX(), i+1);
				final int fromj = Math.max(0, j-1);
				final int toj = Math.min(sizeY(), j+1);
				for(int i1=fromi; i1<toi; i1++) for(int j1=fromj; j1<toj; j1++) if(mask[i1][j1]) neighbours++;
				if(neighbours >= minNeighbours) cleaned[i][j] = true;
			}
			return cleaned;
		}
		
		return mask;
	}

	@Override
	public void convolve(double[][] beam) {
		double[][] beamw = new double[sizeX()][sizeY()];
		data = getConvolved(data, beam, beamw);
		count = getConvolved(count, beam, null);
		weight = beamw;
		this.beam = beam;
	}
	
	@Override
	public void fastConvolve(double[][] beam, int step) {
		double[][] beamw = new double[sizeX()][sizeY()];
		data = getFastConvolved(data, beam, beamw, step);
		count = getFastConvolved(count, beam, null, step);
		weight = beamw;
		this.beam = beam;
	}
	
	
	public void addGaussianAt(double flux, double w, double intTime, double dx, double dy, double FWHM, double beamRadius, float[][] renorm) {		
		FWHM /= delta;
		final double sigma = FWHM / Util.sigmasInFWHM;
		final double A = -0.5 / (sigma*sigma);
		final double i0 = indexOfdX(dx);
		final double j0 = indexOfdY(dy);
		final int di = (int)Math.ceil(beamRadius * FWHM);
		
		final int imin = Math.max(0, (int)Math.floor(i0-di));
		final int imax = Math.min(sizeX()-1, (int)Math.ceil(i0+di));
		final int jmin = Math.max(0, (int)Math.floor(j0-di));
		final int jmax = Math.min(sizeY()-1, (int)Math.ceil(j0+di));
		
		for(int i=imin; i<=imax; i++) {
			double di2 = i - i0;
			di2 *= di2;
			
			for(int j=jmin; j<=jmax; j++) {
				double dj2 = j - j0;
				dj2 *= dj2;
				
				final double a = Math.exp(A * (di2+dj2));
				
				data[i][j] += a * w * flux;
				weight[i][j] += a * w;
				count[i][j] += a * intTime;
				renorm[i][j] += (float) a;
			}
		}
	}
	
	public void addGaussianAt(SkyMap from, int i, int j, double FWHM, double beamRadius, float[][] renorm) {
		addGaussianAt(from.data[i][j], from.weight[i][j], from.count[i][j], from.dXofIndex(i), from.dYofIndex(j), FWHM, 1.0, renorm);
	}
	
	
	public Image getS2NImage() {
		Image image = getRMSImage();
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) image.data[i][j] = data[i][j] / image.data[i][j]; 
		return image;
	}
	
	public Image getRMSImage() {
		double[][] rms = new double[sizeX()][sizeY()];
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) rms[i][j] = getRMS(i,j);
		return getImage(rms);
	}
	
	public Image getWeightImage() { 
		Image image = getImage(weight);
		image.unit = new Unit("weight", 1.0 / (unit.value * unit.value));
		return image;
	}
	
	public Image getTimeImage() { 
		Image image = getImage(count);
		image.unit = Unit.get("s");
		return image;
	}
	
	public Image getImage(double[][] data) {
		Image image = (Image) clone();
		image.data = data;
		return image;		
	}

	public final double getRMS(final MapIndex index) {
		return getRMS(index.i, index.j);
	}
	
	public final double getRMS(final int i, final int j) {
		return 1.0 / Math.sqrt(weight[i][j]);		
	}
	
	public final double getS2N(final MapIndex index) {
		return getS2N(index.i, index.j);
	}
	
	public final double getS2N(final int i, final int j) {
		return data[i][j] / getRMS(i,j);		
	}
	
	public int[][] getSkip(double blankingValue) {
		int[][] skip = (int[][]) copyOf(flag);
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(data[i][j] / getRMS(i,j) > blankingValue) skip[i][j] = 1;
		return skip;
	}

	@Override
	public void filterAbove(double extendedFWHM, int[][] skip) {
		super.filterAbove(extendedFWHM, skip);
		double r = getImageFWHM() / extendedFWHM;
		weightScale(1.0 / (1.0 - r*r));
	}
	
	public void filterAbove(double extendedFWHM, double blankingValue) {
		filterAbove(extendedFWHM, getSkip(blankingValue));
	}

	
	//	 8/20/07 Changed to use blanking 
	//         Using Gaussian taper.
	//         Robust re-levelling at the end.
	public void fftFilterAbove(double FWHM, double blankingValue) {
		// sigma_x sigma_w = 1
		// FWHM_x sigma_w = 2.35
		// FWHM_x * 2Pi sigma_f = 2.35
		// sigma_f = 2.35/2Pi * 1.0/FWHM_x
		// delta_f = 1.0/(Nx * delta_x);
		// sigma_nf = sigma_f / delta_x = 2.35 * Nx * delta_x / (2Pi * FWHM_x)
		
		// Try to get an honest estimate of the extended structures using FFT (while blaning bright sources).
		// Then remove it from the original image...
		
		final double[][] extended = new double[sizeX()][sizeY()];
		final int[][] skip = getSkip(blankingValue);
		
		double avew = 0.0;
		int n=0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			double rms = getRMS(i,j);
			
			if(skip[i][j] > 0) extended[i][j] = 0.0;
			else {
				double w = 1.0 / (rms * rms);
				extended[i][j] = w * data[i][j];
				avew += w;
				n++;
			}
		}
		if(n > 0) avew /= n;	
		
		Complex[][] cdata = ArrayUtil.loadFFT(extended);
		ArrayUtil.uncheckedFFT(cdata, true);

		final int nx = cdata.length;
		final int ny = cdata[0].length;

		double sigmax = Util.sigmasInFWHM *  nx * delta / (2.0*Math.PI * FWHM);
		double sigmay = Util.sigmasInFWHM *  ny * delta / (2.0*Math.PI * FWHM);
		
		final double ax = -0.5/(sigmax*sigmax);
		final double ay = -0.5/(sigmay*sigmay);

		for(int fx=1; fx<nx; fx++) {
			final double axfx2 = ax*fx*fx;
			final Complex[] r1 = cdata[fx];
			final Complex[] r2 = cdata[nx-fx];
			
			for(int fy=1; fy<ny; fy++) {
				final double A = Math.exp(axfx2 + ay*fy*fy);
				final int fy1 = ny - fy;
				
				r1[fy].scale(A);
				r2[fy].scale(A);
				r1[fy1].scale(A);
				r2[fy1].scale(A);
			}
		}
		cdata[0][0].zero();
		cdata[1][0].scale(Math.exp(ax));
		cdata[0][1].scale(Math.exp(ay));
		
		ArrayUtil.uncheckedFFT(cdata, false);
		ArrayUtil.unloadFFT(cdata, extended);
		
		// Scale weights for filtering in the white noise assumption/
		double r = getImageFWHM() / FWHM;
		final double weightScale = 1.0 / (1.0 - r*r);
		
		if(avew > 0.0) for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
			data[i][j] -= extended[i][j] / avew;
			weight[i][j] *= weightScale;
		}

		if(Double.isNaN(extFilterFWHM)) extFilterFWHM = FWHM;
		else extFilterFWHM = 1.0/Math.sqrt(1.0/(extFilterFWHM * extFilterFWHM) + 1.0/(FWHM*FWHM));
		
		filterCorrect(instrument.resolution, skip);
	}
	

	public double getMeanIntegrationTime() {
		double sum = 0.0, sumw = 0.0;
		for(int i=0;i<sizeX(); i++) for(int j=0;j<sizeY(); j++) if(flag[i][j] == 0) {
			sum += count[i][j] * weight[i][j];
			sumw += weight[i][j];
		}
		return sum > 0.0 ? sum/sumw : 0.0;
	}
	
	@Override
	public void scale(final double scalar) {
		if(scalar == 1.0) return;		
		final double scalar2 = scalar*scalar;
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
			data[i][j] *= scalar;
			weight[i][j] /= scalar2;
		}
	}  

	@Override
	public double mean() {
		double sum = 0.0, sumw = 0.0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			final double w = weight[i][j];
			sum += w * data[i][j];
			sumw += w;
		}
		return sum / sumw;
	}
	
	@Override
	public double median() {
		WeightedPoint[] point = new WeightedPoint[countPoints()];

		if(point.length == 0) return 0.0;

		int k=0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			point[k++] = new WeightedPoint(data[i][j], weight[i][j]);
		
		return ArrayUtil.median(point);
	}
	
	public void rmsScale(final double scalar) {
		final double scalar2 = scalar*scalar;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) weight[i][j] /= scalar2;
	}

	public void weightScale(double scalar) {
		if(scalar == 1.0) return;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) weight[i][j] *= scalar;
	}

	public void weight(boolean robust) {
		double weightCorrection = robust ? getRobustChi2() : getChi2();
		weightScale(weightCorrection);
		weightFactor *= weightCorrection;
	}

	// Return to calculated weights...
	public void dataWeight() {
		weightScale(1.0/weightFactor);
		weightFactor = 1.0;
	}

	public double getChi2(boolean robust) {
		return robust ? getRobustChi2() : getChi2();
	}
	
	protected double getRobustChi2() {
		double[] chi2 = new double[countPoints()];

		if(chi2.length == 0) return 0.0;

		for(int i=0, k=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			double s2n = data[i][j] / getRMS(i,j);
			chi2[k++] = s2n * s2n;
		}

		// median(x^2) = 0.454937 * sigma^2 
		return 0.454937 / ArrayUtil.median(chi2);	
	}

	protected double getChi2() {
		double chi2 = 0.0;
		int n=0;
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			double s2n = data[i][j] / getRMS(i,j);
			chi2 += s2n * s2n;
			n++;
		}

		return n / chi2;	
	}
		
	public boolean containsNaN() {
		boolean hasNaN = false;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
			if(Double.isNaN(data[i][j])) hasNaN = true;
			if(Double.isNaN(weight[i][j])) hasNaN = true;
		}
		return hasNaN;
	}


	// derive the MEM correction if MEM modeling
	
	public void MEM(double lambda, boolean forceNonNegative) {
		//double[][] smoothed = getConvolvedTo(data, 5.0*instrument.resolution);
		double[][] smoothed = new double[sizeX()][sizeY()];
		//if(forceNonNegative) for(int x=0; x<sizeX(); x++) for(int y=0; y<sizeY(); y++) if(smoothed[x][y] < 0.0) smoothed[x][y] = 0.0;
		MEM(smoothed, lambda);
	}
		
	public void MEM(double[][] model, double lambda) {
		for(int x=0; x<sizeX(); x++) for(int y=0; y<sizeY(); y++) if(flag[x][y] == 0) {
			double sigma = Math.sqrt(1.0/weight[x][y]);
			double memValue = Math.hypot(sigma, data[x][y]) / Math.hypot(sigma, model[x][y]) ;
			data[x][y] -= Math.signum(data[x][y]) * lambda * sigma * Math.log(memValue);
		}
	}
	
	
	public Fits getFits() throws HeaderCardException, FitsException, IOException {
		FitsFactory.setUseHierarch(true);
		Fits fits = new Fits();	

		calcBeamUnits();
		
		fits.addHDU(makeHDU(fitsUnit.value, fitsUnit.name, "Intensity"));
		fits.addHDU(getImage(count).makeHDU(Unit.sec, "s", "Integration Time"));
		fits.addHDU(getRMSImage().makeHDU(fitsUnit.value, fitsUnit.name, "RMS"));
		fits.addHDU(getS2NImage().makeHDU(1.0, "1", "Signal-to-Noise"));

		editHeader(fits);
		
		for(Scan<?,?> scan : scans) fits.addHDU(scan.getHDU(instrument.options));

		return fits;
	}

	public Fits getFits(String name) throws FitsException, IOException {
		FitsFactory.setUseHierarch(true);
		Fits fits;
	
		try { 
			fits = new Fits(new File(name)); 
			fileName = name;
		}
		catch(Exception e) { 
			fileName = CRUSH.workPath + name;
			fits = new Fits(new File(fileName)); 
		}
		
		return fits;
	}

	public void read(String name) throws HeaderCardException, FitsException, IOException {
		
		// Get the coordinate system information
		BasicHDU[] HDU = getFits(name).read();
		readHeader(HDU[0].getHeader());
		readCrushData(HDU);
	}

	private void readBasicHeaderData(Header header) throws HeaderCardException {
		getCoordinateInfo(header);
		
		if(header.containsKey("INSTRUMEN")) {
			instrument = Instrument.forName(header.getStringValue("INSTRUMEN"));
			if(instrument == null) instrument = new GenericInstrument(header.getStringValue("INSTRUMEN"));
		}
			
		if(header.containsKey("TELESCOP")) ((GenericInstrument) instrument).telescope = header.getStringValue("TELESCOP");
		
		creator = header.getStringValue("CREATOR");
		if(creator == null) creator = "unknown";
		creatorVersion = "unknown";
		sourceName = header.getStringValue("OBJECT");
		
		// get the beam and calculate derived quantities
		if(header.containsKey("BEAM")) 
			instrument.resolution = header.getDoubleValue("BEAM", instrument.resolution / Unit.arcsec) * Unit.arcsec;
		else if(header.containsKey("BMAJ"))
			instrument.resolution =  header.getDoubleValue("BMAJ", instrument.resolution / Unit.deg) * Unit.deg;
		else 
			instrument.resolution = 3.0 * delta;

		smoothFWHM = delta / fwhm2size;
	}


	public void readHeader(Header header) throws HeaderCardException {	
		this.header = header;
		
		weightFactor = 1.0;
		correctingFWHM = Double.NaN;
		filterBlanking = Double.NaN;
		clippingS2N = Double.NaN;
		integrationTime = Double.NaN;
		extFilterFWHM = Double.NaN;
		
		readBasicHeaderData(header);
		readCrushHeaderData(header);
		
		fitsUnit = janskyPerBeam;
		setUnit(header.getStringValue("BUNIT"));
	}

	private void readCrushHeaderData(Header header) throws HeaderCardException {	
		
		weightFactor =  header.getDoubleValue("XWEIGHT", 1.0);
		integrationTime = header.getDoubleValue("INTEGRTN", Double.NaN) * Unit.s;
		
		String correction = header.getStringValue("CORRECTN");
		if(correction == null) correction = "None";
		try { correctingFWHM = Double.parseDouble(correction) * Unit.arcsec; }
		catch(NumberFormatException e) { correctingFWHM = Double.NaN; }
	
		filterBlanking = header.getDoubleValue("FLTRBLNK", header.getDoubleValue("MAP_XBLK", Double.NaN));
		clippingS2N = header.getDoubleValue("CLIPS2N", header.getDoubleValue("MAP_CLIP", Double.NaN));
			
		creatorVersion = header.getStringValue("CRUSHVER");
	
		// get the map resolution
		smoothFWHM = header.getDoubleValue("SMOOTH") * Unit.arcsec;
		if(smoothFWHM < delta / fwhm2size) smoothFWHM = delta / fwhm2size;
		extFilterFWHM = header.getDoubleValue("EXTFLTR", Double.NaN) * Unit.arcsec;		
	}

	private void readCrushData(BasicHDU[] HDU) throws FitsException {
		unit = fitsUnit;
		
		read(HDU[0]);
		getTimeImage().read(HDU[1]);
		
		Image noise = getWeightImage();
		noise.unit = unit;
		noise.read(HDU[2]);
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) weight[i][j] = 1.0 / (weight[i][j] * weight[i][j]); 
	}


	public void write() throws HeaderCardException, FitsException, IOException { write(getFits()); }

	public void write(Fits fits) throws HeaderCardException, FitsException, IOException {
		if(fileName == null) fileName = CRUSH.workPath + File.separator + sourceName + ".fits";  
		BufferedDataOutputStream file = new BufferedDataOutputStream(new FileOutputStream(fileName));
		fits.write(file);
		System.out.println(" Written " + fileName);
	}
	
	@Override
	public void editHeader(Fits fits) throws HeaderCardException, FitsException, IOException {
		super.editHeader(fits);
		
		nom.tam.util.Cursor cursor = fits.getHDU(0).getHeader().iterator();

		// Go to the end of the header cards...
		while(cursor.hasNext()) cursor.next();
	
		// Add the command-line reduction options
		if(commandLine != null) {
			StringTokenizer args = new StringTokenizer(commandLine);
			cursor.add(new HeaderCard("ARGS", args.countTokens(), "The number of arguments passed from the command line."));
			int i=1;
			while(args.hasMoreTokens()) Util.addLongFitsKey(cursor, "ARG" + (i++), args.nextToken(), "Command-line argument.");
		}
		
		cursor.add(new HeaderCard("V2JY", jansky, "The default conversion from V to Jy."));	
		cursor.add(new HeaderCard("INTEGRTN", integrationTime / Unit.s, "The total integration time in seconds."));
		
		if(!Double.isNaN(filterBlanking))
			cursor.add(new HeaderCard("FLTRBLNK", filterBlanking, "The S/N blanking of LSS filter."));
		if(!Double.isNaN(clippingS2N))
			cursor.add(new HeaderCard("CLIPS2N", clippingS2N, "The S/N clipping level used in reduction."));

		cursor.add(new HeaderCard("CORRECTN", 
				(Double.isNaN(correctingFWHM) ? "None" : (correctingFWHM / Unit.arcsec) + ""),
		"Fluxes are corrected for sources of this size (FWHM) in arcsec."));
		
		// The number of scans contributing to this image
		cursor.add(new HeaderCard("SCANS", scans.size(), "The number of scans in this composite image."));			
		
	}
	
	
	public void printHeader() {
		header.dumpHeader(System.out);
	}      

	
	public void despike(double significance) {
		double[][] neighbours = {{ 0, 1, 0 }, { 1, 0, 1 }, { 0, 1, 0 }};
		SkyMap diff = (SkyMap) copy();
		diff.convolve(neighbours);
		
		WeightedPoint point = new WeightedPoint();
		WeightedPoint surrounding = new WeightedPoint();
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			point.value = data[i][j];
			point.weight = weight[i][j];
			surrounding.value = diff.data[i][j];
			surrounding.weight = diff.weight[i][j];
			point.subtract(surrounding);
			if(point.significance() > significance) flag[i][j] = 1;			
		}	
	}
	
	@Override
	public String toString() {
		String info = fileName == null ? "\n" : " Image File: " + fileName + ". ->" + "\n\n" + 
			"  [" + sourceName + "]" + "\n" +
			super.toString() + 
			"  Instrument Beam FWHM: " + Util.f2.format(instrument.resolution / Unit.arcsec) + " arcsec." + "\n" +
			"  Applied Smoothing: " + Util.f2.format(smoothFWHM / Unit.arcsec) + " arcsec." + "\n" +
			"  Image Resolution (FWHM): " + Util.f2.format(getImageFWHM() / Unit.arcsec) + " arcsec. (includes smoothing)" + "\n" +
			"  Noise Estimate from: " + (weightFactor == 1.0 ? "data" : "image (" + Util.f2.format(1.0 / weightFactor) + "x data)") + "\n"; 

		return info;
	}
	
}





