/*=========================================================================

  Program:   Visualization Toolkit
  Module:    grid3D.cxx
  Language:  C++
  Date:      Date: 08/15/2002 17:00:00 
  Version:   Revision: 1.0 

Copyright (c) 1993-2001 Ken Martin, Will Schroeder, Bill Lorensen 
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
   of any contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

 * Modified source versions must be plainly marked as such, and must not be
   misrepresented as being the original software.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

=========================================================================*/

//vtk include files
#include "vtkPolyDataSource.h"
#include "vtkPoints.h"
#include "vtkDataArray.h"
#include "vtkFloatArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"

//custom header files
#include "grid3D.h"

//C++ header files
#include <math.h>

//    const long double DegToRad = vtkMath::Pi() / 180.0;

/************************************************************************

This is basically what the grid will look like mathematically
minus the distortions.  


                Small Grid Representation

            4   .     .     .      .      .      .   

    
            3   .     .     .      .      .      .


(j)         2   .     .     .      .      .      .


            1   ._____.     .      .      .      .
                |     |
                |     |
            0   |_____|     .      .      .      .
    
                0     1     2      3      4      5

                                                    (i)

**In this small representation each Pt represents a cell
                
        EXCEPT the southernmost row j=0.  

Therefore Cell counting should go    (i) 0-179 and (j) 1-136

Pt. (1,1) is the cell shown which would be cell ID #1

The 0th cell would be represented by Pt. (0,1) and be the cell to the left of the starting line.  

I think it would make more sense to label each cell by the NW corner.  But we will see.  


iNumPts = (from NcFile) longitudeDimension->size() = (for now) 180 = iNumCells

jNumPts = (from NcFile) latitudeDimension->size() = (for now) 137
  && jNumCells = (for now) 136 :: one less


************************************************************************/


//----------------------------------------------------------------------------

// Construct sphere with radius = 1.01 and default resolution 8 in both Phi
// and Theta directions. Theta ranges from (0,360) and phi (0,180) degrees.
grid3D::grid3D(long res) : DegToRad(vtkMath::Pi() / 180.0)
                         // initialize list for const
{
    res = res < 4 ? 4 : res;    // check to make sure theat and phi resolutions are not too small

    this->Radius = 1.0;          // overridden in eCubed

    this->Center [0] = 0.0;     //      |
    this->Center [1] = 0.0;     //       }  NO IDEA what these do...they stay 0 the entire time
    this->Center [2] = 0.0;         //      |

    this->ThetaResolution = res;    
    this->StartTheta = 0.0;
    this->EndTheta = 360.0;

    this->PhiResolution = res;
    this->StartPhi = 0.0;
    this->EndPhi = 180.0;


    this->OnRatio = 1; //used to increment in loops
}


//----------------------------------------------------------------------------

