#!/usr/bin/env python

import os
import re
import math

# ----------------------------------------------------------------------------------------------------

import BaseObjectSimTk
import UtilitiesSimTk
import LogSimTk 

endOfLine                  = '\n'
machineEpsilon             = 1.0e-05

# ----------------------------------------------------------------------------------------------------

#   VGrid_readDX code from file apbs/src/mg/vgrid.c python script in tools/python/vgrid
#   VGrid_writeDX also there

#   for (i=0; i<thee->nx; i++) {
#      for (j=0; j<thee->ny; j++) {
#         for (k=0; k<thee->nz; k++) {
#            u = k*(thee->nx)*(thee->ny)+j*(thee->nx)+i;
#            VJMPERR2(1 == Vio_scanf(sock, "%s", tok));
#            VJMPERR1(1 == sscanf(tok, "%lf", &dtmp));
#            (thee->data)[u] = dtmp; 
#         }       
#      }       
#   }
#
#   0       -> [0, 0, 0]
#   1       -> [1, 0, 0]
#   ...
#   nx-1    -> [nx-1, 0, 0]
#   nx      -> [0, 1, 0]
#   nx+1    -> [1, 1, 0]
#   ...
#   nx*ny   -> [0, 0, 1]
#

# ----------------------------------------------------------------------------------------------------

