/*******************************************************************************
 * 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.sourcemodel;

import java.text.*;
import java.util.*;

import util.*;


public class GaussianSource extends CircularRegion {
	public DataPoint peak, FWHM;
	public boolean isCorrected = false;
	
	public GaussianSource(SkyMap map) { super(map); }

	public GaussianSource(SkyMap map, String line) throws ParseException {
		super(map, line);
	}
	
	public GaussianSource(SkyMap map, SphericalCoordinates coords, double r, boolean tweak) {
		super(map, coords, r);
		if(tweak) moveToPeak();
	}
	
	@Override
	public DataPoint finetunePeak() {
		Image rms = map.getRMSImage();
		peak = super.finetunePeak();
		if(peak == null) {
			peak.value = map.valueAtIndex(center);
			peak.setRMS(rms.valueAtIndex(center));
		}
		return peak;
	}
	
	public void scale(double factor) {
		peak.scale(factor);		
	}
	
	public void getChi2(WeightedPoint chi2, double level) {
		chi2.noData();
		getLimits();
		final double sigma = radius.value / Util.sigmasInFWHM / map.delta;
		final double A = -0.5 / (sigma*sigma);

		for(int i=fromi; i<=toi; i++) for(int j=fromj; j<=toj; j++) if(map.flag[i][j] == 0) {
			final double di = i-center.x;
			final double dj = j-center.y;
			final double dev = (map.data[i][j] - level - peak.value * Math.exp(A*(di*di+dj*dj))) / map.getRMS(i, j);
			chi2.value += dev * dev;
			chi2.weight += 1.0;
		}
	}
	
	public void addGaussian(double FWHM, double scaling) {	
		getLimits(3.0, FWHM);
		final double sigma = radius.value / Util.sigmasInFWHM / map.delta;
		final double A = -0.5 / (sigma*sigma);
			
		for(int i=fromi; i<=toi; i++) for(int j=fromj; j<=toj; j++) {
			final double di = i-center.x;
			final double dj = j-center.y;
			map.data[i][j] += scaling * peak.value * Math.exp(A*(di*di+dj*dj));
		}
		
	}
	
	public void add() { add(null); }
	
	public void add(Collection<Region> others) { add(radius.value, 1.0, others); }
	
	public void addPoint() { addPoint(null); }
	
	public void addPoint(Collection<Region> others) { add(map.instrument.resolution, 1.0, others); }	
	
	public void subtract() { subtract(null); }
	
	public void subtract(Collection<Region> others) { add(radius.value, -1.0, others); }
	
	public void subtractPoint() { subtractPoint(null); }
	
	public void subtractPoint(Collection<Region> others) { add(map.instrument.resolution, -1.0, others); }
	
	
	public void add(final double FWHM, final double scaling) {
		add(FWHM, scaling, null);
	}
	
	public void add(final double FWHM, final double scaling, final Collection<Region> others) {			
		double filterFWHM = Double.NaN;
		double imageFWHM = map.getImageFWHM();
	
		if(!Double.isNaN(map.extFilterFWHM)) 
			filterFWHM = Math.sqrt(
					map.extFilterFWHM * map.extFilterFWHM 
					+ FWHM * FWHM 
					- imageFWHM * imageFWHM
			);		
		
		// Remove source only after the filtering corrections are applied to its flux...
		addGaussian(FWHM, scaling);

		// Correct for filtering.
		// Consider that only the tip of the source might escape the filter...
		
		double filterBeamScale = Math.pow(FWHM / filterFWHM, 2.0);
		
		if(!Double.isNaN(map.extFilterFWHM)) {
			filterBeamScale *= Double.isNaN(map.filterBlanking) ? 1.0 : Math.min(1.0, map.filterBlanking / peak.significance());
			addGaussian(filterFWHM, -scaling * filterBeamScale);
		}

		final double sigma = radius.value / Util.sigmasInFWHM / map.delta;
		final double A = -0.5 / (sigma*sigma);		

		// Adjust prior detections.for the filtering around this one.
		if(others != null) for(Region region : others) if(region instanceof GaussianSource) if(region != this) {
			GaussianSource source = (GaussianSource) region;
			final double di = source.center.x - center.x;
			final double dj = source.center.y - center.y;
			
			if(!Double.isNaN(map.extFilterFWHM))
				source.peak.value -= scaling * filterBeamScale * peak.value * Math.exp(A*(di*di+dj*dj));
		}
		
	}
	
	public double getCorrectionFactor(double FWHM) {	
		double correction = 1.0;	
		
		// Correct for filtering.
		// Consider that only the tip of the source might escape the filter...
		if(!Double.isNaN(map.extFilterFWHM)) {
			double filterFraction = Double.isNaN(map.filterBlanking) ? 1.0 : Math.min(1.0, map.filterBlanking / peak.significance());
			double filtering = 1.0 - 1.0 / map.getExtFilterCorrectionFactor(FWHM);;
			correction *= 1.0 / (1.0 - filtering * filterFraction);
		}
		
		return correction;
	}
	
	
	public void correct(double FWHM) {	
		if(isCorrected) throw new IllegalStateException("Source is already corrected.");
		double correction = getCorrectionFactor(FWHM);
		peak.scale(correction);
		isCorrected = true;
	}
	

	public void uncorrect(double FWHM) {
		if(!isCorrected) throw new IllegalStateException("Source is already uncorrected.");
		double correction = getCorrectionFactor(FWHM);
		peak.scale(1.0 / correction);
		isCorrected = false;
	}
	
	
	// Formula from Kovacs et al. (2006)
	public void setSearchRadius(double pointingRMS) {
		double beamSigma = map.instrument.resolution / Util.sigmasInFWHM;
		radius.value = Math.sqrt(4.0 * pointingRMS * pointingRMS - 2.0 * beamSigma * beamSigma * Math.log(1.0 - 2.0 / peak.significance()));
	}

	
	@Override
	public void parse(String line) throws ParseException {
		if(line == null) return;
		if(line.length() == 0) return;
		if(line.charAt(0) == '#' || line.charAt(0) == '!') return;

		StringTokenizer tokens = new StringTokenizer(line);

		if(tokens.countTokens() < 5) return;

		id = tokens.nextToken();

		EquatorialCoordinates coords = new EquatorialCoordinates();
		coords.parse(tokens.nextToken() + " " + tokens.nextToken() + " " + tokens.nextToken());
		//System.err.println(coords.toString());
		if(map.projection.getReference() instanceof Precessing) {
			CoordinateEpoch epoch = ((Precessing) map.projection.getReference()).getEpoch();
			coords.precess(epoch);
		}
			
		Vector2D offset = new Vector2D();
		map.projection.project(coords, offset);
		center.x = map.indexOfdX(offset.x);
		center.y = map.indexOfdY(offset.y);

		radius.value = Double.parseDouble(tokens.nextToken()) * Unit.arcsec;
		if(tokens.hasMoreTokens()) peak.value = Double.parseDouble(tokens.nextToken());
		
		String nextArg = tokens.hasMoreTokens() ? tokens.nextToken() : null;
		if(nextArg != null) {
			if(nextArg.equals("+-")) nextArg = tokens.nextToken();
			try { 
				peak.setRMS(Double.parseDouble(nextArg));
				nextArg = tokens.hasMoreTokens() ? tokens.nextToken() : null;
			}
			catch(NumberFormatException e) {}
		}
		
		Unit unit = nextArg == null ? map.janskyPerBeam : map.getUnit(nextArg);
		while(tokens.hasMoreTokens()) comment += tokens.nextToken() + " ";
		comment.trim();
		
		peak.scale(unit.value);
	}

	@Override
	public String toCrushString() {
		return id + "\t" + super.toCrushString() + "  " + getPeakString(map.unit);
	}
	
	@Override
	public String getComment() {
		return "s/n=" + Util.f2.format(peak.significance()) + " " + super.getComment();
	}
	
	@Override
	public StringTokenizer parseCrush(String line) {
		StringTokenizer tokens = super.parseCrush(line);
		
		peak.value = Double.parseDouble(tokens.nextToken());
		String next = tokens.nextToken();
		if(next.equals("+-")) {
			peak.setRMS(Double.parseDouble(tokens.nextToken()));
			next = tokens.nextToken();
		}
		
		peak.value *= map.getUnit(next).value;
		
		return tokens;
	}
	
	public String getPeakString(Unit unit) {
		return Util.e3.format(peak.value / unit.value) + " +- " + Util.e2.format(peak.rms() / unit.value) + " " + unit.name;
	}
	
	public Vector2D centroid() {
		Vector2D centroid = new Vector2D();
		double sumw = 0.0;
		
		for(int i=fromi; i<=toi; i++) for(int j=fromj; j<=toj; j++) if(map.flag[i][j] == 0) {
			double w = map.getS2N(i,j);
			centroid.x += w * i;
			centroid.y += w * j;
			sumw += Math.abs(w);
		}
		centroid.scale(1.0/sumw);
		return centroid;
	}
	
	public double spread() {
		double sumd2 = 0.0, sumw = 0.0;
		
		for(int i=fromi; i<=toi; i++) for(int j=fromj; j<=toj; j++) if(map.flag[i][i] == 0) {
			double w = map.getS2N(i,j);
			double d = (i - center.x);
			sumd2 += w * d * d;
			d = (j - center.y);
			sumd2 += w * d * d;
			sumw += Math.abs(w);
		}
		
		return Math.sqrt(sumd2 / sumw);
	}
	
}
