package nxm.sys.libg;

import nxm.sys.inc.Constants;
import nxm.sys.lib.Midas;

/** A class to draw Geodetic views.
    <pre>
      Transformation matrix:
        x' = mxx*x + mxb    x = (x'-mxb)/mxx
    </pre>
    @author Jeff Schoen
 */
public class ViewGeo extends View {
  private static final float PI         = (float)Math.PI;              // Constant for fastAtan2
  private static final float PI_LO      = 1.2246467991473531772E-16f;  // Constant for fastAtan2
  private static final float PI_OVER_4  = PI / 4;                      // Constant for fastAtan2
  private static final float PI_OVER_2  = PI / 2;                      // Constant for fastAtan2
  private static final float PI2        = 2 * PI;                      // Constant for fastPix2rwc
  private static final float PI3        = 3 * PI;                      // Constant for fastPix2rwc
  private static final float DEG2RADF   = (float)DEG2RAD;              // Constant for fastPix2rwc

  /** observer latitude */
  public double lat = 0;

  /** observer longitude */
  public double lon = 0;

  /** observer altitude */
  public double alt=6;

  private double[] tsin,tcos,dtsin,dtcos;
  /** sin of negative observer latitude (in radians) */
  private double slat;
  /** cos of negative observer latitude (in radians) */
  private double clat;
  private double zvis;
  private double rvis;
  /** Equatorial radius of the Earth in meters. This is also called the semi-major axis of the Earth. */
  private double equatorialRadius;
  /** inverse of equatorial radius */
  private double inverseEquatorialRad;
  /** observer (relative?) latitude as float*/
  private float rlatO;
  /** observer (relative?) longitude as float*/
  private float rlonO;
  /** observer (relative?) altitude as float*/
  private float raltO;
  /** observer (relative?) longitude as double*/
  private double drlon;
  /** observer (relative?) altitude as double*/
  private double draltO;
  /** number of degrees stored in table lookup */
  private static final int ndeg=725;

  /** Construct a new ViewGeo
      @param MP the MPlot associated with this ViewGeo
   */
  public ViewGeo (MPlot MP) {
    super(MP);
    prepTables();
    setup();
  }

  /** Fill in the sincos lookup tables */
  private void prepTables () {
    tsin = new double[ndeg];
    tcos = new double[ndeg];
    dtsin = new double[ndeg];
    dtcos = new double[ndeg];
    for (int i=0; i<ndeg; i++) {
      tsin[i] = Math.sin( Constants.DEG2RAD * (double)i );
      tcos[i] = Math.cos( Constants.DEG2RAD * (double)i );
    }
    for (int i=0; i<ndeg-1; i++) {
      dtsin[i] = tsin[i+1]-tsin[i];
      dtcos[i] = tcos[i+1]-tcos[i];
    }
  }

  @Override
  public boolean setObserver (double[] arg) {
    lon = arg[0];
    lat = arg[1];
    alt = arg[2];
    while (lat>+180) lat -= 360;
    while (lat<-180) lat += 360;
    while (lon>+180) lon -= 360;
    while (lon<-180) lon += 360;
    setup();
    return true;
  }

  @Override
  public boolean getObserver (double[] arg) {
    arg[0] = lon;
    arg[1] = lat;
    arg[2] = alt;
    return true;
  }

  @Override
  public void setup () {
    rlonO=(float)lon;
    rlatO=(float)lat;
    raltO=(float)alt;
    drlon=lon;
    draltO=alt;
    slat = Math.sin(-lat*DEG2RAD);
    clat = Math.cos(-lat*DEG2RAD);
    zvis = 1.0/alt; // inverse of altitude
    double pang = Math.asin(zvis);
    rvis = (1./(Math.cos(pang)*Math.cos(pang)));
    equatorialRadius  = EARTH_EQUATORIAL_RADIUS;
    inverseEquatorialRad = (1.0/equatorialRadius);
  }

  @Override
  public int rwc2pix (double[] rwc, MPoint pix) {
    while (rwc[0]>+180) rwc[0] -= 360;
    while (rwc[0]<-180) rwc[0] += 360;
    while (rwc[1]>+180) rwc[1] -= 360;
    while (rwc[1]<-180) rwc[1] += 360;
    int vis = rwc2pix (rwc, 0, 3, 1, MP.tpix);
    pix.x = MP.tpix.x[0]; pix.y = MP.tpix.y[0];
    return vis;
  }

