/**
  *  ApPhot extends the StarField class to do simple Aperture Photometry.
  *  It uses an algorithm for peak merging appropriate for uncrowded star fields.
  *  Glen Cowan, RHUL Physics Department
  *  November 2004
  */

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

public class ApPhot extends StarField {
   
// private data fields

  private Vector starVec;
  private Vector apVec;
  private Vector bckApVec;
  private double starThreshold;
  
// constructors

  public ApPhot(){
    super();
    starVec = new Vector();
    apVec = new Vector();
    bckApVec = new Vector();
    starThreshold = 0;
  }

  public ApPhot(int nx, int ny){
    super(nx, ny);
    starVec = new Vector();
    apVec = new Vector();
    bckApVec = new Vector();
    starThreshold = 0;
  }

  public ApPhot(double starThreshold){
    super();
    starVec = new Vector();
    apVec = new Vector();
    bckApVec = new Vector();
    this.starThreshold = starThreshold;
  }

  public ApPhot(int nx, int ny, double starThreshold){
    super(nx, ny);
    starVec = new Vector();
    apVec = new Vector();
    bckApVec = new Vector();
    this.starThreshold = starThreshold;
  }

// public methods

  public Vector getStarVec(){ return starVec; }
  public void setStarVec(Vector starVec){ this.starVec = starVec; }
  public Vector getApVec(){ return apVec; }
  public void setApVec(Vector apVec){ this.apVec = apVec; }
  public Vector getBckApVec(){ return bckApVec; }
  public void setBckApVec(Vector bckApVec){ this.bckApVec = bckApVec; }

// fillStarVec puts the stars into order and inserts them into starVec

  public void fillStarVec(){

    int numClus = clusterVec.size();
    Star[] starArray = new Star[numClus];
    for (int i=0; i<numClus; i++){
      Cluster c = (Cluster)clusterVec.elementAt(i);
      starArray[i] = c.getStar(this);
    }
    Arrays.sort(starArray);

    System.out.println("Sorted star array...");
    for (int i=0; i<numClus; i++){
      System.out.println(i + "  " +  starArray[i].getVal());
    }

// Make threshold cut and pack the stars into starVec

    starVec.clear();
    for (int i=0; i<starArray.length; i++){
      double val = starArray[i].getVal();
      if ( val >= starThreshold ) { starVec.addElement( starArray[i] ); }
    }

    System.out.println ("starVec i, x, y, val...");
    for (int i=0; i<starVec.size(); i++){
      Star s = (Star)starVec.elementAt(i);
      double x   = s.getX();
      double y   = s.getY();
      double val = s.getVal();
      System.out.println(i + "  " + x + "  " + y + "  " + val);
    }

  }                     // end of fillStarVec

  public double getMinDist(){

// returns smallest distance between any two stars

    int numStars = starVec.size();
    double dMin = 0;
    boolean firstPair = true;
    for (int i=0; i<numStars-1; i++){
      Star si = (Star)starVec.elementAt(i);
      for (int j=i+1; j<numStars; j++){
        Star sj = (Star)starVec.elementAt(j);
        double d = starPairDistance(si, sj);      
        if ( firstPair || d < dMin ){
          firstPair = false;
          dMin = d;
        }
      }
    }
    return dMin;
  
  }        // end findMinDist

  static double starPairDistance(Star s1, Star s2){
    double x1 = s1.getX();
    double y1 = s1.getY();
    double x2 = s2.getX();
    double y2 = s2.getY();
    double d = Math.sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
    return d;
  }

  public void fillApVec(double rSig, double rBckInner, double rBckOuter){

// Define circular signal aperture around each star and annular 
// aperture for background.

    apVec.clear();
    bckApVec.clear();
    int numStars = starVec.size();
    for (int i=0; i<numStars; i++){
      Star s = (Star)starVec.elementAt(i);
      double x = s.getX();
      double y = s.getY();
      Aperture aSig = new Aperture(x, y, rSig);       // creates apType = 0
      apVec.addElement(aSig);
      Aperture aBck = new Aperture(x, y, rBckInner, rBckOuter);  // creates apType = 1
      bckApVec.addElement(aBck);
    }
  }       // end of fillApVec

// 1-arg version for signal regions only

  public void fillApVec(double radius){
    this.fillApVec(radius, 0., 0.);
  }

// 2-arg version for background regions only

  public void fillApVec(double rInner, double rOuter){
    this.fillApVec(0., rInner, rOuter);
  }

  public void measureAps(){
    System.out.println("star, sigVal, sigArea, bckVal, bckArea...");
    int numAps = apVec.size();
    for (int i=0; i<numAps; i++){
      Aperture sigAp = (Aperture)apVec.elementAt(i);
      Aperture bckAp = (Aperture)bckApVec.elementAt(i);
      double sigVal  = sigAp.getVal(this);
      int sigArea    = sigAp.getArea(this);
      double bckVal  = bckAp.getVal(this);
      int bckArea    = bckAp.getArea(this);
      System.out.println(i + "  " + sigVal + "  " + sigArea + 
                             "  " + bckVal + "  " + bckArea );
    }
  }