class GridSimTk(BaseObjectSimTk.BaseObjectSimTk): # {
   """ potential"""

   # -------------------------------------------------------------------------------------------------

   def __init__( self ): # {
      """Constructor for GridSimTk"""
     
      BaseObjectSimTk.BaseObjectSimTk.__init__( self )

      self._xOrigin                      = None
      self._yOrigin                      = None
      self._zOrigin                      = None

      self._xDelta                       = None
      self._yDelta                       = None
      self._zDelta                       = None

      self._xGrid                        = None
      self._yGrid                        = None
      self._zGrid                        = None

      self._xOffset                      = None
      self._yOffset                      = None
      self._zOffset                      = None

      self._totalGridPoints              = None
      self._extrema                      = None
      self._boxRange                     = None
      self._chargeDensity                = None
      self._numberDensity                = None
      self._workingDirectory             = None
      self._gridPointList                = None
      self._potentialList                = None
      self._commentList                  = None

      self._gridChargeDensityFileName    = None
      self._gridNumberDensityFileName    = None
      self._gridPotentialDxFileName      = None
      self._gridPotentialFileName        = None

   # } end of __init__(

   # -------------------------------------------------------------------------------------------------

   # set working directory 

   def setWorkingDirectory( self, inputWorkingDirectory ): # {
      """Set accessor for working directory"""
      backslash = re.compile( '/$' )
      if backslash.match( inputWorkingDirectory )== None:
         inputWorkingDirectory = inputWorkingDirectory + '/'
      self._workingDirectory = inputWorkingDirectory

   # } end of setWorkingDirectory

   def getWorkingDirectory( self ): # {
      """Get accessor for working directory"""
      return self._workingDirectory

   # } end of getWorkingDirectory

   # -------------------------------------------------------------------------------------------------

   # helper method for setting full path name

   def _getFullWorkingPathName( self, inputFileName ): # {

      """Helper method for setting full path name; only prepends working directory name
         if working directory is set and inputFileName does not begin w/ '/'
      """

      frontSlashMatch  = re.compile( '^/' )
      workingDirectory = self.getWorkingDirectory()
      if workingDirectory != None and  frontSlashMatch.match( inputFileName ) == None:
         inputFileName = workingDirectory + inputFileName
      return inputFileName

   # } end of _getFullWorkingPathName

   # -------------------------------------------------------------------------------------------------

   # file name accessors for class -- grid potential

   def setGridPotentialFileName( self, inputFileName, includeWorkingDirectory ): # {
      """Set accessor for grid potential file name"""
      if includeWorkingDirectory:
         inputFileName = self._getFullWorkingPathName( inputFileName )
      self._gridPotentialFileName = inputFileName

   # } end of setGridPotentialFileName

   def getGridPotentialFileName( self ): # {
      """Get accessor for grid potential file name"""
      return self._gridPotentialFileName

   # } end of getGridPotentialFileName

   # -------------------------------------------------------------------------------------------------

   # file name accessors for class -- grid potential dx-format

   def setGridPotentialDxFileName( self, inputFileName, includeWorkingDirectory ): # {
      """Set accessor for grid potential file name"""
      if includeWorkingDirectory:
         inputFileName = self._getFullWorkingPathName( inputFileName )
      self._gridPotentialDxFileName = inputFileName

   # } end of setGridPotentialDxFileName

   def getGridPotentialDxFileName( self ): # {
      """Get accessor for grid potential file name"""
      return self._gridPotentialDxFileName

   # } end of getGridPotentialDxFileName

   # -------------------------------------------------------------------------------------------------

   # file name accessors for class -- charge density

   def setGridChargeDensityFileName( self, inputFileName, includeWorkingDirectory ): # {

      """Set accessor for grid charge density file name"""

      if includeWorkingDirectory:
         inputFileName = self._getFullWorkingPathName( inputFileName )
      self._gridChargeDensityFileName = inputFileName

   # } end of setGridChargeDensityFileName

   def getGridChargeDensityFileName( self ): # {

      """Get accessor for grid charge density file name"""

      return self._gridChargeDensityFileName

   # } end of getGridChargeDensityFileName

   # -------------------------------------------------------------------------------------------------

   # file name accessors for class -- number density

   def setGridNumberDensityFileName( self, inputFileName, includeWorkingDirectory ): # {

      """Set accessor for grid number density file name"""

      if includeWorkingDirectory:
         inputFileName = self._getFullWorkingPathName( inputFileName )
      self._gridNumberDensityFileName = inputFileName

   # } end of setGridNumberDensityFileName

   def getGridNumberDensityFileName( self ): # {
      """Get accessor for grid potential file name"""
      return self._gridNumberDensityFileName

   # } end of getGridNumberDensityFileName

   # -------------------------------------------------------------------------------------------------

   # accessors for total no. of grid points in object 

   def setTotalGridPoints( self, total ): # {
      """Set accessor for total no. of grid points"""
      self._totalGridPoints = total

   # } end of setTotalGridPoints

   def getTotalGridPoints( self ): # {
      """Get accessor for total no. of grid points"""
      return self._totalGridPoints

   # } end of setTotalGridPoints

   # -------------------------------------------------------------------------------------------------

   # box range accessors for class

   def setBoxRange( self ): # {
      """Set accessor for box range"""
      self._boxRange = [ [ self._xOrigin, self._xOrigin + (self._xGrid*self._xDelta) ],
                         [ self._yOrigin, self._yOrigin + (self._yGrid*self._yDelta) ], 
                         [ self._zOrigin, self._zOrigin + (self._zGrid*self._zDelta) ]
                       ]       
      return None

   # } end of setBoxRange

   def getBoxRange( self ): # {
      """Get accessor for box range"""
      if self._boxRange == None:
         self.setBoxRange()
      return self._boxRange

   # } end of getBoxRange

   # -------------------------------------------------------------------------------------------------

   def getOriginArray( self ): # {
      """Get accessor for origin array"""
      return [ self._xOrigin, self._yOrigin, self._zOrigin ]

   # } end of getOriginArray 

   def getDeltaArray( self ): # {
      """Get accessor for delta array"""
      return [ self._xDelta, self._yDelta, self._zDelta ]

   # } end of getDeltaArray 

   # -------------------------------------------------------------------------------------------------

   def setPotentialExtrema( self ): # {

      """Set extrema for potential"""

      minV    =  1.0e+200
      maxV    = -1.0e+200

       # should be set in loop, but just in case ...
      
      minIndex = -1
      maxIndex = -1
      for ii in range( self.getTotalGridPoints() ):
         if self._potentialList[ii] > maxV:
            maxV     = self._potentialList[ii]
            maxIndex = ii

         if self._potentialList[ii] < minV:
            minV     = self._potentialList[ii]
            minIndex = ii

      self._extrema = [ [ minIndex, minV ], [ maxIndex, maxV ] ]
      return None

   # } end of setPotentialExtrema

   # -------------------------------------------------------------------------------------------------

   def getPotentialExtrema( self ): # {

      """Get extrema for potential"""

      if self._extrema == None:
         self.setPotentialExtrema( )

      return self._extrema

   # } end of setPotentialExtrema

   # -------------------------------------------------------------------------------------------------

   def getBoxStats( self ): # {

      """Get box stats"""

      message = 'BoxStats '
      logReference = self.getLogReference()
      logReference.info( message + ' start: ' )
      if self.getGridPotentialFileName() != None:
         message += ' File: ' + self.getGridPotentialFileName()
      else:
         message += ' unknown origin'
      message += endOfLine

      boxRange = self.getBoxRange()

      message += ' Total=%d' % self.getTotalGridPoints() 
      message += endOfLine

      message += 'Origin   [ %.5f %.5f %.5f ] ' % (self._xOrigin, self._yOrigin, self._zOrigin )
      message += endOfLine

      message += 'Delta [ %.5f %.5f %.5f ] ' % (self._xDelta, self._yDelta, self._zDelta )
      message += endOfLine

      message += 'BoxRange: [ %.4f %.4f ] [ %.4f, %.4f ] [ %.4f %.4f ] ' % ( boxRange[0][0], boxRange[0][1], boxRange[1][0], boxRange[1][1], boxRange[2][0], boxRange[2][1] )
      message += endOfLine

      logReference.info( message + ' start: ' )
      message = ''

      boxExtrema     = self.getPotentialExtrema()
      minCoordinates = self.getGridCoordinates( boxExtrema[0][0] ) 
      maxCoordinates = self.getGridCoordinates( boxExtrema[1][0] ) 
      message       += 'BoxMin: %d [ %.4f  %.4f, %.4f ] V=%.5f ' % ( boxExtrema[0][0], minCoordinates[0], minCoordinates[1], minCoordinates[1], boxExtrema[0][1] )
      message       += endOfLine
      message       += 'BoxMax: %d [ %.4f  %.4f, %.4f ] V=%.5f ' % ( boxExtrema[1][0], maxCoordinates[0], maxCoordinates[1], maxCoordinates[1], boxExtrema[1][1] )
      message += endOfLine

      # check for small difference in extrema -- exit if too small

      if abs( boxExtrema[0][1] - boxExtrema[1][1] ) < 0.1:
         message += 'Difference in extrema is too small: %.4f %.4f' % ( boxExtrema[0][1], boxExtrema[1][1] )
         logReference.info( message + ' start: ' )
         print message
         return
      