  @Override
  public int pix2rwc (MPoint pix, double[] rwc) {
    double p2,rad,fact,retx,rety,retz,rretx,rrety,rretz;
    retx = (((double)pix.x-MP.mxb)/MP.mxx);
    rety = (((double)pix.y-MP.myb)/MP.myy);
    p2 = draltO*draltO;
    rad = retx*retx + rety*rety;
    fact = quadsolve (p2+rad, -2*p2, p2-1);
    if (fact==0) {  // off earth
      rad = Math.atan2(rety,retx);
      retx = Math.cos(rad);
      rrety = Math.sin(rad);
      rety = rrety*clat;
      retz = rrety*slat;
    } else {
      rretx = retx*fact;
      rrety = rety*fact;
      rretz = Math.sqrt(1-rad*fact*fact);
      retx =  rretx;
      rety =  rrety*clat - rretz*slat;
      retz =  rrety*slat + rretz*clat;
    }
    retx = RAD2DEG*Math.atan2(retx,retz);
    rwc[0] = ((drlon+retx+540) % 360) - 180;
    rwc[1] = RAD2DEG*Math.asin(rety);
    rwc[2] = 0;
    return 1;
  }

  private double quadsolve (double a, double b, double c) {
    double d = b*b - 4*a*c;
    if (d<1e-20 || Math.abs(a)<1e-20) return 0; // off earth
    return (-b-Math.sqrt(d))/(2*a);
  }

  @Override
  public int rwc2pix (float[] buf, int off, int dim, int npts, Pix pix, double xstart, double xdelta) {
    // THIS METHOD IS IDENTICAL TO THE ONE THAT USES A float[] BUFFER FOR INPUT,
    // EXCEPT FOR THE USE OF "float"
    int i,j=off,n=pix.setSize(npts),xoff=off/dim;
    float[] buf2 = new float[n*2];
    // Since 4.0.0 switched to i<n to minimize calculations in the loop, but switched to bufOffset inside the loop to
    // keep the same number of buffer data points and loop iterations
    int bufOffset = 0;
    for (i=0; i<n; j+=dim,i++) {
      double x = xstart + (xoff+i)*xdelta; // Since 4.0.0 NXM-3624 modified x-data calculation to avoid precision loss
      buf2[bufOffset++] = (float)Math.IEEEremainder(x,360.0);
      buf2[bufOffset++] = (float)Math.IEEEremainder(buf[j],360.0);
    }
    return rwc2pix(buf2,0,2,n,pix);
  }

  @Override
  public int rwc2pix (double[] buf, int off, int dim, int npts, Pix pix, double xstart, double xdelta) {
    // THIS METHOD IS IDENTICAL TO THE ONE THAT USES A float[] BUFFER FOR INPUT,
    // EXCEPT FOR THE USE OF "double"
    int i,j=off,n=pix.setSize(npts),xoff=off/dim;
    double[] buf2 = new double[n*2];
    // Since 4.0.0 switched to i<n to minimize calculations in the loop, but switched to bufOffset inside the loop to
    // keep the same number of buffer data points and loop iterations
    int bufOffset = 0;
    for (i=0; i<n; j+=dim,i++) {
      double x = xstart + (xoff+i)*xdelta; // Since 4.0.0 NXM-3624 modified x-data calculation to avoid precision loss
      buf2[bufOffset++] = Math.IEEEremainder(x,360.0);
      buf2[bufOffset++] = Math.IEEEremainder(buf[j],360.0);
    }
    return rwc2pix(buf2,0,2,n,pix);
  }