// Retrieve grid3D elevation values from netCDF file
// this is called manually before Execute and Execute Information
void grid3D::ReadElevationFromFile(char* inp, char* dat, char* lon, char* lat)
{

    cout << "NcFile inp = " << inp << endl;
    cout << "NcVar elevation = " << dat 
         << ", NcVar longitude = " << lon
         << ", NcVar latitude = " << lat << endl;

    // Find existing file
    this->elevationFile = new NcFile(inp, NcFile::ReadOnly);
    if (this->elevationFile->is_valid())
    {
      cerr << "\nReading POP grid3D data from netCDF file " << inp << endl;
    }

    //......................................................................
    // Check everything before retrieving any data

    // stores what will be iNumCells which = number of lines of longitude
    NcDim* longitudeDimension = this->elevationFile->get_dim("longitude");
    cerr << longitudeDimension->name() << " size is "
         << longitudeDimension->size() << " points." << endl;

    // stores what will be jNumCells which is 1 less than # of lines of latitude
    NcDim* latitudeDimension = this->elevationFile->get_dim("latitude");
    cerr << latitudeDimension->name() << " size is "
         << latitudeDimension->size() << " points." << endl;

    NcVar* longitude = this->elevationFile->get_var(lon);

    NcVar* latitude = this->elevationFile->get_var(lat);

    NcVar* elevation = this->elevationFile->get_var(dat);

    //NcAtt* unitsOfElevation = elevation->get_att("units");

    //NcAtt* longNameOfElevation = elevation->get_att("long_name");
    //cerr << longNameOfElevation->as_string(0) << " is in "
    //     << unitsOfElevation->as_string(0) << endl;

    //.................... end of checks ...................................


    /***********************************************************************
     These lengths determine the number of unique points in grid3D.  
     In the E-W direction (i), this is <<equal>> to the # of cells
     b/c of wrap-around.
     In the N-S direction (j), this is <<one more>> than the # of cells. 
     **********************************************************************/

     //.......... convert NetCDF data to local variables ...................

    this->iNumPts = longitudeDimension->size(); // number of lines of longitude
    this->iNumCells = this->iNumPts;            // number of horizontal cells

    this->jNumPts = latitudeDimension->size();  // number of lines of latitude
    this->jNumCells = -1 + this->jNumPts;       // number of vertical cells

    cerr << "NetCDF file " << inp << " has "
         << iNumCells << " E-W cells by "
         << jNumCells << " N-S cells." << endl;


    this->totalNumPts = this->iNumPts*this->jNumPts;
    this->totalNumCells = this->iNumCells*this->jNumCells;


    /**********************************************************************
     Allocate arrays for data values that will be used to construct
     the global grid3D.  These will refer to the 'NE' corners of the
     cells except for the southernmost latitude.  There a SE corner
     location will be calculated and stored for each such south point. 

     There is now a BIG difference between point data and cell data.
     Point data is what is retreived from the NetCDF file.  Cell data
     is what the editor will use to display the grid.  Before, I think
     the association among the two was incorrect.  The cell data was
     being mapped from the bottom right (SE) point of the cell instead
     of the NE point.  This has been fixed now as far as I can
     tell by making very clear the difference among point and cell data.  
        
     **********************************************************************/

     //............... Retrieve elevation data .............................

    // Variables
    long iCell, jCell, iPt, jPt, dummy;


        // ................ Point Arrays ...................
            
            //2-D arrays
                // I think all of the point arrays are 1-D.  

            //1-D arrays

                /*************** work array ********************
                the work array is used as a scratch array into which
                elevation, latitude and longitude from the NetCDF file
                are temporarily set like so 

                ex.
                    elevation->get(&workArray [0] , elevation->edges());

                This allows the values to be transerred to a Point Array
                such as pointElevation whch then can be in turn put into
                the cellElevation array which is needed by the editor.  

                ************************************************/
                float* workArray = new float [this->totalNumPts] ;

                // these arrays help transfer lat and long to cell arrays
                double* thetaValue = new double [this->totalNumPts] ;
                double* phiValue = new double [this->totalNumPts] ;
   
                // These give the lat and long of NE point of each cell
                pointTheta = new float [this->totalNumPts];
                pointPhi = new float [this->totalNumPts];

                // this holds the elevation of each point--not each cell,
                // a big confusion earlier.  This must be correctly
                // mapped to cellElevation
                this->pointElevation = new float [this->totalNumPts] ;


        //................. Cell Arrays......................

            // 2-D arrays

                // hold the unique cell ID for each cell
                // will be made 2-D later on in this function
                this->cellId = new long* [this->iNumCells] ;

                //I guess each point has an Id too but I need to look
                // at this more closely
                // will be made 2-D later on in this function
            //  this->cellPointId = new long* [this->totalNumCells] ;

            // 1-D array
                this->cellIndex = new IdToIndex [this->totalNumCells] ;

                // holds the lat and long of each cell
                cellTheta = new float [this->totalNumCells];
                cellPhi = new float [this->totalNumCells];

                //holds the elevation of a cell based on NE point of cell
                this->cellElevation = new float [this->totalNumCells] ;
    
        
    //............... Generate 'T' cell mesh connectivity ..................

    /***********************************************************************
     This upcoming block of code creates the 2-D array of cell id numbers.

     Given array indices [iPt,jPt] , we will be able to retrieve a unique
     cell id.  If instead given the cell id, we will be able to calculate
     the indices [iCell,jCell].
     This block also associates a subarray of four cell points with each cell.
     The assignment order of the index "dummy" is arbitrary in this loop.  
     However, for clarity, I have preserved the order of the actual data: 
            
              [iPt=0, jPt=1] = [iCell = 0, jCell = 0] = cellId is 0

     iCell, iPt, jCell and jPt all count from zero (0).  However, jCell
     is one less than jPt for any cell [iCell] [jCell] whose NE corner
     is [iPt] [jPt] on the first row.

     On subsequent rows jCell is one less than a factor of INumCells*jPt.  


        **for commenting:  I have changed this around some and think I
        found the crux of the numbering problem.  Everytime the j pt is
        referring to cells it must start at 1 instead of 0 because there
        are cells defined by j=0.  However this is currently how it is
        being mapped so that causes problems for ferret and prob Pop as well.  
        
        
        **this loop sets up the cell ID array counting from bottom to top,
        left to right

       ********************************************************************/
        
        
    for (iCell=0; iCell<iNumCells; iCell++)
    {
        // makes the cellId array 2-D by tacking on j value to each i value
        // after loop is over cellId will be 'cellId [i] [j]'
        this->cellId [iCell] = new long [jNumCells] ;

        for (jCell=0; jCell<jNumCells; jCell++)
        {
            //algorithm that iterates through cells one by one l-r, b-t and 
            // assigns each cell the next counting number, the cellId
            dummy = iCell + jCell*iNumCells;

            // associates the cellId with an [iCell][jCell] index
            this->cellId [iCell]  [jCell] = dummy;
        } //end of for(j=0; j<jNumCells; j++)
    }  //end of for (i=0; i<iNumCells; i++)



    //.................. POINT Elevation of cell in meters .................

    //this sets elevation values from netCDF file into the work array
    //see above for description
    elevation->get(&workArray [0] , elevation->edges());

    // Loops iterate over all cells associating elevation with an Id
    for (iPt=0; iPt<iNumPts; iPt++)
    {
        for (jPt=0; jPt<jNumPts; jPt++)
        {
            dummy = iPt + jPt*iNumPts;
            this->pointElevation [dummy] = workArray [dummy] ;
        }
    }


    //.................... POINT Longitude in degrees .......................

    //this sets longitude(theta) values from netCDF file into the
    //work array see above for description
    longitude->get(&workArray [0] ,longitude->edges());

    // Loops iterate over all cells to assign POINT elevations
    for (jPt=0; jPt<jNumPts; jPt++)
    {
        for (iPt=0; iPt<iNumPts; iPt++)
        {
            dummy = iPt+jPt*iNumPts;
            this->pointTheta [dummy] = workArray [dummy];
        }
    }
    // Now you can retrieve the longitude of any point (i=0:iNumPts,
    // j=0:jNumPts).  Note that the number of j-points is one greater
    // than the number of j-cells.


    //.................... POINT Latitude in degrees .......................

    // Loops iterate over all points to assign each one a latitude
    latitude->get(&workArray [0] ,latitude->edges());

    for (jPt=0; jPt<jNumPts; jPt++)
    {
        for (iPt=0; iPt<iNumPts; iPt++)
        {
            dummy = iPt+jPt*iNumPts;
            this->pointPhi [dummy] = workArray [dummy];

        }
    }


    /************************ POINT to CELL conversions ********************

    I think and hope that this solves the problem.  This last loop
    relates the cell elevation to the point elevation of the point
    above it and to the right (ie NE point of cell).  

    I am not sure if we need a plus 1 in the [i + this->iNumCells].
    We might.

    Lets see, we do not.  The 0-th (i) pt and 1-st (j) pt are cell # 0.
    The top right corner of cell 0 is pt i+iNumCells

    One other problem is the difference between the counter variable i and j and
    the index points i and j.  They DO NOT ALWAYS MATCH UP.  
            
    Therefore the northern boundary problem should now work in ferret. 

    ***************************************************************************/
        
    // associates the cell with its NE point!!
    for(iCell=0; iCell<this->totalNumCells; iCell++)
    {       
        this->cellElevation [iCell] =
            this->pointElevation [iCell + this->iNumPts] ;
        this->cellTheta [iCell] =
            this->pointTheta [iCell + this->iNumPts] ;
        this->cellPhi [iCell] =
             this->pointPhi [iCell + this->iNumPts] ;
    }

    /************************ find min, max CELL elevations ********************
    This information is needed for setting color mapping table
    ***************************************************************************/

    maxElevation = -99999.9;
    minElevation = 99999.9;
    for(iCell=0; iCell<this->totalNumCells; iCell++)
    {       
        dummy = this->cellElevation [iCell];
        if (dummy > maxElevation) maxElevation = dummy;
        if (dummy < minElevation) minElevation = dummy;
    }
    cout << "Maximum data value " << this->maxElevation
         << ", minimum " << this->minElevation << endl;

    this->elevationFile->~NcFile();
    return;
}


