def distance(v0, v1):
    from math import sqrt

    return sqrt((v1[0]-v0[0])**2 + (v1[1]-v0[1])**2)

def length(v):
    from math import sqrt

    return sqrt(v[0]**2 + v[1]**2)

def inner_product(v0, v1):
    from math import sqrt

    #FIXME: Obsolete
    try:
        i = sqrt(v0[0]*v1[0] + v0[1]*v1[1])
    except:
        print v0, v1
        import sys; sys.exit()

    return i    

def angle(v):
    """Compute angle between e1 (the unit vector in the x-direction)
    and the specified vector
    """

    from math import acos, pi

    l = length(v)
    v1 = v[0]/l
    v2 = v[1]/l
        
    theta = acos(v1)

    if v2 < 0:
        #Quadrant 3 or 4
        theta = 2*pi-theta

    return theta


def anglediff(v0, v1):
    """Compute difference between angle of vector x0, y0 and x1, y1.
    This is used for determining the ordering of vertices,
    e.g. for checking if they are counter clockwise. 
    
    Always return a positive value
    """
        
    from math import pi
    
    a0 = angle(v0)
    a1 = angle(v1)

    #Ensure that difference will be positive
    if a0 < a1:
        a0 += 2*pi
            
    return a0-a1


 
def compile(FNs=None, CC=None, LD = None, SFLAG = None, verbose = 1):
  """compile(FNs=None, CC=None, LD = None, SFLAG = None):
  
     Compile FN(s) using compiler CC (e.g. mpicc), 
     Loader LD and shared flag SFLAG.
     If CC is absent use default compiler dependent on platform
     if LD is absent CC is used.
     if SFLAG is absent platform default is used
     FNs can be either one filename or a list of filenames
     In the latter case, the first will be used to name so file.
  """
  import os, string, sys, types
  
  # Input check
  #
  assert not FNs is None, "No filename provided"

  if not type(FNs) == types.ListType:
    FNs = [FNs]


  libext = 'so' #Default extension (Unix)
  libs = ''
  version = sys.version[:3]
  
  # Determine platform and compiler
  #
  if sys.platform == 'sunos5':  #Solaris
    if CC:
      compiler = CC
    else:  
      compiler = 'gcc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'G'
      
  elif sys.platform == 'osf1V5':  #Compaq AlphaServer
    if CC:
      compiler = CC
    else:  
      compiler = 'cc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'shared'    
      
  elif sys.platform == 'linux2':  #Linux
    if CC:
      compiler = CC
    else:  
      compiler = 'gcc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'shared'    
      
  elif sys.platform == 'darwin':  #Mac OS X:
    if CC:
      compiler = CC
    else:  
      compiler = 'cc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'bundle -flat_namespace -undefined suppress'

  elif sys.platform == 'cygwin':  #Cygwin (compilation same as linux)
    if CC:
      compiler = CC
    else:  
      compiler = 'gcc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'shared'
    libext = 'dll'
    libs = '/lib/python%s/config/libpython%s.dll.a' %(version,version)
      
  elif sys.platform == 'win32':  #Windows
    if CC:
      compiler = CC
    else:  
      compiler = 'gcc'
    if LD:
      loader = LD
    else:  
      loader = compiler
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'shared'
    libext = 'dll'

    v = version.replace('.','')
    dllfilename = 'python%s.dll' %(v)
    libs = os.path.join(sys.exec_prefix,dllfilename)
      
      
  else:
    if verbose: print "Unrecognised platform %s - revert to default"\
                %sys.platform
    
    if CC:
      compiler = CC
    else:  
      compiler = 'cc'
    if LD:
      loader = LD
    else:  
      loader = 'ld'
    if SFLAG:
      sharedflag = SFLAG
    else:  
      sharedflag = 'G'

            
  # Find location of include files
  #
  if sys.platform == 'win32':  #Windows
    include = os.path.join(sys.exec_prefix, 'include')    
  else:  
    include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
                           'python'+version)

  # Check existence of Python.h
  #
  headerfile = include + os.sep +'Python.h'
  try:
    open(headerfile, 'r')
  except:
    raise """Did not find Python header file %s.
    Make sure files for Python C-extensions are installed. 
    In debian linux, for example, you need to install a
    package called something like python2.3-dev""" %headerfile
  
  
  # Check filename(s)
  #
  object_files = ''
  for FN in FNs:        
    root, ext = os.path.splitext(FN)
    if ext == '':
      FN = FN + '.c'
    elif ext.lower() != '.c':
      raise Exception, "Unrecognised extension: " + FN
    
    try:
      open(FN,'r')
    except:    
      raise Exception, "Could not open: " + FN

    if not object_files: root1 = root  #Remember first filename        
    object_files += root + '.o '  
  
  
    # Compile
    #
    
    s = "%s -c %s -I%s -o %s.o -Wall -O" %(compiler, FN, include, root)
    if verbose:
      print s
    else:
      s = s + ' 2> /dev/null' #Suppress errors
  
    try:
      err = os.system(s)
      if err != 0:
          raise 'Attempting to compile %s failed - please try manually' %FN 
    except:
      raise 'Could not compile %s - please try manually' %FN  

  
  # Make shared library (*.so or *.dll)
  
  s = "%s -%s %s -o %s.%s %s -lm" %(loader, sharedflag, object_files, root1, libext, libs)
  if verbose:
    print s
  else:
    s = s + ' 2> /dev/null' #Suppress warnings
  
  try:  
    err=os.system(s)
    if err != 0:	
	raise 'Atempting to link %s failed - please try manually' %root1     
  except:
    raise 'Could not link %s - please try manually' %root1
    