  @Override
  public int rwc2pix (float[] buf, int off, int dim, int npts, Pix pix) {
    int i,ilon,ilat,j=off,k=0; int[] pax=pix.x, pay=pix.y; float[] paz=pix.z;
    double alon,alat,plon,plat,cp,sp,ct,st,ctcp;
    double mxx=MP.mxx, myy=MP.myy, mxb=MP.mxb, myb=MP.myb;
    double xx,yy,zz,ff,rr;
    boolean vis, viso=false, chk=(pix.flags&Pix.CHECK)!=0;
    npts = Math.min(npts,pix.max);
    for (i=0; i<npts; i++) {
      alon = buf[j++] - rlonO + 360.0f;
      alat = buf[j++] + 360.0f;
      ilon = (int)alon;
      ilat = (int)alat;
      plon = alon-(float)ilon;
      plat = alat-(float)ilat;
      cp = tcos[ilon] + plon*dtcos[ilon];
      sp = tsin[ilon] + plon*dtsin[ilon];
      ct = tcos[ilat] + plat*dtcos[ilat];
      st = tsin[ilat] + plat*dtsin[ilat];
      ctcp = ct*cp;
      if (dim==2) {
        xx = ct*sp;
        yy = ctcp*slat + st*clat;
        zz = ctcp*clat - st*slat;
        vis = zz > zvis;
        ff = draltO/(draltO-zz);
        xx = xx*ff;
        yy = yy*ff;
      } else {
        rr = (buf[j++]+equatorialRadius)*inverseEquatorialRad;
        xx = rr*(ct*sp);
        yy = rr*(ctcp*slat + st*clat);
        zz = rr*(ctcp*clat - st*slat);
        ff = draltO/(draltO-zz);
        xx = xx*ff;
        yy = yy*ff;
        vis = zz<draltO && (zz>zvis || xx*xx+yy*yy>=rvis);
      }
      if (vis) {
        pax[k] = (int)(xx*mxx + mxb);
        pay[k] = (int)(yy*myy + myb);
        paz[k] = (float)zz;
        k++;
      }
      else if (chk) { pax[k]=0;  pay[k]=0; paz[k]=(float)zz; k++; }
      else if (viso) break;
      viso = vis;
    }
    pix.off=i;
    pix.n = i;
    if (!chk) k = checkrange(pax,pay,k);
    return k;
  }