  public void subtractBck(){

// fills nStar and sigmaNStar into the stars in starVec

    System.out.println("star, nStar, sigmaNStar...");
    double gain = this.getGain();
    int numStars = starVec.size();
    for (int i=0; i<numStars; i++){

      Star s         = (Star)starVec.elementAt(i);
      Aperture sigAp = (Aperture)apVec.elementAt(i);
      Aperture bckAp = (Aperture)bckApVec.elementAt(i);

      double sigVal  = sigAp.getVal(this);
      int sigArea    = sigAp.getArea(this);
      double bckVal  = bckAp.getVal(this);
      int bckArea    = bckAp.getArea(this);

      double w       = (double)sigArea / (double)bckArea;
      double Ns      = sigVal * gain;                     // Ns and Nb in units of
      double Nb      = bckVal * gain;                     // photoelectrons
      s.nStar        = Ns - w*Nb;
      s.sigmaNStar   = Math.sqrt(Ns + w*w*Nb);

      System.out.println(i + "  " + s.nStar  + "  " + s.sigmaNStar);
    }

  }         // end of subtractBck

// deltaM and sigmaDeltaM give difference of magnitudes and error on diff.

  public double deltaM(int i, int j){
    Star si = (Star)starVec.elementAt(i);
    Star sj = (Star)starVec.elementAt(j);
    double nStari = si.getNStar();
    double nStarj = sj.getNStar();
    double mij = -2.5 * Math.log(nStari/nStarj) / Math.log(10.) ;
    return mij;    
  }

  public double sigmaDeltaM(int i, int j){
    Star si = (Star)starVec.elementAt(i);
    Star sj = (Star)starVec.elementAt(j);
    double nStari = si.getNStar();
    double sigmaNStari = si.getSigmaNStar();
    double nStarj = sj.getNStar();
    double sigmaNStarj = sj.getSigmaNStar();
    double ri = sigmaNStari / nStari;
    double rj = sigmaNStarj / nStarj;
    double sigma = (2.5 / Math.log(10.)) * Math.sqrt(ri*ri + rj*rj);
    return sigma;    
  }

}         // end of class ApPhot

//-------------------------------------------------------------------------

class Aperture extends Annulus {

  private double val;         // in ADU
  private int area;           // in pixels
  private int apType;         // 0 for signal, 1 for background

  Aperture(){
    super();
    this.val = -1;
    this.area = -1;
    this.apType = 0;
  }

  Aperture(double x, double y, double radius){
    super(x, y, radius);
    this.val = -1;
    this.area = -1;
    this.apType = 0;
  }

  Aperture(double x, double y, double rInner, double rOuter){
    super(x, y, rInner, rOuter);
    this.val = -1;
    this.area = -1;
    this.apType = 1;
  }

  public int getApType(){ return this.apType; }
  public void setApType(int apType){ this.apType = apType; }

  public double getVal(ApPhot ap){
    if ( this.val < 0 ) {
      this.measure(ap);
    }
    return this.val;
  }

  public int getArea(ApPhot ap){
    if ( this.area < 0 ) {
      this.measure(ap);
    }
    return this.area;
  }

  void measure(ApPhot ap){

    int nx = ap.getNx();
    int ny = ap.getNy();
    Pixel[][] pixelArray = ap.getPixelArray();
    
// Loop over pixels, sum if contained in aperture, but exclude if pixel is
// within the outer radius of any other aperture of the same apType.
// For now, loop is over all pixels.  Can be made faster to limit sum over a
// rectangle that bounds the aperture.

    Vector apVec = null;
    if ( this.getApType() == 0 ) {
      apVec = ap.getApVec();
    }
    else {
      apVec = ap.getBckApVec();
    }

    int numStars = apVec.size();
    double valSum = 0.;
    int areaSum = 0;
    for(int i=0; i<nx; i++){
      for (int j=0; j<ny; j++){

        Pixel p = pixelArray[i][j];

        boolean excludePixel = false, inOtherAp = false;
        for (int k=0; k<numStars; k++){
          Aperture a = (Aperture)apVec.elementAt(k);
          inOtherAp = a.inOuterAperture(p) && !(this == a);
          if ( inOtherAp ) {
             excludePixel = true;
          }
        }
        if ( this.inAperture(p) && !excludePixel ){
          valSum += (double)p.getVal();
          areaSum++;
        }

      } 
    }

    this.val = valSum;
    this.area = areaSum;

  }      // end of measure 

  boolean inAperture(Pixel p){

    final double epsilon = 1e-6;       // rInner any smaller and it's a circle

    double xAp = this.getX();
    double yAp = this.getY();
    double rInner = this.getRInner();
    double rOuter = this.getROuter();

    double xPix = (double)p.getI();
    double yPix = (double)p.getJ();

    double distSquared = (xAp-xPix)*(xAp-xPix) + (yAp-yPix)*(yAp-yPix);
    boolean insideROuter = distSquared <= rOuter*rOuter;
    boolean outsideRInner = true;
    if ( rInner > epsilon ) {
      outsideRInner = distSquared >= rInner*rInner;
    }
    boolean inAp = outsideRInner && insideROuter;

    return inAp;
  }

  boolean inOuterAperture(Pixel p){
    double xAp = this.getX();
    double yAp = this.getY();
    double rOuter = this.getROuter();
    double xPix = (double)p.getI();
    double yPix = (double)p.getJ();
    double distSquared = (xAp-xPix)*(xAp-xPix) + (yAp-yPix)*(yAp-yPix);
    boolean inAp = distSquared <= rOuter*rOuter;
    return inAp;
  }


}             //  end of class Aperture