* Created on Mar 10, 2005

*

*/

*

*  contains methods translated from HEALPix Fortran90

*  with increased map resolution in comparison to original Fortran code.

*

* @author N Kuropatkin

* 

*/

  private static double twothird = 2. / 3.;

  private static double PI = Math.PI;

  private static double TWOPI = 2. * PI;

  private static double FOURPI = 4. * PI;

  private static double HALFPI = PI / 2.0;

    private static int ns_max = 1048576; // 2^20

    private static int xmax = 16384;

    private static int pixmax = 131072;

    private static int xmid = 4096;

   

  private static long[] x2pix = new long[xmax+1];

  private static long[] y2pix = new long[xmax+1];

  private static long[] pix2x = new long[pixmax+1];

  private static long[] pix2y = new long[pixmax+1];

  private Vector3d pixVect = new Vector3d(0., 0., 0.);

  private double[][] pixVertex = new double[3][4];

  private BitManipulation bm = new BitManipulation();

  /**

   * defauld constructor

   *

   * 

   */

  public PixTools() {

    for (int i = 0; i <= xmax; i++) {

      x2pix[i] = 0;

      y2pix[i] = 0;

    }

    for (int i = 0; i <= pixmax; i++) {

      pix2x[i] = 0;

      pix2y[i] = 0;

    }

  }

  /**

   * finds pixels having a colatitude (measured from North pole) : theta1 <

   * colatitude < theta2 with o <= theta1 < theta2 <= Pi if theta2 < theta1

   * then pixels with 0 <= colatitude < theta2 or theta1 < colatitude < Pi are

   * returned

   *

   * @param nside

   *            long the map resolution parameter

   * @param theta1

   *            lower edge of the colatitude

   * @param theta2

   *            upper adge of the colatitude

   * @param nest

   *            long if = 1 result is in NESTED scheme

   * @return  ArrayList of  pixel numbers (long)

   * @throws IllegalArgumentException

   */

  public ArrayList query_strip(long nside, double theta1, double theta2,

      long nest) throws Exception {

    ArrayList res = new ArrayList();

    ArrayList listir = new ArrayList();

    long npix, nstrip;

    long iz, ip,  irmin, irmax;

        int is;

    double phi0, dphi;

    double[] colrange = new double[4];

    boolean nest_flag = false;

    String SID = " QUERY_STRIP";

    /* ---------------------------------------- */

    npix = Nside2Npix(nside);

    if (nest == 1)

      nest_flag = true;

    if (npix < 0) {

      throw new IllegalArgumentException(SID + " Nside should be power of 2");

    }

    if ((theta1 < 0.0 || theta1 > PI) || (theta2 < 0.0 || theta2 > PI)) {

      throw new IllegalArgumentException(SID + " Illegal value of theta1, theta2");

    }

    if (theta1 <= theta2) {

      nstrip = 1;

      colrange[0] = theta1;

      colrange[1] = theta2;

    } else {

      nstrip = 2;

      colrange[0] = 0.0;

      colrange[1] = theta2;

      colrange[2] = theta1;

      colrange[3] = PI;

    }

    /* loops on strips */

    for (is = 0; is < nstrip; is++) {

      irmin = RingNum(nside, Math.cos(colrange[2 * is]));

      irmax = RingNum(nside, Math.cos(colrange[2 * is + 1]));

      /* loop on ring number */

      for (iz = irmin; iz <= irmax; iz++) {

        phi0 = 0.;

        dphi = PI;

        listir = InRing(nside, iz, phi0, dphi, nest_flag);

        res.addAll(listir);

      }

    }

    return res;

  }

  /**

   * finds pixels that lie within a CONVEX polygon defined by its vertex on

   * sphere

   *

   * @param nside

   *            the map resolution

   * @param vlist

   *            ArrayList of vectors defining the polygin vertices

   * @param nest

   *            if set to 1 use NESTED scheme

   * @param inclusive

   *            if set 1 returns all pixels crossed by polygon boundaries

   * @return  ArrayList of pixels

   *

   * algorithm: the polygon is divided into triangles vertex 0 belongs to all

   * triangles

   * @throws IllegalArgumentException

   */

  public ArrayList query_polygon(long nside, ArrayList vlist, long nest,

      long inclusive) throws Exception {

    ArrayList res = new ArrayList();

    int nv = vlist.size();

    Vector3d vp0, vp1, vp2;

    Vector3d vo;

    ArrayList vvlist = new ArrayList();

    double surface, fsky, hand;

    double[] ss = new double[nv];

    int n_in_trg, ilist, ntl;

        long npix;

    int ix = 0;

   

    int n_remain, np, nm, nlow;

    String SID = "QUERY_POLYGON";

    //    System.out.println("Start polygon");

    for (int k = 0; k < nv; k++)

      ss[k] = 0.;

    /* -------------------------------------- */

    n_remain = nv;

    if (n_remain < 3) {

      throw new IllegalArgumentException(SID + " Number of vertices should be >= 3");

    }

    /*---------------------------------------------------------------- */

    /* Check that the poligon is convex or has only one concave vertex */

    /*---------------------------------------------------------------- */

    int i0 = 0;

    int i2 = 0;

    if (n_remain > 3) { // a triangle is always convex

      for (int i1 = 1; i1 <= n_remain - 1; i1++) { // in [0,n_remain-1]

        i0 = (int) bm.MODULO(i1 - 1, n_remain);

        i2 = (int) bm.MODULO(i1 + 1, n_remain);

        vp0 = (Vector3d) vlist.get(i0); // select vertices by 3

                        // neighbour

        vp1 = (Vector3d) vlist.get(i1);

        vp2 = (Vector3d) vlist.get(i2);

        // computes handedness (v0 x v2) . v1 for each vertex v1

        vo = new Vector3d(crossProduct(vp0, vp2));

        hand = dotProduct(vo, vp1);

        if (hand >= 0.) {

          ss[i1] = 1.0;

        } else {

          ss[i1] = -1.0;

        }

      }

      np = 0; // number of vert. with positive handedness

      for (int i = 0; i < nv; i++) {

        if (ss[i] > 0.)

          np++;

      }

      nm = n_remain - np;

      nlow = Math.min(np, nm);

      if (nlow != 0) {

        if (nlow == 1) { // only one concave vertex

          if (np == 1) { // ix index of the vertex in the list

            for (int k = 0; k < nv - 1; k++) {

              if (Math.abs(ss[k] - 1.0) <= 1.e-12) {

                ix = k;

                break;

              }

            }

          } else {

            for (int k = 0; k < nv - 1; k++) {

              if (Math.abs(ss[k] + 1.0) <= 1.e-12) {

                ix = k;

                break;

              }

            }

          }

          // rotate pixel list to put that vertex in #0

          int n_rot = vlist.size() - ix;

          int ilast = vlist.size() - 1;

          for (int k = 0; k < n_rot; k++) {

            Vector3d temp = new Vector3d((Vector3d) vlist

                .get(ilast));

            vlist.remove(ilast);

            vlist.add(0, (Object) temp);

          }

        }

        if (nlow > 1) { // more than 1concave vertex

          System.out

              .println(" The polygon has more than one concave vertex");

          System.out.println(" The result is unpredictable");

        }

      }

    }

    /* fill the poligon, one triangle at a time */

    npix = (long) Nside2Npix(nside);

    while (n_remain >= 3) {

      vp0 = (Vector3d) vlist.get(0);

      vp1 = (Vector3d) vlist.get(n_remain - 2);

      vp2 = (Vector3d) vlist.get(n_remain - 1);

      /* find pixels within the triangle */

      ArrayList templist = new ArrayList();

      templist = query_triangle(nside, vp0, vp1, vp2, nest, inclusive);

      vvlist.addAll(templist);

      n_remain--;

    }

    /* make final pixel list */

    npix = vvlist.size();

    long[] pixels = new long[(int)npix];

    for (int i = 0; i < npix; i++) {

      pixels[i] = ((Long) vvlist.get(i)).longValue();

    }

    Arrays.sort(pixels);

    int k = 0;

    res.add(k, new Long(pixels[0]));

    for (int i = 1; i < pixels.length; i++) {

      if (pixels[i] > pixels[i - 1]) {

        k++;

        res.add(k, new Long(pixels[i]));

      }

    }

    return res;

  }

  /**

   * generates a list of pixels that lie inside a triangle defined by

   * the three vertex vectors

   *

   * @param nside

   *            long map resolution parameter

   * @param v1

   *            Vector3d defines one vertex of the triangle

   * @param v2

   *            Vector3d anothe vertex

   * @param v3

   *            Vector3d yet another one

   * @param nest

   *            long 0 (default) RING numbering scheme, if set to 1 the NESTED

   *            scheme will be used.

   * @param inclusive

   *            long 0 (default) only pixels whose centers are inside the

   *            triangle will be listed, if set to 1 all pixels overlaping the

   *            triangle will be listed

   * @return ArrayList with pixel numbers

   * @throws healpixException if the triangle is degenerated

   * @throws IllegalArgumentException

   */

  public ArrayList query_triangle(long nside, Vector3d v1, Vector3d v2,

      Vector3d v3, long nest, long inclusive) throws Exception {

    ArrayList res;

    res = new ArrayList();

    ArrayList listir;

    long npix, ilist, iz, irmin, irmax, n12, n123a, n123b, ndom = 0;

    boolean test1, test2, test3;

    boolean test1a, test1b, test2a, test2b, test3a, test3b;

    double dth1, dth2, determ, sdet;

    double zmax, zmin, z1max, z1min, z2max, z2min, z3max, z3min;

    double z, zz, tgth, st, offset, sin_off;

    double phi_pos, phi_neg;

    Vector3d[] vv = new Vector3d[3];

    Vector3d[] vo = new Vector3d[3];

    double[] sprod = new double[3];

    double[] sto = new double[3];

    double[] phi0i = new double[3];

    double[] tgthi = new double[3];

    double[] dc = new double[3];

    double[][] dom = new double[3][2];

    double[] dom12 = new double[4];

    double[] dom123a = new double[4];

    double[] dom123b = new double[4];

    double[] alldom = new double[6];