//----------------------------------------------------------------------------
void grid3D::ExecuteInformation()
{
    long numTris, numPts;
    unsigned long size;

    // ignore poles
    numPts = this->ThetaResolution * (this->PhiResolution + 1);
    numTris = this->ThetaResolution * this->PhiResolution * 2;
    size = numPts * 3 * sizeof(float);
    size += numTris * 4 * sizeof(long);

    // convert to kilobytes
    size = (size / 1000) + 1;

}



//----------------------------------------------------------------------------

// Execute is called automatically as soon as an object is created
void grid3D::Execute()
{
    // counter variables
    long i, j;

    //this->totalNumPts is TotalNumPts, this->totalNumPts is TotalNumCells
    // this->totalNumPts = this->iNumPts*jNumPts;
    //this->totalNumCells = this->iNumCells*jNumCells;

    //vtk objects to display points and cells
    vtkPoints *newPoints;
    vtkCellArray *newPolys;

    float  x [3] , n [3] ,  norm, radius, phi, theta;

    long numPoles=0;

    vtkIdType pts [4] ;

    vtkPolyData *output = this->GetOutput();

    long piece = output->GetUpdatePiece();

    long numPieces = output->GetUpdateNumberOfPieces();

//..........................................................................
// All are now derived from netCDF.
    long iphys, jphys, iphys_beg, iphys_end, jphys_beg, jphys_end;

    iphys = this->iNumPts;
    iphys_beg = 0;
    iphys_end = this->iNumPts-1;

    jphys = this->jNumPts;
    jphys_beg = 0;
    jphys_end = this->jNumPts-1;

    cout << "C Convention 0-based indexing\n";
    cout << "Number of I, J " << this->iNumPts << ' ' << this->jNumPts << endl;
    cout << "I Range includes " << iphys_beg << " to " << iphys_end << endl;
    cout << "J Range includes " << jphys_beg << " to " << jphys_end << endl;

    // Allocate 1-D arrays for data values
    double* thetaValue = new double [this->totalNumPts] ;
    double* phiValue = new double [this->totalNumPts] ;
    double* thetaMetric = new double [this->totalNumPts] ;
    double* phiMetric = new double [this->totalNumPts] ;

    // The total number of values is jphys*iphys, but the number of
    // visualization points and cells is smaller.

    //long* indexLon = new long [this->totalNumCells] ; // Used only here
    //long* indexLat = new long [this->totalNumCells] ; // Used only here
    // Begin at southernmost latitude
    for (j=0; j<this->jNumPts; j++)  // Step north thru 'latitude'
    {
        int dummy = 0;
        i=0;  // Begin at westernmost point
        for (i=0; i<this->iNumPts; i++)  // Step east thru 'longitude'
        {
            dummy = i+j*this->iNumPts;

            theta = this->pointTheta [dummy];
            thetaValue [dummy] = theta;

            phi = this->pointPhi [dummy];
            phiValue [dummy] = phi;
        }
    }

//---------------------------------------------------------------------

    // Set things up; allocate memory
    //
    vtkDebugMacro("SphereReader Executing");

    newPoints = vtkPoints::New();
    newPoints->Allocate(this->totalNumPts);


    long maxPtsPerCell = 4;
    newPolys = vtkCellArray::New();
    newPolys->Allocate(newPolys->EstimateSize(this->totalNumCells,
                       maxPtsPerCell));

    long totalPoints = 0;
    long totalCells = 0;

    // Create sphere

    // First create auxiliary arrays to hold ID numbers and
    // other data to make it easier to work with the grid.

    this->pointId = new long* [this->iNumPts] ;
    if (!this->pointId) exit(1);

    pointIndex = new IdToIndex [this->totalNumPts] ;
    if (!pointIndex) exit(1);

    pointElevation = new float [this->totalNumPts] ;
    if (!pointElevation) exit(1);

    for (i=0; i<this->iNumPts; i++)
    {
        this->pointId [i] = new long [this->jNumPts] ;
        if (!this->pointId [i] ) exit(1);
        for (j=0; j<jNumPts; j++)
        {
            this->pointId [i]  [j] = -99;
            pointElevation [i*j] = -9999.0;
            pointTheta [i*j] = -9999.0;
            pointPhi [i*j] = -9999.0;
        }
    }


    // Intermediate points
    for (i=0; i<this->iNumPts; i++)  // Go from west to east
    {
        for (j=0; j<jNumPts; j++, totalPoints++)  // Go from south to north
        {
            long ii = j*this->iNumPts + i;
            //phi = pointPhi [totalPoints];
            //theta = pointTheta [totalPoints];
            phi   = phiValue [ii] ;
            theta = thetaValue [ii] ;

            radius = fabs(this->Radius * cos((double)(DegToRad*phi)));
            n [0] = radius * cos((double)(DegToRad*theta));
            n [1] = radius * sin((double)(DegToRad*theta));
            n [2] = this->Radius * sin((double)(DegToRad*phi));
            x [0] = n [0] + this->Center [0] ;
            x [1] = n [1] + this->Center [1] ;
            x [2] = n [2] + this->Center [2] ;
            newPoints->InsertNextPoint(x);

            this->pointId [i]  [j] = totalPoints;
            pointPhi [totalPoints] = phi;
            pointTheta [totalPoints] = theta;
            float elevation = 0.0;
            pointElevation [totalPoints] = elevation;

            if ( (norm = vtkMath::Norm(n)) == 0.0 )
            {
                norm = 1.0;
            }
            n [0] /= norm; n [1] /= norm; n [2] /= norm;
        }
    } 
        


//********** Generate 'T' cell mesh connectivity between the poles *********

  
   long x1, x2, y1, y2;
   for (j=0; j<this->jNumCells; j+=OnRatio)
   {
      for (i=0; i<this->iNumCells; i+=OnRatio)
      {
          if (j >= this->jNumCells)// Along the northern seam ** took out -1
          {
              if(i > (this->iNumCells/2)-1)  // Return side
              {
                  // Off-seam points for current cell
                  x1 = i;
                  x2 = (i+OnRatio) % this->iNumCells;
                  y1 = j;

                  pts [0] = this->pointId [ x1 ]  [ y1  ] ;
                  pts [3] = this->pointId [ x2 ]  [ y1 ] ;

                  // On-seam points from opposing cell
                  x1 = this->iNumCells-i;
                  x2 = (this->iNumCells-(i+OnRatio)) % this->iNumCells;
                  y2 = j+OnRatio;
                  if (y2 > this->jNumCells) y2 = this->jNumCells;

                  pts [1] = this->pointId [ x1 ]  [ y2  ] ;
                  pts [2] = this->pointId [ x2 ]  [ y2 ] ;

              } //end of if(i > (this->iNumCells/2)-1)

              else  // Starting side
              {
                  x1 = i;
                  x2 = (i+OnRatio) % this->iNumCells;
                  y1 = j;
                  y2 = j+OnRatio;
                  if (y2 > this->jNumCells) y2 = this->jNumCells;

                  pts [0] = this->pointId [ x1 ]  [ y1  ] ;
                  pts [1] = this->pointId [ x1 ]  [ y2  ] ;
                  pts [2] = this->pointId [ x2 ]  [ y2 ] ;
                  pts [3] = this->pointId [ x2 ]  [ y1 ] ;

              } //end of else

          } //end of if (j >= this->jNumCells)

          else // anywhere else on the globe
          {
              x1 = i;
              x2 = (i+OnRatio) % this->iNumCells;
              y1 = j;
              y2 = j+OnRatio;
              if (y2 > this->jNumCells) y2 = this->jNumCells;

              pts [0] = this->pointId [ x1 ]  [ y1  ] ;
              pts [1] = this->pointId [ x1 ]  [ y2  ] ;
              pts [2] = this->pointId [ x2 ]  [ y2 ] ;
              pts [3] = this->pointId [ x2 ]  [ y1 ] ;

          } //end of else

          newPolys->InsertNextCell(4, pts);
          cellId [i]  [j] = totalCells;
          cellIndex [cellId [i]  [j] ] .I = i;
          cellIndex [cellId [i]  [j] ] .J = j;

          /*********************** PROBLEM ******************************

          either this is old code or something but it should not be 
          commented out.  This is the part that associates the four
          corner points of each cell with that cell and right now it 
          is not behaving.  I have commented it out and apparently it 
          doesn't affect the editor too much.  
                            
          I don't know if this info is vital for POP or the NetCDF file
          but now it is not working.  

          **************************************************************/

         /*   for (k=0; k<4; k++)
            {
                long dummyId = cellId [i]  [j] ;
               // cellPointId [dummyId]  [k] = pts [k] ;
            } */


          totalCells++;
       } //end of for (i=0; i<this->iNumCells; i+=OnRatio)

    } //end of for (j=0; j<this->jNumCells; j+=OnRatio)

    cout << "Reader Total Points, Cells: "
         << totalPoints << ", "
         << totalCells << endl;

    // Update ourselves and release memory
    //

    //this is the same as resetScalarMapper
    //maybe a function call here would be better instead of code twice.  
    vtkFloatArray *newElevations;
    newElevations = vtkFloatArray::New();
    newElevations->SetNumberOfComponents(1);
    newElevations->SetNumberOfTuples(totalCells);
    newElevations->SetName("Elevations");
    for (i = 0; i < totalCells; i++)
    {
        newElevations->SetTuple1(i,cellElevation[i]);
    }

    newPoints->Squeeze();
    output->SetPoints(newPoints);
    output->GetCellData()->SetScalars(newElevations);
    newPoints->Delete();

    output->SetPolys(newPolys);
    newPolys->Delete();
}






