import java.util.*;
import cds.aladin.*;
import cds.astro.*;
import cds.tools.Util;

final public class ProbabiliteIdentification extends AladinPlugin {
   
   public String menu()   { return "X-match probability"; }
   public String author() { return "Pierre Fernique [CDS]"; }
   public String description() {
      return "This plugin has been developped for Hight Energy Team demonstration\n\n" +
            "It determines a probability of a xmatch function.\n" +
            "1) The user does a xmatch between two catalogs\n" +
            "2) The user have to select the most dense catalog plane and launch the plugin. " +
            "It will determine the local density of brighter stars around each xtamched object" +
            "and thus, computes the probability of this xmatch based on the function:\n" +
            "       Proba = 1 - sigma*error*exp(-sigma*error)";
   }
   public String category() { return "Prototype"; }
   public String version()  { return "2.0 - January 2007"; }
   public String url()      { return "http://aladin.u-strasbg.fr/java/Plugins/CatalogCreationPlug.java"; }
   public boolean inSeparatedThread() { return true; }
               
   static final int MAX = 200;
   public void exec() {
      
      try {
         AladinData xPlan = aladin.getAladinData("XMatch*");
         AladinData gPlan = aladin.getAladinData("GSC*");
         int xN           = xPlan.getNbObj();
         int gN           = gPlan.getNbObj();
         Obj [] xObj      = xPlan.seeObj();
         Obj [] gObj      = gPlan.seeObj();
         double magX []   = findXmatchMag(xObj,xN);         // Mags du cat. xmatch
         int bMag         = gObj[0].indexOf("phot.mag*;*B"); // Indice colonne des Mag du catalogue
         int rMag         = gObj[0].indexOf("phot.mag*;*R"); // Indice colonne des Mag du catalogue
         int iPosErr      = xObj[0].indexOf("PosErr*");      // Indice colonne des PosErr dans Xmatch
         
         if( bMag==-1 || rMag==-1 || iPosErr==-1 ) throw new Exception("Bad column indexes !");

         // Indexation Qbox
         int level = 8;
         Qbox.setLevel(level);
         int size = (int)Math.pow(4, level) * 6;
         ArrayList al [] = new ArrayList[size];

         for( int i=0; i<gN; i++ ) {
            Qbox q = new Qbox(new Coo(gObj[i].raj,gObj[i].dej));
            int index = q.index();
            if( al[index]==null ) al[index] = new ArrayList();
            al[index].add(gObj[i]);
         }

         // Ajout de deux nouvelles colonnes pour mémoriser la densité et la proba
         int pos = xObj[0].getSize();
         xObj[0].setColumn(pos,   "sigma", null, null,6);
         xObj[0].setColumn(pos+1, "proba", null, null,6);

         // Recherche des N objets aussi brillants autour de chaque Xmatch
         for( int j=0; j<xN; j++ ) {
            TreeSet ts = new TreeSet();
            boolean ok [] = new boolean[size];

            // On regarde autour en agrandissant peu à peu le cercle
            for( double rayon=0.5; ts.size()<MAX && rayon<10; rayon+=1 ) {
               Enumeration en = Qbox.circle(new Coo(xObj[j].raj,xObj[j].dej),rayon);
               while( en.hasMoreElements()) {
                  Qbox q = (Qbox)en.nextElement();
                  int index = q.index();
                  if( ok[index] ) continue;
                  ok[index]=true;

                  ArrayList a = al[index];
                  if( a==null ) continue;

                  // Ajout dans un arbre binaire borné
                  for( int i=0; i<a.size(); i++ ) {
                     Obj x = (Obj)a.get(i);
                     try { 
                        if( Double.parseDouble( x.getValues()[bMag].trim() )
                              > magX[j] ) continue;
                     } catch( Exception e ) { 
                        try { 
                           if( Double.parseDouble( x.getValues()[rMag].trim() )
                                 > magX[j] ) continue;
                        } catch( Exception e1 ) { continue; }
                     }

                     double dist = x.getDistance(xObj[j]);

                     Case c = new Case(x,dist);
                     if( ts.isEmpty() ) ts.add(c);
                     else {
                        Case last = (Case)ts.last();
                        if( ts.size()==MAX-1 && c.dist>last.dist ) continue;
                        if( c.dist<last.dist ){
                           ts.add(c);
                           if( ts.size()==MAX ) ts.remove(last);
                        }
                     }
                  }
               }
            }

            // Ajout de la valeur résultante dans la colonne additionnelle
            if( !ts.isEmpty() ) {
               try {
                  double r = ((Case)ts.last()).dist;
                  double sigma = ts.size()/Math.PI*r*r;
                  double deltaS = Double.parseDouble(xObj[j].getValues()[iPosErr]);
                  deltaS = deltaS/3600.;
                  double proba = 1 - sigma*deltaS * Math.exp( -sigma*deltaS );

                  xObj[j].setValue(pos,Util.myRound(""+sigma,5));
                  xObj[j].setValue(pos+1,Util.myRound(""+proba,5));

                  xObj[j].setSelected(true);
                  xPlan.objModified();
               } catch( Exception e ) { e.printStackTrace(); }
            }
         }    

         // Création d'un filtre pour afficher des cercles proportionnels à la proba
         aladin.execCommand("filter f1 { { draw fillcircle(1-${proba}) }};" +
                            "md Folder;" +
                            "mv XMatch* Folder;" +
                            "mv f1 Folder;" +
                            "filter f1 on");
         
      } catch( Exception e ) {
         aladin.warning("Plugin error: "+e.getMessage()+"\n" +
                        "Usage: load GSC2.2, load ROSAT and X-match them\n" +
                        "before launching this plugin.");
         return;
      }
   }
   
   /** Extration des magnitudes GSC des objets X-matchés */
   private double [] findXmatchMag(Obj obj[], int n) {
      double mag [] = new double[n];
      int bMag = obj[0].indexOf("phot.mag*;*B");
      int rMag = obj[0].indexOf("phot.mag*;*R");
      
      for( int i=0; i<n; i++ ) {
         try { mag[i] = Double.parseDouble( obj[i].getValues()[bMag].trim() ); }
         catch( Exception e ) { mag[i]=-1; }
         if( mag[i]==-1 ) 
            try { mag[i] = Double.parseDouble( obj[i].getValues()[rMag].trim() ); }
            catch( Exception e ) { mag[i]=-1; }
      }
      return mag;
   }

   
   final protected class Case implements Comparable {
      protected Obj obj;
      protected double dist;
      Case(Obj o,double d) { obj=o; dist=d; }
      void set(Obj o,double d) { obj=o; dist=d;}
      public int compareTo(Object arg0) {
         return ((Case)arg0).dist == dist ? 0 : ((Case)arg0).dist > dist ? -1 : 1;
      }           
   }
}