def can_use_C_extension(filename):
    """Determine whether specified C-extension
    can and should be used.
    """

    from config import use_extensions

    from os.path import splitext

    root, ext = splitext(filename)
    
    C=False
    if use_extensions:
        try:
            s = 'import %s' %root
            #print s
            exec(s)
        except:
            try:
                open(filename)
            except:
                msg = 'C extension %s cannot be opened' %filename
                print msg                
            else:    
                print '------- Trying to compile c-extension %s' %filename
	    
                compile(filename)		    
                try:
                    compile(filename)
                except:
                    print 'WARNING: Could not compile C-extension %s'\
                          %filename
                else:
                    try:
                        exec('import %s' %root)
                    except:
                        msg = 'C extension %s seems to compile OK, '
                        msg += 'but it can still not be imported.'
                        raise msg
                    else:
                        C=True
        else:
            C=True
            
    if not C:
        pass
        print 'NOTICE: C-extension %s not used' %filename

    return C





####################################################################
#Python versions of function that are also implemented in util_ext.c
#

def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
    """
    """
    
    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
    a /= det

    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
    b /= det            

    return a, b



def rotate_python(q, normal, direction = 1):
    """Rotate the momentum component q (q[1], q[2])
    from x,y coordinates to coordinates based on normal vector.

    If direction is negative the rotation is inverted.
    
    Input vector is preserved

    This function is specific to the shallow water wave equation
    """

    #FIXME: Needs to be tested

    from Numeric import zeros, Float
    
    assert len(q) == 3,\
           'Vector of conserved quantities must have length 3'\
           'for 2D shallow water equation'
    try:
        l = len(normal)
    except:
        raise 'Normal vector must be an Numeric array'

    assert l == 2, 'Normal vector must have 2 components'

 
    n1 = normal[0]
    n2 = normal[1]    
    
    r = zeros(len(q), Float) #Rotated quantities
    r[0] = q[0]              #First quantity, height, is not rotated

    if direction == -1:
        n2 = -n2


    r[1] =  n1*q[1] + n2*q[2]
    r[2] = -n2*q[1] + n1*q[2]
    
    return r




##############################################
#Initialise module

import util
if util.can_use_C_extension('util_ext.c'):
    from util_ext import gradient, rotate
else:
    gradient = gradient_python
    rotate = rotate_python




if __name__ == "__main__":
    pass