//----------------------------------------------------------------------------
// Calculate angle A in degrees between vector b and true east c on earth
// given points in longitude degrees East and latitudes in degrees N/S (+/-).
// First (lonA, latA) pair is the common origin of the two vectors.
// Second pair (lonC, latC) is the end point of vector b.  The end point
// (lonB, latB) of vector c is taken as (lonC, latA) to form a right spherical
// triangle where B is the 90-degree angle.
//
/*
                            C
                           /|
                         /  |
                   b   /    |
                     /      | a
                   /        |
                 /          |
               /            |
             A-------------->B
                     c
*/

long double
grid3D::SphericalAngle(double lon0, double lat0,
                               double lon1, double lat1)
{
    long double aValue, lonA, lonB, latB, latC; // lonC, latA,

    if (lon0 > lon1)
    {
        lonA = DegToRad * (lon0 - (long double)360.0);
    }
    else
    {
        lonA = DegToRad * lon0;
    }

    //lonC = DegToRad * lon1;
    lonB = DegToRad * lon1;

    //latA = DegToRad * lat0;
    latB = DegToRad * lat0;
    latC = DegToRad * lat1;

    // No matter what the grid3D, c lies along the east reference direction,
    // a line of true latitude.  We want the angle between it and the
    // arbitrary b vector.  Even though b lies along a grid line in
    // the transformed grid3D system, inside this routine it is merely
    // an ordinary vector whose true coordinates are known.  Since we
    // choose b's end point to lie on the same true meridan as c's end
    // point, the a vector lies along that true meridian.
    long double c = (lonB - lonA) * cos(latB);
    long double a = latC - latB;
    long double b = sqrt(a*a + c*c);  // Approximate

    // However, we don't have a right spherical triangle because c does
    // not lie along a great circle.  So the following formula is slightly
    // off.
    long double sinA = sin(a) / sin(b); // sinB = 1
    aValue = asin(sinA);

    return aValue;
}