#      radius                = 2.5
#      message              += 'Points surrounding minimum:' + endOfLine
#      surroundingPointsList = self.getGridIndicesWithinPointAndRadius( minCoordinates, radius )
#      yJump          = 1.0e+200
#      logReference.info( message + ' start: ' )
#      message = ''
#      for ii in surroundingPointsList:
#         coordinates = self.getGridCoordinates( ii )
#
#         if abs( coordinates[1] - yJump ) > 1.0e-03 and yJump < 1.0e+200:
#            message    += endOfLine
#
#         message    += 'V=%.5f at [ %.4f %.4f %.4f ] index=%d' % ( self._potentialList[ii], coordinates[0], coordinates[1], coordinates[2], ii )
#         yJump       = coordinates[1]
#         message    += endOfLine
#
#      message       += 'Points surrounding maximum:' + endOfLine
#      surroundingPointsList = self.getGridIndicesWithinPointAndRadius( maxCoordinates, radius )
#      for ii in surroundingPointsList:
#         coordinates = self.getGridCoordinates( ii )
#         message    += 'V=%.5f at [ %.4f %.4f %.4f ]' % ( self._potentialList[ii], coordinates[0], coordinates[1], coordinates[2] )
#         message    += endOfLine

      message += 'Sample potential values'
      message += endOfLine
      printMax = 5
      for ii in range( printMax ):
         message      += str(ii)
         message      += '=%.3f ' % self._potentialList[ii]
      message += endOfLine
   
      totalGridPoints = self.getTotalGridPoints()
      for ii in range( totalGridPoints  - printMax, totalGridPoints ):
         message += str(ii) 
         message += '=%.3f ' % self._potentialList[ii]

      return message 

   # } end of getBoxStats

   # -------------------------------------------------------------------------------------------------

   def getDxSortedPotentialList( self ): # {
      """Get grid indices given point:
           xIndex = int( ( x - xOrigin )/xDelta )
           ...

      """

      sortedPotential = self._gridPointList[:]
      sortedPotential.sort( gridPointCompare )

      diagnostics     = 1

      logReference    = self.getLogReference()
      szList          = len( sortedPotential )
      logReference.info( 'In getDxSortedPotentialList ' + str( szList ) + '\n')
      print 'In getDxSortedPotentialList ' + str( szList ) + ' diagn=' + str( diagnostics ) + '\n'

      if diagnostics > 0:
         logReference = self.getLogReference()
         szList       = len( sortedPotential )
         logReference.info( 'In getDxSortedPotentialList ' + str( szList ) + '\n')
         toPrint      = 50
         if szList < toPrint:
            startRange = szList
         else:
            startRange = toPrint

         endRange = szList - toPrint
         if endRange < 0:
            endRange = 0

         for ii in range( 0, startRange ):
            gridPoint       = sortedPotential[ii]
            originalIndex   = gridPoint.getOriginalGridIndex() 
            x_coordinate    = gridPoint.getXCoordinate() 
            y_coordinate    = gridPoint.getYCoordinate() 
            z_coordinate    = gridPoint.getZCoordinate() 
            potential       = gridPoint.getPotential() 
   
            message         = '%5d %5d [ %.4f %.4f %.4f ] %.4f' % ( ii, originalIndex, x_coordinate, y_coordinate, z_coordinate, potential )
            logReference.info( message )

         logReference.info( "\n\nBottom portion" )
         for ii in range( endRange, szList  ):
            gridPoint       = sortedPotential[ii]
            originalIndex   = gridPoint.getOriginalGridIndex() 
            x_coordinate    = gridPoint.getXCoordinate() 
            y_coordinate    = gridPoint.getYCoordinate() 
            z_coordinate    = gridPoint.getZCoordinate() 
            potential       = gridPoint.getPotential() 
   
            message         = '%5d %5d [ %.4f %.4f %.4f ] %.4f' % ( ii, originalIndex, x_coordinate, y_coordinate, z_coordinate, potential )
            logReference.info( message )

      return sortedPotential

   # } end of getRawGridIndicesGivenPoint

   # -------------------------------------------------------------------------------------------------

   def writeGridPotentialDxFile( self, isimPotentialDxFileName ): # {
      """A method to write ISIM a dx-formatted potential
   
      """
    
      # This procedure reads in a ISIM potential grid
   
      # log info
   
      logReference = self.getLogReference()

      parsedPath   = UtilitiesSimTk.parseFileName( __file__ )
      methodName   = parsedPath[1] + '::writeGridPotentialDxFile';
      message      = methodName + '\n   writing file=<' + isimPotentialDxFileName + '>'
      logReference.info( message )
   
      # ----------------------------------------------------------------------------------------------
   
      # open file -- log error if file can not be opened and return None
   
      try:
         isimPotentialFile = open( isimPotentialDxFileName, "w" )
      except IOError:
         logReference.error( methodName + ' file=<' + isimPotentialDxFileName + '> could not be opened.' );
         return None
    
      # ----------------------------------------------------------------------------------------------
   
      # # Data from APBS 0.3.2
      # #
      # # POTENTIAL (kT/e)
      # #
      # object 1 class gridpositions counts 97 97 97
      # origin -1.062995e+02 -1.061945e+02 -1.062995e+02
      # delta 2.214573e+00 0.000000e+00 0.000000e+00
      # delta 0.000000e+00 2.214573e+00 0.000000e+00
      # delta 0.000000e+00 0.000000e+00 2.214573e+00
      # object 2 class gridconnections counts 97 97 97
      # object 3 class array type double rank 0 items 912673 data follows
      # -2.711670e+00 -2.730499e+00 -2.749317e+00
      # -2.768112e+00 -2.786872e+00 -2.805584e+00
        
      # write header comments
      
      isimPotentialFile.write( '#\n# Isim potential\n#\n#\n' )
   
      # write header

      # assume dimensions of grid are cubic -- equal length in all three dimensions

      xCount = int( pow( self.getTotalGridPoints(), 1.0/3.0 ) + machineEpsilon )

      isimPotentialFile.write( 'object 1 class gridpositions counts ' + str( xCount ) + ' ' + str( xCount ) + ' ' + str( xCount ) + '\n' )

      # origin

      originArray  = self.getOriginArray()
      originString = 'origin %.6e %.6e %.6e\n' % (originArray[0], originArray[1], originArray[2] )
      isimPotentialFile.write( originString )

      # delta

      deltaArray  = self.getDeltaArray()
      deltaString = ''
      for ii in range( 3 ):
         deltaString += 'delta '
         for jj in range( 3 ):
            if ii == jj:
               deltaString += '%.6e ' % deltaArray[ii]
            else:
               deltaString += '%.6e ' % 0
         deltaString += '\n'
      isimPotentialFile.write( deltaString )
         
      # grid connections

      gridConnectionsString  = 'object 2 class gridconnections counts %d %d %d\n' % ( xCount, xCount, xCount )
      gridConnectionsString += 'object 3 class array type double rank 0 items ' + str( self.getTotalGridPoints() ) + ' data follows\n'
      isimPotentialFile.write( gridConnectionsString )
      
      # grid potential values
  
      # my understanding is that data are ordered so that z-coordinate increases most quickly, followed by y, followed by x

      dxSortedCoordinates = self.getDxSortedPotentialList()
      numberOfLines       = int( self.getTotalGridPoints()/3.0 + machineEpsilon )
      maxII               = self.getTotalGridPoints() - 2
      ii                  = 0
      while ii < maxII:  
         potentialString  = '%.6e %.6e %.6e\n' % ( dxSortedCoordinates[ii].getPotential(),
                                                   dxSortedCoordinates[ii+1].getPotential(),
                                                   dxSortedCoordinates[ii+2].getPotential() )
         isimPotentialFile.write( potentialString )
         ii += 3

      # stragglers

      ii     = 3*numberOfLines
      offset = self.getTotalGridPoints() - ii
      message = 'Total points=' + str( self.getTotalGridPoints() ) + ' ii=' + str( ii ) + ' offset=' + str( offset )
      logReference.info( message )

      if offset == 2:
         potentialString  = '%.6e %.6e\n' % ( dxSortedCoordinates[ii].getPotential(),
                                             dxSortedCoordinates[ii+1].getPotential() ) 
         isimPotentialFile.write( potentialString )
      elif offset == 1:
         potentialString  = '%.6e\n' % ( dxSortedCoordinates[ii].getPotential() )
         isimPotentialFile.write( potentialString )
 
      # tail

      isimPotentialFile.write( 'attribute "dep" string "positions"\n' )
      isimPotentialFile.write( 'object "regular positions regular connections" class field\n' )
      isimPotentialFile.write( 'component "positions" value 1\n' )
      isimPotentialFile.write( 'component "connections" value 2\n' )
      isimPotentialFile.write( 'component "data" value 3' )

      # close file

      isimPotentialFile.close( )

      # ----------------------------------------------------------------------------------------------

      # log messsage

      message = 'Wrote dx-formatted file to ' + isimPotentialDxFileName + '\n'
      logReference.info( message )

      # ----------------------------------------------------------------------------------------------
   
   # -------------------------------------------------------------------------------------------------

   # } end of writeGridPotentialDxFile

   # -------------------------------------------------------------------------------------------------

   def isValidGridIndex( self, gridIndex ): # {
      """Method returns 0 if grid index in invalid; otherwise it returns 1
   
      """
      # ----------------------------------------------------------------------------------------------

      diagnostics  = 0
      logReference = self.getLogReference()

      # validate input gridIndex

      methodName = 'isValidGridIndex: '
      if gridIndex < 0:
         message = methodName + 'grid index ' + str( gridIndex ) + ' is negative.'
         logReference.error( message )
         return 0 

      if gridIndex >= self.getTotalGridPoints():
         message = methodName + 'grid index ' + str( gridIndex ) + ' is too big: max=' + str( self.getTotalGridPoints() )
         logReference.error( message )
         return 0 

      return 1


   # } end of method isValidGridIndex

   # -------------------------------------------------------------------------------------------------

   def addGridPoint( self, gridIndex, xCoord, yCoord, zCoord, V, chargeDensity, numberDensity ): # {
      """Method adds grid point
   
      """
      # ----------------------------------------------------------------------------------------------

      if self._gridPointList == None:
         self._gridPointList = []

      gridPoint = GridPointSimTk()

      gridPoint.setOriginalGridIndex( gridIndex )
      gridPoint.setXCoordinate( xCoord )
      gridPoint.setYCoordinate( yCoord )
      gridPoint.setZCoordinate( zCoord )
      gridPoint.setPotential( V )
      gridPoint.setChargeDensity( chargeDensity )
      gridPoint.setNumberDensity( numberDensity )

      self._gridPointList.append( gridPoint )

      return 1


   # } end of method addGridPoint

   # -------------------------------------------------------------------------------------------------

   def getGridCoordinates( self, gridIndex ): # {
      """A method to convert grid index to grid coordinates
   
      """
    
      # This procedure converts a grid index to grid coordinates
   
      # log info
   
      diagnostics = 0
      logReference = self.getLogReference()

      # ----------------------------------------------------------------------------------------------

      # validate input gridIndex

      methodName = 'getGridCoordinates: '
      if self.isValidGridIndex( gridIndex ) == 0:
         return None

      # ----------------------------------------------------------------------------------------------

      xCoord             = self._xCoordinateList[gridIndex]
      yCoord             = self._yCoordinateList[gridIndex]
      zCoord             = self._zCoordinateList[gridIndex]

      return [ xCoord, yCoord, zCoord ]

   # } end of method getGridCoordinates

   # -------------------------------------------------------------------------------------------------

   # accessors for comment list

   def setCommentList( self, commentList ): # {

      """Set accessor for comment list"""

      self._commentList = commentList

   # } end of setCommentList

   def getCommentList( self ): # {

      """Get accessor for comment list"""

      return self._commentList

   # } end of getCommentList

   # -------------------------------------------------------------------------------------------------