  /** {@inheritDoc}
      <pre>
      Note: These calculations are based on Geocentric non-earth flattening coordinates (ECEF/ECR)

      BASE CASE OF OBSERVER AT 0,0,0

      P - point to convert to pixels
      x - axis that points towards 0 degrees latitude and 90 degrees longitude
      y - axis that points towards north pole
      z - axis that points in the direction of 0 degrees latitude 0 degrees longitude
                                           y
                                    P   . .|. .                               .      .
                                     \     |      .
                                   .  \    |        .      &lt;- earth
                                  .    \   |         .
                                 .      \  |          .
                                 .       \ |          .
                                 .        \|-------------- x
                                 .        /           .
                                 .       /            .
                                  .     /            .
                                   .   /            .
                                     ./           .
                                     /   . . . .
                                    z
        r = (altP + equatorialRadius)/equatorialRadius
        x = r*cos(latP)*sin(lonP)
        y = r*sin(latP)
        z = r*cos(latP)*cos(latP)
        where latP, lonP, altP are the latitude, longitude, and altitude of the Point

        VIEW FROM OBSERVER - ROTATION ABOUT X-AXIS
        (Note: the rotation due to Observer longitude is handled by adjusting the Point
        longitude: alon = buf[j++] - lon + 360.0f)

        |x'|   | 1     0        0        | |x|
        |y'| = | 0  cos(-latO) sin(latO) | |y|
        |z'|   | 0 -sin(-latO) cos(latO) | |z|

        Earth is rotated to line up with observer

                                                  y
                                           y'   ' |
      cos(latO) contribution to y' from y ~ \'90  |*
      (from this pt. to 0,0,0)               \    |*
                                  sin(latO) ~.\   |*
                  (from this pt. to 0,0,0) . 90\  | *
                                          .     \a| . &lt;- -sin(latO) contribution to z' from y
                   Observer            z ________\|.        (from this pt. to 0,0,0)
                                           .     a/  x = x'
                                            .    /
                                             .90/
                                              ./ &lt;- cos(latO) contribution to z' from z
                                              /        (from this pt. to 0,0,0)
                                             z'

        where a / latO is the latitude angle of the Observer
        and 90 is a right angle

        CALCULATING VIEWABLE AREAS AND RADIAL LINE ADJUSTMENT

                                      ....
                                 .            .
                               \                 .
                           .    \                   .
                         .       \                    .        &lt;- earth
                       .          \                     .
                      .            \                     .
                     .              \                     .
                    .                \                     .
                   .                  \                     .
                   .                   \           1        .
                   .                    \ _________________ .
                   .                   / \               90 .' &lt;- horizon (angle between radial
                    .                 /90 \             *  . '    line to horizons and line to center
                     .               /- aO \ ~ 1/altO*   .   '    is 90 degrees)
                      .          1  /       \    *      .    '
                       .           /      90 \*       .      '
                         .        /       *   \     .        '
                           .     / aO  *       \ .           '
                     horizon ~  /. *          . \            '
                                ` 90   ...       \           '
                                  `               \~ altO    '
                                    `              \         '
                                      `             \        ' &lt;- radial line to the horizon
                                        `            \       '    (distance from Observer to horizon
                                          `           \      '    is altO*cos(aO))
                                            `          \     '
                                              `         \    '
                                                `        \   '
                                                   `      \  '
                                                      `  aO\ '
                                                         `  \'
                                                            Observer

     where altO is the altitude of the Observer (referred to as raltO or draltO in equations)
     and 1/altO is the distance from the center to the line connecting the horizon points
     where units are earth radius (as seen by the 1 as the distance from the center to the horizon points)
     and aO is the altitude angle of the Observer (referred to in setup as pang)

     radialLineAdjustment = altO/(altO-zz)
     zvis = 1/altO
     aO = arcsin(zvis)
     rvis = 1/(cos(aO)*cos(aO))  - Note: in setup() this is 1/(cos(pang)*cos(pang))

     A radial line adjustment is needed for the Point's x and y so object is the same
     line-of-sight are at the same pixel

     xx = xx*radialLineAdjustment;
     yy = yy*radialLineAdjustment;

    CHECKING VISIBILITY OF THE POINT
    Is it blocked by the earth from the observers point of view

    to be visible the Point needs to be:
    in front of the observer: zz &gt; altO
    and either in front of the earth: zz &gt; zvis
    or visible along a radial line to the side of the horizon: xx*xx + yy*yy &gt;= rvis
    so: vis = (zz &lt; altO) &amp;&amp; ((zz &gt; zvis) || ((xx*xx + yy*yy) &gt;= rvis))

     </pre>
   */
  @Override
  public int rwc2pix (double[] buf, int off, int dim, int npts, Pix pix) {
    // If incoming data is invalid, print Midas warning
    if ((off + dim*npts) > buf.length) {
      Midas.log("Data is incompatible with View Geo", Midas.WARNING);
      return 0;
    }
    int i,ilon,ilat,j=off,k=0; int[] pax=pix.x, pay=pix.y; float[] paz=pix.z;
    double alon,alat,plon,plat,cp,sp,ct,st,ctcp;
    double mxx=MP.mxx, myy=MP.myy, mxb=MP.mxb, myb=MP.myb;
    double xx,yy,zz,radialLineAdjustment,rr;
    boolean vis, viso=false, chk=(pix.flags&Pix.CHECK)!=0;
    boolean keep = (pix.flags&Pix.KEEP) != 0;
    boolean fullkeep = (pix.flags&Pix.FULLKEEP) != 0; // Part of image in view, so keep this coordinate NXM-3992
    npts = Math.min(npts,pix.max);
    for (i=0; i<npts; i++) {
      alon = buf[j++] - lon + 360.0f; // longitute of point relative to observer (add 360 to make sure degrees are positive)
      alat = buf[j++] + 360.0f; // latitude of point (as positive number)
      ilon = (int)alon;
      ilat = (int)alat;
      plon = alon-ilon; // removed float cast 3.9.1
      plat = alat-ilat; // removed float cast 3.9.1

      // Bug 2769 - Prevent exception when view and data are out of sync
      if (ilon<0 || ilon>=tcos.length || ilon>=dtcos.length || ilon>tsin.length
          || ilon>dtsin.length || ilat<0 || ilat>=tcos.length
          || ilat>=dtcos.length || ilat>tsin.length || ilat>dtsin.length )
        continue;

      cp = tcos[ilon] + plon*dtcos[ilon];
      sp = tsin[ilon] + plon*dtsin[ilon];
      ct = tcos[ilat] + plat*dtcos[ilat];
      st = tsin[ilat] + plat*dtsin[ilat];
      ctcp = ct*cp;  // cos(latitude) * cos(longitude)
      if (dim==2) {
        xx = ct*sp;
        yy = ctcp*slat + st*clat;
        zz = ctcp*clat - st*slat;
        vis = zz > zvis;
        radialLineAdjustment = draltO/(draltO-zz);
        xx = xx*radialLineAdjustment;
        yy = yy*radialLineAdjustment;
      } else {
        // get point's altitude in radius form relative to the earth's radius
        rr = (buf[j++]+equatorialRadius)*inverseEquatorialRad;
        xx = rr*(ct*sp);               // rr * (cos(lat)*sin(lon))
        yy = rr*(ctcp*slat + st*clat); // rr * ((cos(lat)*cos(lon)*sin(-observerlat) + sin(lat)*cos(-observerlat))
        zz = rr*(ctcp*clat - st*slat); // rr * ((cos(lat)*cos(lon)*cos(-observerlat) - sin(lat)*sin(-observerlat))
        radialLineAdjustment = draltO/(draltO-zz);
        xx = xx*radialLineAdjustment;
        yy = yy*radialLineAdjustment;
        vis = zz<draltO && (zz>zvis || xx*xx+yy*yy>=rvis);
      }
      if (vis || fullkeep) {
        pax[k] = (int)(xx*mxx + mxb);
        pay[k] = (int)(yy*myy + myb);
        paz[k] = (float)zz;
        k++;
      }
      else if (chk) { pax[k]=0; pay[k]=0; paz[k]=(float)zz; k++; }
      else if (viso) break;
      viso = vis;
    }
    pix.off=i;
    pix.n = i;
    if (!chk && !keep) k = checkrange(pax,pay,k);
    return k;
  }