//----------------------------------------------------------------------------
// Calculate distance in meters between two points on idealized spherical earth
// given longitudes in degrees East and latitudes in degrees N/S (+/-).
double
grid3D::SphericalDistance(double lon1, double lat1, double lon2, double lat2)
{
    double aValue;

    long double lonA = DegToRad * lon1;
    long double latA = DegToRad * lat1;
    long double lonB = DegToRad * lon2;
    long double latB = DegToRad * lat2;

    long double dLon = fabs( (long double)(lonA - lonB) );
    long double cosDLon = cos(dLon);
    long double cosLatA = cos(latA);
    long double cosLatB = cos(latB);
    long double sinLatA = sin(latA);
    long double sinLatB = sin(latB);

    aValue =
        (
        (cosDLon * cosLatA * cosLatB)
        +
        (sinLatA * sinLatB)
        );

    aValue = acos(aValue);

    aValue *= earthRadius;

    return aValue;
}






//----------------------------------------------------------------------------

/**********************************************************************

this is a part gleaned from grid::Execute that maps the color topography
by a scalar value, the ocean depth.  At first this only occurred once,
the first time the object was created.  This was fine until we wanted
immediate updates.  The only time then that the program would update is
when you would exit and then come back in.  Once this small but vital
chunk of code was isolated here, any topography update method could call
this function, for it was made public, and update the topography.  

This function takes all of the new cell elevations, after some point and
click or elevation change and remapps the grid according to a preset
color scheme in topography.cxx.  this part merely says it wants to set
the Scalars as the newElevations.  Later the scalars are color-mapped.  

*********************************************************************/
void grid3D::ResetScalarMapper()
{
    vtkPolyData *output = this->GetOutput();

    //computes all of the cells
    long totalCells = (this->jNumCells)*(this->iNumCells);

    vtkFloatArray *newElevations;
    newElevations = vtkFloatArray::New();
    newElevations->SetNumberOfComponents(1);
    newElevations->SetNumberOfTuples(totalCells);
    newElevations->SetName("Elevations");

    // this loop puts all of the cellElevations into newElevations
    // which is then set to scalars
    for (long i = 0; i < totalCells; i++)
    {
        newElevations->SetTuple1(i,cellElevation[i]);
    }

    // the line that makes newElevations the scalar value to be mapped.  
    output->GetCellData()->SetScalars(newElevations);
}



