source: inundation/ga/storm_surge/pyvolution/util.py @ 628

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

    from math import acos, pi, sqrt
    from Numeric import sum, array

    l = sqrt( sum (array(v)**2))
    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 mean(x):
    from Numeric import sum
    return sum(x)/len(x)




#FIXME: Maybe move to util as this is quite general
class File_function:
    """Read time series from file and return a callable object:
    f(t,x,y)
    which will return interpolated values based on the input file.

    The file format is assumed to be either two fields separated by a comma:

        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...

    or four comma separated fields

        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...

    In either case, the callable obejct will return a tuple of
    interpolated values, one each value stated in the file.

    
    E.g

      31/08/04 00:00:00, 1.328223 0 0
      31/08/04 00:15:00, 1.292912 0 0

    will provide a time dependent function f(t,x=None,y=None),
    ignoring x and y, which returns three values per call.

    
    NOTE: At this stage the function is assumed to depend on
    time only, i.e  no spatial dependency!!!!!
    When that is need we can use the least_squares interpolation.
    """

    
    def __init__(self, filename, domain=None):
        """Initialise callable object from a file.
        See docstring for class File_function

        If domain is specified, model time (domain,starttime)
        will be checked and possibly modified

        ALl times are assumed to be in UTC
        """

        import time, calendar
        from Numeric import array
        from config import time_format

        assert type(filename) == type(''),\
               'First argument to File_function must be a string'


        try:
            fid = open(filename)
        except Exception, e:
            msg = 'File "%s" could not be opened: Error="%s"'\
                  %(filename, e)
            raise msg


        line = fid.readline()
        fid.close()
        fields = line.split(',')
        msg = 'File %s must have the format date, value0 value1 value2 ...'
        msg += ' or date, x, y, value0 value1 value2 ...'
        mode = len(fields)
        assert mode in [2,4], msg

        try:
            starttime = calendar.timegm(time.strptime(fields[0], time_format))
        except ValueError:    
            msg = 'First field in file %s must be' %filename
            msg += ' date-time with format %s.\n' %time_format
            msg += 'I got %s instead.' %fields[0] 
            raise msg


        #Split values
        values = []
        for value in fields[mode-1].split():
            values.append(float(value))

        q = array(values)
        
        msg = 'ERROR: File must contain at least one independent value'
        assert len(q.shape) == 1, msg
            
        self.number_of_values = len(q)

        self.filename = filename
        self.starttime = starttime


        if domain is not None:
            if domain.starttime is None:
                domain.starttime = starttime
            else:
                msg = 'WARNING: Start time as specified in domain (%s)'\
                      %domain.starttime
                msg += ' is earlier than the starttime of file %s: %s.'\
                         %(self.filename, self.starttime)
                msg += 'Modifying starttime accordingly.'
                if self.starttime > domain.starttime:
                    #FIXME: Print depending on some verbosity setting
                    #print msg
                    domain.starttime = self.starttime #Modifying model time

        if mode == 2:        
            self.read_times() #Now read all times in.
        else:
            raise 'x,y dependecy not yet implemented'

        
    def read_times(self):
        """Read time and values 
        """
        from Numeric import zeros, Float, alltrue
        from config import time_format
        import time, calendar
        
        fid = open(self.filename)        
        lines = fid.readlines()
        fid.close()
        
        N = len(lines)
        d = self.number_of_values

        T = zeros(N, Float)       #Time
        Q = zeros((N, d), Float)  #Values
        
        for i, line in enumerate(lines):
            fields = line.split(',')
            realtime = calendar.timegm(time.strptime(fields[0], time_format))

            T[i] = realtime - self.starttime

            for j, value in enumerate(fields[1].split()):
                Q[i, j] = float(value)

        msg = 'File %s must list time as a monotonuosly ' %self.filename
        msg += 'increasing sequence'
        assert alltrue( T[1:] - T[:-1] > 0 ), msg

        self.T = T     #Time
        self.Q = Q     #Values
        self.index = 0 #Initial index

        
    def __repr__(self):
        return 'File function'

        

    def __call__(self, t, x=None, y=None):
        """Evaluate f(t,x,y)

        FIXME: x, y dependency not yet implemented
        """

        from math import pi, cos, sin, sqrt
        
        
        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
              %(self.filename, self.T[0], self.T[1], t)
        if t < self.T[0]: raise msg
        if t > self.T[-1]: raise msg        

        oldindex = self.index

        #Find slot
        while t > self.T[self.index]: self.index += 1
        while t < self.T[self.index]: self.index -= 1

        #t is now between index and index+1
        ratio = (t - self.T[self.index])/\
                (self.T[self.index+1] - self.T[self.index])

        #Compute interpolated values for values
        q = self.Q[self.index,:] +\
            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 

        #Return vector of interpolated values
        return q
    



####################################################################
#Python versions of function that are also implemented in util_gateway.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





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

import compile
if compile.can_use_C_extension('util_ext.c'):
    from util_ext import gradient
else:
    gradient = gradient_python


if __name__ == "__main__":
    pass