  private int checkrange (int[] pax, int[] pay, int kpts) {
    if (kpts==2) {
      if (pax[0]>MP.ix2+100 && pax[1]>MP.ix2+100) kpts=0;
      if (pax[0]<MP.ix1-100 && pax[1]<MP.ix1-100) kpts=0;
      if (pay[0]>MP.iy2+100 && pay[1]>MP.iy2+100) kpts=0;
      if (pay[0]<MP.iy1-100 && pay[1]<MP.iy1-100) kpts=0;
    }
    else if (kpts==1) {
      if (pax[0]>MP.ix2+100) kpts=0;
      if (pax[0]<MP.ix1-100) kpts=0;
      if (pay[0]>MP.iy2+100) kpts=0;
      if (pay[0]<MP.iy1-100) kpts=0;
    }
    return kpts;
  }

  /** A fast version of {@link #pix2rwc(MPoint, double[])}. This method is used by {@link LayerNMS}
      to project imagery onto a GEODETIC view of the Earth. {@link LayerNMS} will call this method
      (on the order of) 100000 times inside a tight loop. For this situation speed is much more
      important than accuracy. Positive radians were used to eliminate two adds and three multiplies
      (between this method and {@link LayerNMS}).
             <b>*** WARNING: lat and lon will be POSITIVE and returned in RADIANS. ***</b>
      @param x   Pixel x-position.
      @param y   Pixel y-position.
      @param rwc Array of length two to receive output lat and lon.
      @return true if the point is off of the Earth,
              false if the point is a valid point on the Earth.
   */
  public boolean fastPix2rwc (int x, int y, float[] rwc) {
    boolean offEarth = false;
    float retz,rrety,rretz;
    float retx = (float)((x-MP.mxb)/MP.mxx);
    float rety = (float)((y-MP.myb)/MP.myy);
    float p2   = raltO*raltO;
    float rad  = retx*retx + rety*rety;
    float fact = (float)quadsolve (p2+rad, -2*p2, p2-1);

    if (fact==0f) {  // off earth
      offEarth = true;
    } else {
      rrety = rety*fact;
      rretz = (float)Math.sqrt(1-rad*fact*fact);
      retx =  retx*fact;
      rety =  (float)(rrety*clat - rretz*slat);
      retz =  (float)(rrety*slat + rretz*clat);

      float rlonRad = DEG2RADF*rlonO;
      retx   = fastAtan2(retx,retz);
      rwc[0] = (rlonRad + retx + PI3) % PI2;
      rwc[1] = PI - (float)Math.asin(rety);
    }
    return offEarth;
  }

