/*******************************************************************************
 * Copyright (c) 2013 Attila Kovacs <attila_kovacs[AT]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[AT]post.harvard.edu> - initial API and implementation
 ******************************************************************************/
package crush.devel;

import jnum.Unit;
import jnum.math.specialfunctions.CumulativeNormalDistribution;

public class DetectCountEstimator {

	double fluxRatio = 1.0;
	
	public static void main(String[] args) {
		double area = 20 * Unit.arcmin * 20 * Unit.arcmin;
		double rms = 0.35; // mJy
		double from = 1.4; // mJy
		double to = 2.0; // mJy
		
		DetectCountEstimator estimator = new DetectCountEstimator();
		System.err.println(estimator.getCount(area, rms, from, to) / Unit.deg2);
		
	}
	
	// Flux limit in mJy, area in normal units.
	public double getCount(double area, double rms, double fluxlimit, double upperlim) {
		double dF = 0.1 * rms;
		
		double sum = 0.0;
		
		for(double F=0.0; F<upperlim; F += dF) {
			double dev = (fluxlimit - F) / rms;
			double completeness = CumulativeNormalDistribution.inverseComplementAt(dev);
			if(dev < 0.0) completeness = 1.0 - completeness;
			//System.err.println(("c(" + Util.f3.format(F / rms)) + ") = " + Util.f3.format(completeness));
			
			sum += completeness * dN(F);
		}
		
		return area * dF * sum;
	}
	
	public double dN(double S) {
		S /= fluxRatio;
		
		// LESS:
		//return S < 0.5 ? 0.0 : 150.0 * Math.pow(S/5.0, -3.2);
		
		// SHADES:
		//return S < 0.5 ? 0.0 : 0.5*(189.0+136.0)*Math.pow(S/5.0, -2.95);
		
		//return 1600*Math.pow(S/3.3, -1.0) * Math.exp(-S/3.3);		
		//return 382.7 / (0.005 + Math.pow(S, 4.30));	
		//return 2800 * Math.pow(S, -2.39) * Math.exp(-0.62 * S);
		
		
		// 2-mm counts from Eli, (straightedned by eye)
		return S < 0.01 ? 0.0 : 5 * 16.0 * Math.pow(S, -2.8);
	}

}