# } end of class GridSimTk

# ----------------------------------------------------------------------------------------------------

class GridPointSimTk(BaseObjectSimTk.BaseObjectSimTk): # {
   """ grid point """

   # -------------------------------------------------------------------------------------------------

   def __init__( self ): # {

      """Constructor for GridPointSimTk"""
     
      BaseObjectSimTk.BaseObjectSimTk.__init__( self )

      self._originalIndex       = None
      self._xCoord              = None
      self._yCoord              = None
      self._zCoord              = None
      self._potential           = None
      self._chargeDensity       = None
      self._numberDensity       = None

   # } end of __init__(

   # -------------------------------------------------------------------------------------------------

   def setOriginalGridIndex( self, inputOriginalGridIndex ): # {
      """Set accessor for original grid index"""
      self._originalGridIndex = inputOriginalGridIndex

   # } end of setOriginalGridIndex

   def getOriginalGridIndex( self ): # {
      """Get accessor for original grid index"""

      return self._originalGridIndex

   # } end of getOriginalGridIndex

   # -------------------------------------------------------------------------------------------------

   def setPotential( self, inputPotential ): # {
      """Set accessor for potential"""
      self._potential = inputPotential

   # } end of setPotential

   def getPotential( self ): # {
      """Get accessor for potential"""

      return self._potential

   # } end of getPotential

   # -------------------------------------------------------------------------------------------------

   def setChargeDensity( self, inputChargeDensity ): # {
      """Set accessor for chargeDensity"""
      self._chargeDensity = inputChargeDensity

   # } end of setChargeDensity

   def getChargeDensity( self ): # {
      """Get accessor for chargeDensity"""

      return self._chargeDensity

   # } end of getChargeDensity

   # -------------------------------------------------------------------------------------------------

   def setNumberDensity( self, inputNumberDensity ): # {
      """Set accessor for numberDensity"""
      self._numberDensity = inputNumberDensity

   # } end of setNumberDensity

   def getNumberDensity( self ): # {
      """Get accessor for numberDensity"""

      return self._numberDensity

   # } end of getNumberDensity

   # -------------------------------------------------------------------------------------------------

   def setXCoordinate( self, inputXCoordinate ): # {
      """Set accessor for x-coordinate"""
      self._xCoordinate = inputXCoordinate

   # } end of setXCoordinate

   def getXCoordinate( self ): # {
      """Get accessor for x-coordinate"""

      return self._xCoordinate

   # } end of getXCoordinate

   # -------------------------------------------------------------------------------------------------

   def setYCoordinate( self, inputYCoordinate ): # {
      """Set accessor for y-coordinate """
      self._yCoordinate = inputYCoordinate

   # } end of setYCoordinate

   def getYCoordinate( self ): # {
      """Get accessor for y-coordinate"""

      return self._yCoordinate

   # } end of getYCoordinate

   # -------------------------------------------------------------------------------------------------

   def setZCoordinate( self, inputZCoordinate ): # {
      """Set accessor for z-coordinate"""
      self._zCoordinate = inputZCoordinate

   # } end of setZCoordinate

   def getZCoordinate( self ): # {
      """Get accessor for z-coordinate"""

      return self._zCoordinate

   # } end of getZCoordinate

   # -------------------------------------------------------------------------------------------------

   def getCoordinateString( self ): # {
      """Get string containing grid coordinates"""
      return '[ %12.6f %12.6f %12.6f ]' % ( self._xCoordinate, self._yCoordinate, self._zCoordinate )

   # } end of getCoordinateString

   def getZCoordinate( self ): # {
      """Get accessor for z-coordinate"""

      return self._zCoordinate

   # } end of getZCoordinate


   # -------------------------------------------------------------------------------------------------