  /** A fast approximation of the atan2 function used by fastPix2rwc. This method started
      out as a C method titled __ieee754_atan2 written by SunSoft that is available from
      http://www.netlib.org/fdlibm/e_atan2.c.
      <pre>
        Original Copyright:
          Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.

          Developed at SunSoft, a Sun Microsystems, Inc. business. Permission to use, copy, modify,
          and distribute this software is freely granted, provided that this notice is preserved.

        Original Comments:
          __ieee754_atan2(y,x)
          Method:
            1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
            2. Reduce x to positive by (if x and y are unexceptional):
                ARG (x+iy) = arctan(y/x)           ... if x &gt; 0,
                ARG (x+iy) = pi - arctan[y/(-x)]   ... if x &lt; 0,
          Special cases:
                ATAN2((anything), NaN ) is NaN;
                ATAN2(NAN , (anything) ) is NaN;
                ATAN2(+-0, +(anything but NaN)) is +-0  ;
                ATAN2(+-0, -(anything but NaN)) is +-pi ;
                ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
                ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
                ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
                ATAN2(+-INF,+INF ) is +-pi/4 ;
                ATAN2(+-INF,-INF ) is +-3pi/4;
                ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;

        Changes Made:
          (1) Converted from C to Java. (11 APR 2006)
          (2) Removed most of the special cases since they were not being hit when this
              method was called by fastPix2rwc. (11 APR 2006)
          (3) Changed from using double precision to using float. (11 APR 2006)
      </pre>
   */
  private float fastAtan2(float y, float x) {
    int yy = Float.floatToRawIntBits(y);
    int xx = Float.floatToRawIntBits(x);
    int ix = xx&0x7fffffff;
    int iy = yy&0x7fffffff;
    int m  = ((yy>>31)&1)|((xx>>30)&2); /* 2*sign(x)+sign(y) */
    float z;

    /* compute y/x */
    int k = (iy-ix)>>20;
    if(k > 60) z=PI_OVER_2+0.5f*PI_LO;      /* |y/x| >  2**60 */
    else if(xx<0&&k<-60) z=0.f;             /* |y|/x < -2**60 */
    else z=(float)Math.atan(Math.abs(y/x)); /* safe to do y/x */

    int zz = Float.floatToRawIntBits(z);
    switch (m) { /* can be 0, 1, 2, or 3 */
      case 0:  return z;                                     /* atan(+,+) */
      case 1:  return Float.intBitsToFloat(zz ^ 0x80000000); /* atan(-,+) */
      case 2:  return PI-(z-PI_LO);                          /* atan(+,-) */
      default: return (z-PI_LO)-PI;                          /* atan(-,-) */
    }
  }

  //************************************************************************************************
  //*************** new 3.9.1 high precision double pixel methods **********************************

  // Assumes 1 3D point (overridden in 4.1.3)
  @Override
  int rwc2dpix (double[] buf, DPix pix) {
    return rwc2dpix(buf,0,3,1,pix);
  }

  @Override
  public int rwc2dpix (double[] buf, int off, int dim, int npts, DPix pix) {
    int i,j=off,k=0;
    double[] pax=pix.x, pay=pix.y, paz=pix.z;
    double alon,alat,plon,plat,cp,sp,ct,st,ctcp;
    double mxx=MP.mxx, myy=MP.myy, mxb=MP.mxb, myb=MP.myb;
    double xx,yy,zz,ff,rr;
    boolean vis, viso=false, chk=(pix.flags&Pix.CHECK)!=0;
    boolean keep = (pix.flags&Pix.KEEP) != 0;
    boolean fullkeep = (pix.flags&Pix.FULLKEEP) != 0;
    npts = Math.min(npts,pix.max);
    for (i=0; i<npts; i++) {
      alon = buf[j++] - lon + 360.0f;
      alat = buf[j++] + 360.0f;

      int ilon = (int)alon;
      int ilat = (int)alat;
      plon = alon-ilon; // removed float cast 3.9.1
      plat = alat-ilat; // removed float cast 3.9.1

      // Bug 2769 - Prevent exception when view and data are out of sync
      if (ilon<0 || ilon>=tcos.length || ilon>=dtcos.length || ilon>tsin.length
          || ilon>dtsin.length || ilat<0 || ilat>=tcos.length
          || ilat>=dtcos.length || ilat>tsin.length || ilat>dtsin.length )
        continue;

      cp = tcos[ilon] + plon*dtcos[ilon];
      sp = tsin[ilon] + plon*dtsin[ilon];
      ct = tcos[ilat] + plat*dtcos[ilat];
      st = tsin[ilat] + plat*dtsin[ilat];
      ctcp = ct*cp;  // cos(latitude) * cos(longitude)
      /** cos(latitude) * cos(longitude) */
      ctcp = ct*cp;
      if (dim==2) {
        xx = ct*sp;
        yy = ctcp*slat + st*clat;
        zz = ctcp*clat - st*slat;
        vis = zz > zvis;
        ff = draltO/(draltO-zz);
        xx = xx*ff;
        yy = yy*ff;
      } else {
        rr = (buf[j++]+equatorialRadius)*inverseEquatorialRad;
        xx = rr*(ct*sp);
        yy = rr*(ctcp*slat + st*clat);
        zz = rr*(ctcp*clat - st*slat);
        ff = draltO/(draltO-zz);
        xx = xx*ff;
        yy = yy*ff;
        vis = zz<draltO && (zz>zvis || xx*xx+yy*yy>=rvis);
      }
      if (vis || fullkeep) {
        pax[k] = xx*mxx + mxb;
        pay[k] = yy*myy + myb;
        paz[k] = zz;
        k++;
      }
      else if (chk) { pax[k]=0; pay[k]=0; paz[k]=zz; k++; }
      else if (viso) break;
      viso = vis;
    }
    pix.off=i;
    pix.n = i;
    if (!chk && !keep) k = checkrange(pax,pay,k);
    return k;
  }

