public static int jauEpv00(final double date1, final double date2,
                      double pvh[][], double pvb[][])
         {
         /*
         * Matrix elements for orienting the analytical model to DE405.
         *
         * The corresponding Euler angles are:
         *
         *                       d  '  "
         *   1st rotation    -  23 26 21.4091 about the x-axis  (obliquity)
         *   2nd rotation    +         0.0475 about the z-axis  (RA offset)
         *
         * These were obtained empirically, by comparisons with DE405 over
         * 1900-2100.
         */
            final double am12 =  0.000000211284,
                                am13 = -0.000000091603,
                                am21 = -0.000000230286,
                                am22 =  0.917482137087,
                                am23 = -0.397776982902,
                                am32 =  0.397776982902,
                                am33 =  0.917482137087;
            
      
        
        
    /* Pointers to coefficient arrays, in x,y,z sets */
       final double ce0[][] = { Ephemeris.e0x, Ephemeris.e0y, Ephemeris.e0z },
                           ce1[][] = { Ephemeris.e1x, Ephemeris.e1y, Ephemeris.e1z },
                           ce2[][] = { Ephemeris.e2x, Ephemeris.e2y, Ephemeris.e2z },
                           cs0[][] = { Ephemeris.s0x, Ephemeris.s0y, Ephemeris.s0z },
                           cs1[][] = { Ephemeris.s1x, Ephemeris.s1y, Ephemeris.s1z },
                           cs2[][] = { Ephemeris.s2x, Ephemeris.s2y, Ephemeris.s2z };
       /* Numbers of terms for each component of the model, in x,y,z sets */
       final int ne0[] = {Ephemeris.e0x.length/3,
               Ephemeris.e0y.length/3,
               Ephemeris.e0z.length/3 },
                        ne1[] = {Ephemeris.e1x.length/3,
               Ephemeris.e1y.length/3,
               Ephemeris.e1z.length/3 },
                        ne2[] = {Ephemeris.e2x.length/3,
               Ephemeris.e2y.length/3,
               Ephemeris.e2z.length/3 },
                        ns0[] = {Ephemeris.s0x.length/3,
               Ephemeris.s0y.length/3,
               Ephemeris.s0z.length/3 },
                        ns1[] = {Ephemeris.s1x.length/3,
               Ephemeris.s1y.length/3,
               Ephemeris.s1z.length/3 },
                        ns2[] = {Ephemeris.s2x.length/3,
               Ephemeris.s2y.length/3,
               Ephemeris.s2z.length/3 };
       int nterms;

    /* Miscellaneous */
       int jstat, i, j;
       double t, t2, xyz, xyzd, a, b, c, ct, p, cp,
              ph[] = new double[3], vh[] = new double[3], pb[] = new double[3], vb[] = new double[3], x, y, z;

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

    /* Time since reference epoch, Julian years. */
       t = ((date1 - DJ00) + date2) / DJY;
       t2 = t*t;

    /* Set status. */
       jstat = abs(t) <= 100.0 ? 0 : 1;

    /* X then Y then Z. */
       for (i = 0; i < 3; i++) {

       /* Initialize position and velocity component. */
          xyz = 0.0;
          xyzd = 0.0;

       /* ------------------------------------------------ */
       /* Obtain component of Sun to Earth ecliptic vector */
       /* ------------------------------------------------ */

       /* Sun to Earth, T^0 terms. */
          nterms = ne0[i];
          int idx;
          for (j = 0, idx=0; j < nterms; j++) {
             a = ce0[i][idx++];
             b = ce0[i][idx++];
             c = ce0[i][idx++];
             p = b + c*t;
             xyz  += a*cos(p);
             xyzd -= a*c*sin(p);
          }

       /* Sun to Earth, T^1 terms. */
          nterms = ne1[i];
          for (j = 0, idx= 0; j < nterms; j++) {
             a = ce1[i][idx++];
             b = ce1[i][idx++];
             c = ce1[i][idx++];
             ct = c*t;
             p = b + ct;
             cp = cos(p);
             xyz  += a*t*cp;
             xyzd += a*( cp - ct*sin(p) );
          }

       /* Sun to Earth, T^2 terms. */
          nterms = ne2[i];
          for (j = 0, idx = 0; j < nterms; j++) {
             a = ce2[i][idx++];
             b = ce2[i][idx++];
             c = ce2[i][idx++];
             ct = c*t;
             p = b + ct;
             cp = cos(p);
             xyz  += a*t2*cp;
             xyzd += a*t*( 2.0*cp - ct*sin(p) );
          }

       /* Heliocentric Earth position and velocity component. */
          ph[i] = xyz;
          vh[i] = xyzd / DJY;

       /* ------------------------------------------------ */
       /* Obtain component of SSB to Earth ecliptic vector */
       /* ------------------------------------------------ */

       /* SSB to Sun, T^0 terms. */
          nterms = ns0[i];
          for (j = 0, idx = 0; j < nterms; j++) {
             a = cs0[i][idx++];
             b = cs0[i][idx++];
             c = cs0[i][idx++];
             p = b + c*t;
             xyz  += a*cos(p);
             xyzd -= a*c*sin(p);
          }

       /* SSB to Sun, T^1 terms. */
          nterms = ns1[i];
          for (j = 0, idx = 0; j < nterms; j++) {
             a = cs1[i][idx++];
             b = cs1[i][idx++];
             c = cs1[i][idx++];
             ct = c*t;
             p = b + ct;
             cp = cos(p);
             xyz  += a*t*cp;
             xyzd += a*(cp - ct*sin(p));
          }

       /* SSB to Sun, T^2 terms. */
          nterms = ns2[i];
          for (j = 0, idx = 0; j < nterms; j++) {
             a = cs2[i][idx++];
             b = cs2[i][idx++];
             c = cs2[i][idx++];
             ct = c*t;
             p = b + ct;
             cp = cos(p);
             xyz  += a*t2*cp;
             xyzd += a*t*(2.0*cp - ct*sin(p));
         }

       /* Barycentric Earth position and velocity component. */
         pb[i] = xyz;
         vb[i] = xyzd / DJY;

       /* Next Cartesian component. */
       }

    /* Rotate from ecliptic to BCRS coordinates. */

       x = ph[0];
       y = ph[1];
       z = ph[2];
       pvh[0][0] =      x + am12*y + am13*z;
       pvh[0][1] = am21*x + am22*y + am23*z;
       pvh[0][2] =          am32*y + am33*z;

       x = vh[0];
       y = vh[1];
       z = vh[2];
       pvh[1][0] =      x + am12*y + am13*z;
       pvh[1][1] = am21*x + am22*y + am23*z;
       pvh[1][2] =          am32*y + am33*z;

       x = pb[0];
       y = pb[1];
       z = pb[2];
       pvb[0][0] =      x + am12*y + am13*z;
       pvb[0][1] = am21*x + am22*y + am23*z;
       pvb[0][2] =          am32*y + am33*z;

       x = vb[0];
       y = vb[1];
       z = vb[2];
       pvb[1][0] =      x + am12*y + am13*z;
       pvb[1][1] = am21*x + am22*y + am23*z;
       pvb[1][2] =          am32*y + am33*z;

    /* Return the status. */
       return jstat;

        }