# } end of class GridPointSimTk

# -------------------------------------------------------------------------------------------------

printCount = 0
def gridPointCompare( gridPoint1, gridPoint2 ):

   """Comparison method used for sorting grid points into dx-order"""

   returnValue = 0

   # message in apbs python script that reads dx files:

   # WARNING!  THIS PROGRAM ASSUMES YOU HAVE THE DATA ORDERED SUCH
   # THAT Z INCREASES MOST QUICKLY, Y INCREASES NEXT MOST QUICKLY,
   # AND X INCREASES MOST SLOWLY.

   # see comments at top of this file for how data is indexed

   if abs( gridPoint1.getZCoordinate() - gridPoint2.getZCoordinate() ) < machineEpsilon:  
      if abs( gridPoint1.getYCoordinate() - gridPoint2.getYCoordinate() ) < machineEpsilon:  

         if gridPoint1.getXCoordinate() > gridPoint2.getXCoordinate():
            returnValue = 1
         elif gridPoint1.getXCoordinate() < gridPoint2.getXCoordinate():
            returnValue = -1
         else:
            returnValue = 0

      elif gridPoint1.getYCoordinate() > gridPoint2.getYCoordinate():
         returnValue = 1
      else:
         returnValue = -1

   elif gridPoint1.getZCoordinate() > gridPoint2.getZCoordinate():
      returnValue = 1
   else:
      returnValue = -1

   # diagnostics

   global printCount
   if printCount < 100:
       point1 = gridPoint1.getCoordinateString()
       point2 = gridPoint2.getCoordinateString()
       if gridPoint1.getXCoordinate() > gridPoint2.getXCoordinate():
           diffX = 1
       else:
           diffX = 0

       if gridPoint1.getYCoordinate() > gridPoint2.getYCoordinate():
           diffY = 1
       else:
           diffY = 0

       if gridPoint1.getZCoordinate() > gridPoint2.getZCoordinate():
           diffZ = 1
       else:
           diffZ = 0

       print 'Cmp=' + str( returnValue ) + ' ' +  \
             point1 + ' ' + point2 + ' [' + str(diffX) + ' ' + str( diffY ) + \
             ' ' + str( diffZ ) + ']'
       printCount += 1

   return returnValue

# } end of method gridPointCompare