  private int checkrange (double[] pax, double[] pay, int kpts) {
    if (kpts==2) {
      if (pax[0]>MP.ix2+100 && pax[1]>MP.ix2+100) kpts=0;
      if (pax[0]<MP.ix1-100 && pax[1]<MP.ix1-100) kpts=0;
      if (pay[0]>MP.iy2+100 && pay[1]>MP.iy2+100) kpts=0;
      if (pay[0]<MP.iy1-100 && pay[1]<MP.iy1-100) kpts=0;
    }
    else if (kpts==1) {
      if (pax[0]>MP.ix2+100) kpts=0;
      if (pax[0]<MP.ix1-100) kpts=0;
      if (pay[0]>MP.iy2+100) kpts=0;
      if (pay[0]<MP.iy1-100) kpts=0;
    }
    return kpts;
  }

  //************************************************************************************************
  //*************** The following pixel methods are not designed for ViewGeo ***********************

  @Override
  int rwc2pix (double[] rwc, MPoint pix, boolean rangeAdjust) {
    Midas.log("rwc2pix with range adjustment designed for latloncontinuous is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2pix(rwc,pix);
  }

  @Override
  int rwc2pix (double[] buf, int off, int dim, int npts, Pix pix, boolean isSpecializedFill) {
    // We do not want the implementation that is used by non geodetic views, use the old way since 3.9.2
    return rwc2pix(buf,off,dim,npts,pix);
  }

  @Override
  int rwc2pixClip (float[] buf, int off, int dim, int npts, Pix pix, boolean mer) {
    Midas.log("rwc2pixClip is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2pix(buf,off,dim,npts,pix);
  }

  @Override
  int rwc2pixClip (double[] buf, int off, int dim, int npts, MPlot MP, Pix pix) {
    Midas.log("rwc2pixClip is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2pix(buf,off,dim,npts,pix);
  }

  @Override
  int rwc2pixClip (double[] buf, int off, int dim, int npts, Pix pix, boolean mer) {
    Midas.log("rwc2pixClip is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2pix(buf,off,dim,npts,pix);
  }

  @Override
  int rwc2pixClipAdjustForFill (double[] buf, int off, int dim, int npts, Pix pix, boolean mer) {
    Midas.log("rwc2pixClipAdjustForFill is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2pix(buf,off,dim,npts,pix);
  }

  int rwc2dpixClip (double[] buf, int off, int dim, int npts, MPlot MP, DPix pix) {
    Midas.log("rwc2dpixClip is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2dpix(buf, off, dim, npts, pix);
  }

  int rwc2dpixClip (double[] buf, int off, int dim, int npts, DPix pix, boolean mer) {
    Midas.log("rwc2dpixClip is NOT designed for ViewGeo", Midas.WARNING);
    return rwc2dpix(buf, off, dim, npts, pix);
  }
}
