code: /*******************************************************************************
   * Copyright (c) 2013 Attila Kovacs <attila[AT]sigmyne.com>.
   * 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[AT]sigmyne.com> - initial API and implementation
   ******************************************************************************/
  package crush.tools;
  
  import java.io.IOException;
  
  import jnum.Unit;
  import jnum.Util;
  import jnum.data.*;
  
  
  public class SimpleDeboosting {
  	static String fileName = "/home/pumukli/data/gismo/counts/dn_ds_beth.txt";
  	//static String fileName = null;
  	
  	int order = 2;
  	SimpleInterpolator table; 
  	
  	public static void main(String[] args) {
  
  		SimpleDeboosting debooster = new SimpleDeboosting();		
  		
  		if(fileName != null) {
  			try { 
  				debooster.table = new SimpleInterpolator(fileName); 
  				
  				for(Interpolator.Data point : debooster.table) {
  					point.ordinate /= 1e-3;
  					point.value *= 1e-3 * Unit.sqdeg / Unit.rad2;
  					point.value = Math.log(point.value);
  					
  					//System.err.println("> " + Util.f3.format(point.ordinate) + " : " + Util.f3.format(point.value));
  				}
  					
  				System.err.println("### Got " + debooster.table.size() + " entries.");
  			}
  			catch(IOException e) { e.printStackTrace(); }
  		}
  		
  		double S1 = Double.parseDouble(args[0]);
  		double sigma = Double.parseDouble(args[1]);			
  		double beamFWHM = 17.5 * Unit.arcsec;
  		//double beamFWHM = 2.0 * Unit.arcsec;
  		
  		System.err.println(" Deboosted: " + debooster.deboost(S1, sigma, beamFWHM).toString(Util.s3));
  	}
  	
  
  	public WeightedPoint deboost(double Sobs, double sigma, double beamFWHM) {		
  		double suma = 0.0, suma2 = 0.0, sumb = 0.0;
  		dS = 0.03 * sigma;
  		double A = -0.5 / (sigma * sigma);
  		double beamArea = 1.13 * beamFWHM * beamFWHM;	// beam-smoothed map
  		double nBeams = Unit.deg2 / beamArea;
  		double confusionLimit = nBeams;
  
  		int toN = (int) Math.ceil((Sobs + 4.0 * sigma) / dS);
  		
  		System.err.println(" dS = " + Util.e2.format(dS) + " mJy");
  		
  		// Around the confusion limit, every beam pretty much has sources in it
  		// So the sources make just a continuous bg, and add to the noise...
  		// Around the confusion limit, every beam pretty much has sources in it
  		// So the sources make just a continuous bg, and add to the noise...
  		int fromN = 1;
  		while(dN(fromN * dS) <= 0.0) fromN++;
  		
  		System.err.println(" no counts below " + Util.e2.format(fromN * dS) + " mJy");
  			
  		// First calculate probability that a beam
  		// has an individual unconfused source in it at each level...
  		double[] C = new double[toN + 1];
  		int maxN = toN;
  		while(dN(maxN * dS) * (maxN*dS) > 0.001 * confusionLimit) maxN *= 2;
  			
  		System.err.println("max: " + Util.e3.format(maxN * dS) + " (" + maxN + ")");
  		
  		boolean isConfused = false;
  		double sumN = 0.0;
  		for(int i=maxN; i >= fromN; i--) {
  			sumN += dN(i * dS) * dS;
  			if(!isConfused) if(sumN > confusionLimit) {
  				System.err.println(" confused below " + Util.e2.format(i * dS) + " mJy");
  				isConfused = true;
  			}
  			
  			if(i <= toN) {
  				C[i] = Math.exp(-sumN / confusionLimit);
  				//System.err.println(Util.f5.format(i*dS) + ": " + Util.f5.format(C[i]));
  			}
  		}
  		
  		
  		for( ; fromN < toN; fromN++) if(C[fromN] > 0.001) break;
  		System.err.println(" completely confused below " + Util.e3.format(fromN * dS));
  		
  		
  		if(C[toN] < 0.1) System.err.println("WARNING! Counts yield no unconfused sources for input flux.");	
  		
  		
  		
  		
  		double[] P = new double[(order+1) * (toN+1)];
  		
  			
  		
  		for(int i=toN; i>=fromN; i--) {
  			double Sa = i * dS;
  			double dN0a = dN(Sa) * dS;
  			
  			// Near the confusion limit sources appear as fewer individual ones on top of a noisy background...
  			double pa = C[i] * dN0a / nBeams; // probability of a beam with individual source.
  					
  			
  			double D = Sa - Sobs;
  			double ea = Math.exp(A*D*D);
  			double Pa = ea * pa;
  			suma += Pa * Sa; 
  			suma2 += Pa * Sa*Sa;
  			sumb += Pa;				// This is the probability that the brightest component is has brightness Sa
  
  			P[i] += Pa;
  			
  			//System.out.println(Util.e3.format(Sa) + "\t" + Util.e3.format(Pa));
  			
  			// Add overlap fluxes...
  			// Do not double count by counting just fainter
  			
  			
  			if(order > 0) for(int j=i; j>=fromN; j--) {
  				double Sb = j * dS;
  				double dN0b = dN(Sb) * dS;
  				double pb = C[j] * dN0b / nBeams; // probability of a beam with individual source.
  				
  
  				double S2 = Sa + Sb;
  				D = S2 - Sobs;
  				double eb = Math.exp(A*D*D);
  				double Pb = (Pa * pb * eb); // Probability of flux arising from Sa+Sb and not just Sa
  				
  				// (Pa - Pb) * Sa + Pb (Sa + Sb) = [Pa Sa] + Pb Sb
  				suma += Pb * Sb;
  				
  				// (Pa - Pb) * Sa^2 + Pb * (Sa+Sb)^2 = [Pa Sa^2] + 2 Pb Sa Sb + Pb Sb^2
  				// = [Pa Sa^2] + Pb Sb (2 Sa + Sb) = [Pa Sa^2] + Pb Sb (Sa + S2)
  				suma2 += Pb * Sb * (Sa + S2); 
  				
  				P[i] -= Pb;
  				P[i+j] += Pb;
  						
  				
  				// Add double overlap.
  				// Do not double count by counting just fainter
  				if(order > 1) for(int k=j; k>=fromN; k--) {
  					double Sc = k*dS;
  					double dN0c = dN(Sc) * dS;
  					double pc = C[k] * dN0c / nBeams; // probability of a beam with individual source.
  
  					double S3 = S2 + Sc;
  					D = S3 - Sobs;
  					double ec = Math.exp(A*D*D);
  					double Pc = (Pb * pc * ec); // probability of flux arising from Sa+Sb+Sc and not just Sa+Sb
  					
  					suma += Pc * Sc;
  					suma2 += Pc * Sc * (S2 + S3);
  					
  					P[i+j] -= Pc;
  					P[i+j+k] += Pc;
  				}
  				
  			}	
  		}
  
  		double mean = suma / sumb;
  		double var = suma2 / sumb - mean*mean;	
  			
  		WeightedPoint result = new WeightedPoint(mean, 1.0/var);
  			
  		double sumP = 0.0;
  		for(int i=1; i<P.length; i++) sumP += P[i];
  		double norm = 1.0 / (sumP * dS);
  		
  		DataPoint S0 = new DataPoint(Sobs, sigma);
  		
  		System.out.println("# observed: " + S0);
  		System.out.println("# deboosted: " + result.toString(Util.getDecimalFormat(DataPoint.significanceOf(result))));
  		System.out.println("# counts: " + fileName);
  		System.out.println("# overlap order: " + order);
  		System.out.println("#");
  		System.out.println("# Si\tp");
  		
  		for(int i=1; i<P.length; i++) {
  			System.out.println(Util.e3.format(i*dS) + "\t" + Util.e3.format(norm * P[i]));
  		}
  		
  		
  		return result;
  	}
  	
  	double dS;
  
  	
  	// 2-mm to 850um ratio = 10.56
  	// (assuming T/(1+z) = 10K, beta=1.5)
  	//double fluxRatio = 1.0 / 10.56;
  
  	// 1.2mm to 850um ratio = 1.75
  	//double fluxRatio = 1.0 / 1.75;
  
  	// in band...
  	double fluxRatio = 1.0;
  
  
  	public double dN(double S) {
  		S /= fluxRatio;
  
  		if(table != null) {
  			
  			// The table is per rad^2 per Jy vs. Jy fluxes.
  			try { return Math.exp(table.getValue(S)); }
  			catch(Exception e) { 
  				System.err.print("!");
  				
  				// Lapi extension to below 100 uJy.
  				//if(S < 0.1) return 3.16e6 * Math.pow(S, -2.262) * Unit.sqdeg / Unit.rad2;
  				
  				
  				
  				//System.err.println(Util.e3.format(S));
  				//System.exit(1);
  				
  				return 0.0; 
  			}
  		}
  		
  		// LESS:
  		//return S < 0.5 ? dN(0.5*fluxRatio) : 150.0 * Math.pow(S/5.0, -3.2);
  
  		// SHADES single:
  		//else return S < 0.5 ? dN(0.5*fluxRatio) : 0.5*(189.0+136.0) * Math.pow(S/5.0, -2.95);
  
  		// SHADES broken:
  		return S < 0.5 ? dN(0.5*fluxRatio) : (735.0/9.0) / (Math.pow(S/9.0, 2.0) + Math.pow(S/9.0, 5.8));
  		
  		//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 : 16.0 * Math.pow(S, -2.8);
  	}
  
  }
  