//----------------------------------------------------------------------------
// Store grid3D elevation values in current netCDF file
void
grid3D::WriteElevationToFile(char* inp, char* dat)
{

    cout << "NcFile inp = " << inp << endl;
    cout << "NcVar elevation = " << dat << endl;

    /***********************************************************************

     WE STILL HAD THE SAME PROBLEM WE HAD IN READING BEFORE---FIXED

    It has now been fixed. The problem was that we needed to associate each
    cell with its NE corner.  This was done by setting the first iNumCells
    from pointElevation back into the workArray and then putting the rest
    of the cellElevation into the workArray.  This means that the bottom
    row of points which cannot be changed stay put while all the other
    points reflect the value of the cell to their southwest.  

    With this setup I initially thought there would be a problem at the
    northern boundary.  Common sense told me that any point on the
    northern boundary would determine two cells one on each side of the
    boundary.  This however does not seem to be the case and is quite baffling.

    I thought that if I changed the elevation of one cell on the boundary,
    it would then go and set its NE corner point to that elevation (cause
    the NetCDF only stores point data).  This would in turn cause a cell
    on the opposite end of the boundary to also when being mapped, change.
    This however was not the case at all.  hmmmmm

    ************************************************************************/
    float* workArray;

    // longitudeValue, latitudeValue, tempValue;
    long i, dummy;
    workArray = new float [(iNumPts)*(jNumPts)] ;
    dummy = 0;
    this->elevationFile = new NcFile(inp, NcFile::Write);
    //NcDim* longitudeDimension = this->elevationFile->get_dim("longitude");
    //NcDim* latitudeDimension = this->elevationFile->get_dim("latitude");
    NcVar* elevation = this->elevationFile->get_var(dat);

    cout << "Current netCDF file " << inp
        << " has " << iNumPts << ", " << jNumPts << endl;

   // Elevation of cell in meters
   // Iterate over all cells.

     /***********************************************************************
       These last loops put all of the altered cell elevations back to the
       correct pts so that they can be properly stored in the NetCDF file.
       The first iNumCells points cannot be altered as they are associated
       with no cell; therefore they merely remain.  
     ************************************************************************/

    //writes unchangeable points into workArray
    for(i=0; i<iNumCells; i++)
    {
        workArray[i]= this->pointElevation[i];
    }

    // writes all changed (or unchanged) cellElevations back into workArray
    for(i=0; i<(this->iNumCells*this->jNumCells); i++)
    {       
        workArray [i + this->iNumPts] = this->cellElevation [i] ;
    }

    // writes the workArray, now pointElevations, back into NetCDF format.
    // 'put' ends define mode.
    elevation->put(&workArray [0] ,elevation->edges());

    this->elevationFile->~NcFile();
    return;
}

