//
//   Computes "rho" derivative of a 2D array
//   
//   Uses Centered differences in interior
//   Uses backward/forward differences along boundaries
//
//   Copyright 1998, Steve  Liebling, The University of Texas at Austin
//


// shared root typedef struct {
//    long                *nDim           "Num Dimensions";
//    long*               *dims[nDim]     "Dimensions Array";
//    cxData(nDim, dims)  *data           "Data Structure";
//    cxCoord(cDim, dims) *coord          "Coord Structure";
// } cxLattice;

double drho;
int    nrho;
int    nz;
int    i;
int    j;

nz   = shape(First_In)[0];
nrho = shape(First_In)[1];
drho = (double) ( First_In_coord[0][1] - First_In_coord[0][0] ) / ( nrho - 1 );

write( list( drho, nrho, nz ) );


First_Out := allocate( (double) 1.1, [nz,nrho,1],fill_zero);

for (i = 0; i < nrho; i += 1) {
   for (j = 0; j < nz; j += 1) {
      if ( i != 0 && i != nrho-1 )
      // centered difference
         First_Out[j][i][0] =
            ( First_In[j][i+1][0] - First_In[j][i-1][0] ) / ( 2.0*drho );
      else if ( i == 0 )
      // forwards difference
         First_Out[j][i][0] =
            ( -3.0 * First_In[j][i  ][0]
              +4.0 * First_In[j][i+1][0]
              -      First_In[j][i+2][0] ) / ( 2.0*drho );
      else if ( i == nrho-1 )
      // backwards difference
         First_Out[j][i][0] =
            (  3.0 * First_In[j][i  ][0]
              -4.0 * First_In[j][i-1][0]
              +      First_In[j][i-2][0] ) / ( 2.0*drho );
      //write( list(i,j,First_In[j][i][0],First_Out[j][i][0] ) );
   };
};

First_Out_lat := make_lat(First_In_coord,First_Out);