source: lmdz_wrf/trunk/tools/geometry_tools.py @ 2728

# Python tools to manage netCDF files.
# L. Fita, CIMA. March 2019
# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
#
# pyNCplot and its component geometry_tools.py comes with ABSOLUTELY NO WARRANTY. 
# This work is licendes under a Creative Commons 
#   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
#
## Script for geometry calculations and operations as well as definition of different 
###    standard objects and shapes

import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import os
import generic_tools as gen
import numpy.ma as ma
import module_ForSci as fsci

errormsg = 'ERROR -- error -- ERROR -- error'
infmsg = 'INFORMATION -- information -- INFORMATION -- information'

####### Contents:
# add_secpolygon_list: Function to add a range of points of a polygon into a list
# angle_vectors2D: Angle between two vectors with sign
# cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the [x/y]-axis
# cut_[x/y]polygon: Function to cut a polygon from a given value of the [x/y]-axis
# deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
# displace_objdic_2D: Function to displace 2D plain the vertices of all polygons of an object
# dist_points: Function to provide the distance between two points
# join_circ_sec: Function to join aa series of points by circular segments
# join_circ_sec_rand: Function to join aa series of points by circular segments with 
#   random perturbations
# max_coords_poly: Function to provide the extremes of the coordinates of a polygon
# mirror_polygon: Function to reflex a polygon for a given axis
# mod_vec: Function to compute the module of a vector
# position_sphere: Function to tranform fom a point in lon, lat deg coordinates to 
#   cartesian coordinates over an sphere
# read_join_poly: Function to read an ASCII file with the combination of polygons
# rm_consecpt_polygon: Function to remove consecutive same point of a polygon
# rotate_2D: Function to rotate a vector by a certain angle in the plain
# rotate_objdic_2D: Function to rotate 2D plain the vertices of all polygons of an object
# rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon
# rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a 
#   certain angle in the plain
# rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y 
#   coordinates by a certain angle in the plain
# spheric_line: Function to transform a series of locations in lon, lat coordinates 
#   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
# val_consec_between: Function to provide if a given value is between two consecutive ones
# write_join_poly: Function to write an ASCII file with the combination of polygons

## Shapes/objects
# circ_sec: Function union of point A and B by a section of a circle
# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
# p_angle_triangle: Function to draw a triangle by an initial point and two 
#   consecutive angles and the first length of face. The third angle and 2 and 3rd 
#   face will be computed accordingly the provided values
# p_doubleArrow: Function to provide an arrow with double lines
# p_circle: Function to get a polygon of a circle
# p_cross_width: Function to draw a cross with arms with a given width and an angle
# p_prism: Function to get a polygon prism
# p_reg_polygon: Function to provide a regular polygon of Nv vertices
# p_reg_star: Function to provide a regular star of Nv vertices
# p_sinusoide: Function to get coordinates of a sinusoidal curve
# p_square: Function to get a polygon square
# p_spiral: Function to provide a polygon of an Archimedean spiral
# p_triangle: Function to provide the polygon of a triangle from its 3 vertices
# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates

## Plotting
# draw_secs: Function to draw an object according to its dictionary
# paint_filled: Function to draw an object filling given sections
# plot_sphere: Function to plot an sphere and determine which standard lines will be 
#   also drawn

def deg_deci(angle):
    """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
      angle: list of [deg, minute, sec] to pass
    >>> deg_deci([41., 58., 34.])
    0.732621346072
    """
    fname = 'deg_deci'

    deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600.

    if angle[0] < 0.: deg = -deg*np.pi/180.
    else: deg = deg*np.pi/180.

    return deg

def position_sphere(radii, alpha, beta):
    """ Function to tranform fom a point in lon, lat deg coordinates to cartesian  
          coordinates over an sphere
      radii: radii of the sphere
      alpha: longitude of the point
      beta: latitude of the point
    >>> position_sphere(10., 30., 45.)
    (0.81031678432964027, -5.1903473778327376, 8.5090352453411846
    """
    fname = 'position_sphere'

    xpt = radii*np.cos(beta)*np.cos(alpha)
    ypt = radii*np.cos(beta)*np.sin(alpha)
    zpt = radii*np.sin(beta)

    return xpt, ypt, zpt

def spheric_line(radii,lon,lat):
    """ Function to transform a series of locations in lon, lat coordinates to x,y,z 
          over an 3D space
      radii: radius of the sphere
      lon: array of angles along longitudes
      lat: array of angles along latitudes
    """
    fname = 'spheric_line'

    Lint = lon.shape[0]
    coords = np.zeros((Lint,3), dtype=np.float)

    for iv in range(Lint):
        coords[iv,:] = position_sphere(radii, lon[iv], lat[iv])

    return coords

def rotate_2D(vector, angle):
    """ Function to rotate a vector by a certain angle in the plain
      vector= vector to rotate [y, x]
      angle= angle to rotate [rad]
    >>> rotate_2D(np.array([1.,0.]), np.pi/4.)
    [ 0.70710678 -0.70710678]
    """
    fname = 'rotate_2D'

    rotmat = np.zeros((2,2), dtype=np.float)

    rotmat[0,0] = np.cos(angle)
    rotmat[0,1] = -np.sin(angle)
    rotmat[1,0] = np.sin(angle)
    rotmat[1,1] = np.cos(angle)

    rotvector = np.zeros((2), dtype=np.float)

    vecv = np.zeros((2), dtype=np.float)

    # Unifying vector
    modvec = vector[0]**2+vector[1]**2
    if modvec != 0 and vector[0] != gen.fillValue: 
        vecv[0] = vector[1]/modvec
        vecv[1] = vector[0]/modvec

        rotvec = np.matmul(rotmat, vecv)
        rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec)

        rotvector[0] = rotvec[1]*modvec
        rotvector[1] = rotvec[0]*modvec
    else:
        rotvector = vector + 0.

    return rotvector

def rotate_polygon_2D(vectors, angle):
    """ Function to rotate 2D plain the vertices of a polygon
      line= matrix of vectors to rotate
      angle= angle to rotate [rad]
    >>> square = np.zeros((4,2), dtype=np.float)
    >>> square[0,:] = [-0.5,-0.5]
    >>> square[1,:] = [0.5,-0.5]
    >>> square[2,:] = [0.5,0.5]
    >>> square[3,:] = [-0.5,0.5]
    >>> rotate_polygon_2D(square, np.pi/4.)
    [[-0.70710678  0.        ]
     [ 0.         -0.70710678]
     [ 0.70710678  0.        ]
     [ 0.          0.70710678]]
    """
    fname = 'rotate_polygon_2D'

    rotvecs = np.zeros(vectors.shape, dtype=np.float)

    mavec = False
    if type(vectors) == type(gen.mamat):
        mavec = True
        vectors = ma.filled(vectors,gen.fillValueF)

    Nvecs = vectors.shape[0]
    for iv in range(Nvecs):
        rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)

    if mavec: 
        rotvecs = ma.masked_equal(rotvecs, gen.fillValueF)

    return rotvecs

def displace_objdic_2D(objdic, distance):
    """ Function to displace 2D plain the vertices of all polygons of an object
      objdic= dictionary with all the polygons of the object
      distance= distance to displace [ydist, xdist]
    """
    fname = 'displace_objdic_2D'

    disobjdic = dict(objdic)

    for secn in objdic.keys():
        objv = objdic[secn]
        vectors = objv[0]
        lt = objv[1]
        lc = objv[2]
        lw = objv[3]

        disvecs = np.zeros(vectors.shape, dtype=np.float)    
        disvecs = vectors + distance
        disobjdic[secn] = [disvecs, lt, lc, lw]

    return disobjdic

def rotate_objdic_2D(objdic, angle):
    """ Function to rotate 2D plain the vertices of all polygons of an object
      objdic= dictionary with all the polygons of the object
      angle= angle to rotate [rad]
    """
    fname = 'rotate_objdic_2D'

    rotobjdic = dict(objdic)

    for secn in objdic.keys():
        objv = objdic[secn]
        vectors = objv[0]
        lt = objv[1]
        lc = objv[2]
        lw = objv[3]

        rotvecs = np.zeros(vectors.shape, dtype=np.float)

        Nvecs = vectors.shape[0]
        for iv in range(Nvecs):
            rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
        rotobjdic[secn] = [rotvecs, lt, lc, lw]

    return rotobjdic

def rotate_line2D(line, angle):
    """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain 
          angle in the plain
      line= line to rotate as couple of points [[y0,x0], [y1,x1]]
      angle= angle to rotate [rad]
    >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.)
    [[ 0.          0.        ]
     [0.70710678  -0.70710678]]
    """
    fname = 'rotate_2D'

    rotline = np.zeros((2,2), dtype=np.float)
    rotline[0,:] = rotate_2D(line[0,:], angle)
    rotline[1,:] = rotate_2D(line[1,:], angle)

    return rotline

def rotate_lines2D(lines, angle):
    """ Function to rotate multiple lines given by mulitple pars of x,y coordinates  
          by a certain angle in the plain
      line= matrix of N couples of points [N, [y0,x0], [y1,x1]]
      angle= angle to rotate [rad]
    >>> square = np.zeros((4,2,2), dtype=np.float)
    >>> square[0,0,:] = [-0.5,-0.5]
    >>> square[0,1,:] = [0.5,-0.5]
    >>> square[1,0,:] = [0.5,-0.5]
    >>> square[1,1,:] = [0.5,0.5]
    >>> square[2,0,:] = [0.5,0.5]
    >>> square[2,1,:] = [-0.5,0.5]
    >>> square[3,0,:] = [-0.5,0.5]
    >>> square[3,1,:] = [-0.5,-0.5]
    >>> rotate_lines2D(square, np.pi/4.)
    [[[-0.70710678  0.        ]
      [ 0.         -0.70710678]]

     [[ 0.         -0.70710678]
      [ 0.70710678  0.        ]]

     [[ 0.70710678  0.        ]
      [ 0.          0.70710678]]

     [[ 0.          0.70710678]
      [-0.70710678  0.        ]]]
    """
    fname = 'rotate_lines2D'

    rotlines = np.zeros(lines.shape, dtype=np.float)

    Nlines = lines.shape[0]
    for il in range(Nlines):
        line = np.zeros((2,2), dtype=np.float)
        line[0,:] = lines[il,0,:]
        line[1,:] = lines[il,1,:]

        rotlines[il,:,:] = rotate_line2D(line, angle)

    return rotlines

def dist_points(ptA, ptB):
    """ Function to provide the distance between two points
      ptA: coordinates of the point A [yA, xA]
      ptB: coordinates of the point B [yB, xB]
    >>> dist_points([1.,1.], [-1.,-1.])
    2.82842712475
    """
    fname = 'dist_points'

    dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2)

    return dist

def mod_vec(vec):
    """ Function to compute the module of a vector
      vec: vector [y, x]
    >>> mod_vec([1., 1.])
    1.41421356237
    """
    fname = 'mod_vec'

    v = np.array(vec, dtype=np.float)
    vv = v*v
    mod = np.sqrt(np.sum(vv[:]))

    return mod

def angle_vectors2D(veca, vecb):
    """ Angle between two vectors with sign
        FROM: https://stackoverflow.com/questions/5188561/signed-angle-between-two-3d-vectors-with-same-origin-within-the-same-plane
      veca: angle A [ya, xa]
      vecb: angle B [yb, xb]
      NOTE: angle from A to B
    >>> angle_vectors2D([1.,0.], [0.,1.])
    1.57079632679
    >>> angle_vectors2D([0.,1.], [1.,0.])
    -1.57079632679
    """
    fname = 'angle_vectors2D'

    v1 = np.array(veca, dtype=np.float)
    v2 = np.array(vecb, dtype=np.float)

    moda = mod_vec(v1)
    modb = mod_vec(v2)
    modab = mod_vec(v1*v2)

    vc = np.cross(v1,v2)
    theta = np.arcsin(vc/(moda*modb))

    # Without sign
    #alpha = np.arccos(modab/(moda*modb))

    return theta

def max_coords_poly(polygon):
    """ Function to provide the extremes of the coordinates of a polygon
      polygon: coordinates [Nvertexs, 2] of a polygon
    >>> square = np.zeros((4,2), dtype=np.float)
    >>> square[0,:] = [-0.5,-0.5]
    >>> square[1,:] = [0.5,-0.5]
    >>> square[2,:] = [0.5,0.5]
    >>> square[3,:] = [-0.5,0.5]
    >>> max_coords_poly(square)
    [-0.5, 0.5], [-0.5, 0.5], [0.5, 0.5], 0.5
    """
    fname = 'max_coords_poly'

    # x-coordinate min/max
    nx = np.min(polygon[:,1])
    xx = np.max(polygon[:,1])

    # y-coordinate min/max
    ny = np.min(polygon[:,0])
    xy = np.max(polygon[:,0])

    # x/y-coordinate maximum of absolute values
    axx = np.max(np.abs(polygon[:,1]))
    ayx = np.max(np.abs(polygon[:,0]))

    # absolute maximum
    xyx = np.max([axx, ayx])

    return [nx, xx], [ny, xy], [ayx, axx], xyx

def mirror_polygon(polygon,axis):
    """ Function to reflex a polygon for a given axis
      polygon: polygon to mirror
      axis: axis at which mirror is located ('x' or 'y')
    """
    fname = 'mirror_polygon'

    reflex = np.zeros(polygon.shape, dtype=np.float)

    N = polygon.shape[0]
    if axis == 'x':
        for iv in range(N):
            reflex[iv,:] = [-polygon[iv,0], polygon[iv,1]]
    elif axis == 'y':
        for iv in range(N):
            reflex[iv,:] = [polygon[iv,0], -polygon[iv,1]]

    return reflex

def join_circ_sec(points, radfrac=3., arc='short', side='left', N=200):
    """ Function to join aa series of points by circular segments
      points: main points of the island (clockwise ordered, to be joined by circular 
        segments of radii as the radfrac factor of the distance between 
        consecutive points)
      radfrac: multiplicative factor of the distance between consecutive points to 
        draw the circular segment (3., default)
      arc: type of arc ('short', default)
      pos: position of arc ('left', default)
      N: number of points (200, default)
    """
    fname = 'join_circ_sec'

    jcirc_sec = np.ones((N,2), dtype=np.float)

    # main points
    lpoints = list(points)
    Npts = len(lpoints)
    Np = int(N/(Npts+1))
    for ip in range(Npts-1):
        p1 = lpoints[ip]
        p2 = lpoints[ip+1]
        dps = dist_points(p1, p2)
        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1,p2,dps*radfrac, arc, side, Np)

    Np2 = N - (Npts-1)*Np
    p1 = lpoints[Npts-1]
    p2 = lpoints[0]
    dps = dist_points(p1, p2)
    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, side, Np2)

    return jcirc_sec

def join_circ_sec_rand(points, radfrac=3., Lrand=0.1, arc='short', pos='left', N=200):
    """ Function to join aa series of points by circular segments with random 
        perturbations
      points: main points of the island (clockwise ordered, to be joined by circular 
        segments of radii as the radfrac factor of the distance between 
        consecutive points)
      radfrac: multiplicative factor of the distance between consecutive points to 
        draw the circular segment (3., default)
      Lrand: maximum length of the random perturbation to be added perpendicularly to 
        the direction of the union line between points (0.1, default)
      arc: type of arc ('short', default)
      pos: position of arc ('left', default)
      N: number of points (200, default)
    """
    import random
    fname = 'join_circ_sec_rand'

    jcirc_sec = np.ones((N,2), dtype=np.float)

    # main points
    lpoints = list(points)
    Npts = len(lpoints)
    Np = int(N/(Npts+1))
    for ip in range(Npts-1):
        p1 = lpoints[ip]
        p2 = lpoints[ip+1]
        dps = dist_points(p1, p2)
        angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, arc, pos, Np)
        drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
        for iip in range(Np*ip,Np*(ip+1)):
            jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.)

    Np2 = N - (Npts-1)*Np
    p1 = lpoints[Npts-1]
    p2 = lpoints[0]
    dps = dist_points(p1, p2)
    angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, pos, Np2)
    drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
    for iip in range(Np*(Npts-1),N):
        jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.)

    return jcirc_sec

def write_join_poly(polys, flname='join_polygons.dat'):
    """ Function to write an ASCII file with the combination of polygons
      polys: dictionary with the names of the different polygons
      flname: name of the ASCII file
    """
    fname = 'write_join_poly'

    of = open(flname, 'w')

    for polyn in polys.keys():
        vertices = polys[polyn]
        Npts = vertices.shape[0]
        for ip in range(Npts):
            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')

    of.close()

    return

def read_join_poly(flname='join_polygons.dat'):
    """ Function to read an ASCII file with the combination of polygons
      flname: name of the ASCII file
    """
    fname = 'read_join_poly'

    of = open(flname, 'r')

    polys = {}
    polyn = ''
    poly = []
    for line in of:
        if len(line) > 1: 
            linevals = line.replace('\n','').split(' ')
            if polyn != linevals[0]:
                if len(poly) > 1:
                    polys[polyn] = np.array(poly)
                polyn = linevals[0]
                poly = []
                poly.append([np.float(linevals[2]), np.float(linevals[1])])
            else:
                poly.append([np.float(linevals[2]), np.float(linevals[1])])

    of.close()
    polys[polyn] = np.array(poly)

    return polys

def val_consec_between(valA, valB, val):
    """ Function to provide if a given value is between two consecutive ones
      valA: first value
      valB: second value
      val: value to determine if it is between
      >>> val_consec_between(0.5,1.5,0.8)
      True
      >>> val_consec_between(0.5,1.5.,-0.8)
      False
      >>> val_consec_between(0.5,1.5,0.5)
      True
      >>> val_consec_between(-1.58, -1.4, -1.5)
      True
      >>> val_consec_between(-1.48747753212, -1.57383530044, -1.5)
      False
    """
    fname = 'val_consec_between'

    btw = False
    diffA = valA - val
    diffB = valB - val
    absdA = np.abs(diffA)
    absdB = np.abs(diffB)
    #if (diffA/absdA)* (diffB/absdB) < 0.: btw = True
#    if valA < 0. and valB < 0. and val < 0.:
#        if (valA >= val and valB < val) or (valA > val and valB <= val): btw =True
#    else:
#        if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True
    if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True

    return btw

def add_secpolygon_list(listv, iip, eep, polygon):
    """ Function to add a range of points of a polygon into a list
      listv: list into which add values of the polygon
      iip: initial value of the range
      eep: ending value of the range
      polygon: array with the points of the polygon
    """
    fname = 'add_secpolygon_list'

    if eep > iip:
        for ip in range(iip,eep): listv.append(polygon[ip,:])
    else:
        for ip in range(iip,eep,-1): listv.append(polygon[ip,:])

    return

def rm_consecpt_polygon(polygon):
    """ Function to remove consecutive same point of a polygon
      poly: polygon
    >>> poly = np.ones((5,2), dtype=np.float)
    >>> poly[2,:] = [2., 1.]
    rm_consecpt_polygon(poly)
    [[ 1.  1.]
     [ 2.  1.]
     [ 1.  1.]]
    """
    fname = 'rm_consecpt_polygon'

    newpolygon = []
    prevpt = polygon[0,:]
    newpolygon.append(prevpt)
    for ip in range(1,polygon.shape[0]):
        if polygon[ip,0] != prevpt[0] or polygon[ip,1] != prevpt[1]:
            prevpt = polygon[ip,:]
            newpolygon.append(prevpt)

    newpolygon = np.array(newpolygon)

    return newpolygon

def cut_ypolygon(polygon, yval, keep='below', Nadd=20):
    """ Function to cut a polygon from a given value of the y-axis
      polygon: polygon to cut
      yval: value to use to cut the polygon
      keep: part to keep from the height ('below', default)
         'below': below the height
         'above': above the height
      Nadd: additional points to add to draw the line (20, default)
    """
    fname = 'cut_ypolygon'

    N = polygon.shape[0]
    availkeeps = ['below', 'above']

    if not gen.searchInlist(availkeeps, keep):
        print errormsg
        print '  ' + fname + ": wring keep '" + keep + "' value !!"
        print '    available ones:', availkeeps
        quit(-1)

    ipt = None
    ept = None

    # There might be more than 1 cut...
    Ncuts = 0
    icut = []
    ecut = []
    ipt = []
    ept = []

    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
      type(gen.mamat.mask[1]):
        # Assuming clockwise polygons
        for ip in range(N-1):
            if not polygon.mask[ip,0]:
                eep = ip + 1
                if eep == N: eep = 0
      
                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
                    icut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])

                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                    ecut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
                    Ncuts = Ncuts + 1
    else:
        # Assuming clockwise polygons
        for ip in range(N-1):
            eep = ip + 1
            if eep == N: eep = 0
      
            if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
                icut.append(ip)
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = yval - polygon[ip,0]
                ipt.append([yval, polygon[ip,1]+dx*dd/dy])

            if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                ecut.append(ip)
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = yval - polygon[ip,0]
                ept.append([yval, polygon[ip,1]+dx*dd/dy])
                Ncuts = Ncuts + 1

    # Looking for repeated
    newicut = icut + []
    newecut = ecut + []
    newipt = ipt + []
    newept = ept + []
    newNcuts = Ncuts
    for ic in range(newNcuts-1):
        for ic2 in range(ic+1,newNcuts):
            if newipt[ic] == newipt[ic2]:
                Ncuts = Ncuts-1
                icut.pop(ic2)
                ecut.pop(ic2)
                ipt.pop(ic2)
                ept.pop(ic2)
                newNcuts = Ncuts + 0

    if ipt is None or ept is None or Ncuts == 0:
        print errormsg
        print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
    else:
        print '  ' + fname + ': found ', Ncuts, ' Ncuts'
        if Ncuts > 1 and keep == 'below':
            # Re-shifting cuts by closest distance.
            xis = []
            xes = []
            for ic in range(Ncuts):
                xp = ipt[ic]
                xis.append(xp[1])
                xp = ept[ic]
                xes.append(xp[1])
            xs = xis + xes
            xs.sort()
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            icut = []
            ecut = []
            ipt = []
            ept = []
            for xv in xs:
                ic = xis.count(xv)
                if ic != 0: 
                    icc = xis.index(xv)
                    if len(icut) > len(ecut): 
                        ecut.append(newicut[icc])
                        ept.append(newipt[icc])
                    else: 
                        icut.append(newicut[icc])
                        ipt.append(newipt[icc])
                else:
                    icc = xes.index(xv)
                    if len(icut) > len(ecut): 
                        ecut.append(newecut[icc])
                        ept.append(newept[icc])
                    else: 
                        icut.append(newecut[icc])
                        ipt.append(newept[icc])

#            # Re-shifting cuts. 1st icut --> last ecut; 1st ecut as 1st icut; 
#            #    2nd icut --> last-1 ecut, ....
#            newicut = icut + []
#            newecut = ecut + []
#            newipt = ipt + []
#            newept = ept + []
#            for ic in range(Ncuts-1):
#                ecut[ic] = newecut[Ncuts-ic-1]
#                ept[ic] = newept[Ncuts-ic-1]
#                icut[ic+1] = newecut[ic]
#                ipt[ic+1] = newept[ic]

#            ecut[Ncuts-1] = newicut[Ncuts-1]
#            ept[Ncuts-1] = newipt[Ncuts-1]

##        print '    yval=', yval, 'cut, ip; ipt ep; ept ________'
##        for ic in range(Ncuts):
##            print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 

    # Length of joining lines
    Nadds = []
    if Ncuts > 1:
        Naddc = (Nadd-Ncuts)/(Ncuts)
        if Naddc < 3:
            print errormsg
            print '  ' + fname + ': too few points for jioning lines !!'
            print '    increase Nadd at least to:', Ncuts*3+Ncuts
            quit(-1)
        for ic in range(Ncuts-1):
            Nadds.append(Naddc)

        Nadds.append(Nadd-Naddc*(Ncuts-1))
    else:
        Nadds.append(Nadd)

    # Total points cut polygon
    Ntotpts = 0
    Ncpts = []
    for ic in range(Ncuts):
        if keep == 'below':
            if ic == 0:
                dpts = icut[ic] + Nadds[ic] + (N - ecut[ic])
            else:
                dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
 
            # Adding end of the polygon in 'left' keeps
            if ic == Ncuts - 1: dpts = dpts + N-ecut[ic]
        else:
            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1

        Ntotpts = Ntotpts + dpts
        Ncpts.append(ecut[ic] - icut[ic])

    cutpolygon = np.ones((Ntotpts+Ncuts,2), dtype=np.float)*gen.fillValue

    iipc = 0
    for ic in range(Ncuts):
        dcpt = Ncpts[ic]
        if keep == 'below':
            if ic == 0:
                cutpolygon[0:icut[ic],:] = polygon[0:icut[ic],:]
                iipc = icut[ic]
            else:
                cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:]
                iipc = iipc + dcpt -1
        else:
            cutpolygon[iipc,:] = ipt[ic]
            cutpolygon[iipc:iipc+dcpt-1,:]=polygon[icut[ic]+1:ecut[ic],:]
            iipc = iipc+dcpt-1

        # cutting line
        cutline = np.zeros((Nadds[ic],2), dtype=np.float)
        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
        cutline[0,:] = ipt[ic]
        for ip in range(1,Nadds[ic]-1):
            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
        cutline[Nadds[ic]-1,:] = ept[ic]
        if keep == 'below':
            if ic == 0: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline
            else: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
            iipc = iipc+Nadds[ic]
            if ic == 0:
                cutpolygon[iipc:iipc+N-ecut[ic]-1,:] = polygon[ecut[ic]+1:N,:]
                iipc = iipc + N-ecut[ic]-1
                cutpolygon[iipc,:] = polygon[0,:]
                iipc = iipc + 1
        else:
            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
            iipc = iipc+Nadds[ic]
        iipc = iipc + 1

    rmpolygon = []
    Npts = cutpolygon.shape[0]
    if keep == 'below':
        for ip in range(Npts):
            if cutpolygon[ip,0] > yval:
                rmpolygon.append([gen.fillValueF, gen.fillValueF])
            else:
                rmpolygon.append(cutpolygon[ip,:])
    else:
        for ip in range(Npts):
            if cutpolygon[ip,0] < yval:
                rmpolygon.append([gen.fillValueF, gen.fillValueF])
            else:
                rmpolygon.append(cutpolygon[ip,:])
    Npts = len(rmpolygon)
    cutpolygon = np.array(rmpolygon)
    cutpolygon = rm_consecpt_polygon(cutpolygon)
    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)

    return Npts, cutpolygon

def cut_xpolygon(polygon, xval, keep='left', Nadd=20):
    """ Function to cut a polygon from a given value of the x-axis
      polygon: polygon to cut
      yval: value to use to cut the polygon
      keep: part to keep from the value ('left', default)
         'left': left of the value
         'right': right of the value
      Nadd: additional points to add to draw the line (20, default)
    """
    fname = 'cut_xpolygon'

    N = polygon.shape[0]
    availkeeps = ['left', 'right']

    if not gen.searchInlist(availkeeps, keep):
        print errormsg
        print '  ' + fname + ": wring keep '" + keep + "' value !!"
        print '    available ones:', availkeeps
        quit(-1)

    ipt = None
    ept = None

    # There might be more than 1 cut ...
    icut = []
    ecut = []
    ipt = []
    ept = []
    Ncuts = 0
    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
      type(gen.mamat.mask[1]):
        # Assuming clockwise polygons
        for ip in range(N-1):
            if not polygon.mask[ip,1]:
                eep = ip + 1
                if eep == N: eep = 0
      
                if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
                    icut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ipt.append([polygon[ip,0]+dy*dd/dx, xval])

                if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
                    ecut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ept.append([polygon[ip,0]+dy*dd/dx, xval])
                    Ncuts = Ncuts + 1
    else:
        # Assuming clockwise polygons
        for ip in range(N-1):
            eep = ip + 1
            if eep == N: eep = 0
      
            if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
                icut.append(ip)
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = xval - polygon[ip,1]
                ipt.append([polygon[ip,0]+dy*dd/dx, xval])

            if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
                ecut.append(ip)
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = xval - polygon[ip,1]
                ept.append([polygon[ip,0]+dy*dd/dx, xval])
                Ncuts = Ncuts + 1

    # Looking for repeated
    newicut = icut + []
    newecut = ecut + []
    newipt = ipt + []
    newept = ept + []
    newNcuts = Ncuts
    for ic in range(newNcuts-1):
        for ic2 in range(ic+1,newNcuts):
            if newipt[ic] == newipt[ic2]:
                Ncuts = Ncuts-1
                icut.pop(ic2)
                ecut.pop(ic2)
                ipt.pop(ic2)
                ept.pop(ic2)
                newNcuts = Ncuts + 0

    if ipt is None or ept is None or Ncuts == 0:
        print errormsg
        print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
    else:
        ##print '  ' + fname + ': found ', Ncuts, ' Ncuts'
        if Ncuts >= 1 and keep == 'left':
            # Re-shifting cuts by closest heigth.
            yis = []
            yes = []
            for ic in range(Ncuts):
                yp = ipt[ic]
                yis.append(yp[0])
                yp = ept[ic]
                yes.append(yp[0])
            ys = yis + yes
            ys.sort()
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            icut = []
            ecut = []
            ipt = []
            ept = []
            for yv in ys:
                ic = yis.count(yv)
                if ic != 0: 
                    icc = yis.index(yv)
                    if len(icut) > len(ecut): 
                        ecut.append(newicut[icc])
                        ept.append(newipt[icc])
                    else: 
                        icut.append(newicut[icc])
                        ipt.append(newipt[icc])
                else:
                    icc = yes.index(yv)
                    if len(icut) > len(ecut): 
                        ecut.append(newecut[icc])
                        ept.append(newept[icc])
                    else: 
                        icut.append(newecut[icc])
                        ipt.append(newept[icc])
        #print '    xval=', xval, 'cut, ip; ipt ep; ept ________'
        #for ic in range(Ncuts):
        #    print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic] 

    # Length of joining lines
    Nadds = []    
    if Ncuts > 1:
        Naddc = (Nadd-Ncuts)/(Ncuts)
        if Naddc < 3:
            print errormsg
            print '  ' + fname + ': too few points for jioning lines !!'
            print '    increase Nadd at least to:', Ncuts*3+Ncuts
            quit(-1)
        for ic in range(Ncuts-1):
            Nadds.append(Naddc)

        Nadds.append(Nadd-Naddc*(Ncuts-1))
    else:
        Nadds.append(Nadd)

   # Total points cut polygon
    Ntotpts = 0
    Ncpts = []
    for ic in range(Ncuts):
        if keep == 'left':
            if ic == 0: 
                dpts = icut[ic] + Nadds[ic] + (N - ecut[ic])
            else:
                dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1

            # Adding end of the polygon in 'left' keeps
            if ic == Ncuts - 1: dpts = dpts + N-ecut[ic]
        else:
            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1

        Ntotpts = Ntotpts + dpts
        Ncpts.append(ecut[ic] - icut[ic])

    cutpolygon = []
    iipc = 0
    for ic in range(Ncuts):
        dcpt = Ncpts[ic]
        cutpolygon.append(ipt[ic])
        if keep == 'left':
            if ic == 0:
                add_secpolygon_list(cutpolygon,icut[ic]+1,N,polygon)
                add_secpolygon_list(cutpolygon,0,ecut[ic],polygon)
                iipc = icut[ic]
            else:
                add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
        else:
            add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
            iipc = iipc+dcpt-1
        # cutting line
        cutline = np.zeros((Nadds[ic],2), dtype=np.float)
        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
        cutline[0,:] = ipt[ic]
        for ip in range(1,Nadds[ic]-1):
            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
        cutline[Nadds[ic]-1,:] = ept[ic]
        if keep == 'left':
            for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:])
            iipc = iipc+Nadds[ic]
            if ic == 0:
                add_secpolygon_list(cutpolygon,ecut[ic],N,polygon)
                cutpolygon.append(polygon[0,:])
                iipc = iipc + 1
        else:
            for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:])
            iipc = iipc+Nadds[ic]
        cutpolygon.append([gen.fillValueF, gen.fillValueF])
        iipc = iipc + 1

    cutpolygon = np.array(cutpolygon)
    rmpolygon = []
    Npts = cutpolygon.shape[0]
    if keep == 'left':
        for ip in range(Npts):
            if cutpolygon[ip,1] > xval:
                rmpolygon.append([gen.fillValueF, gen.fillValueF])
            else:
                rmpolygon.append(cutpolygon[ip,:])
    else:
        for ip in range(Npts):
            if cutpolygon[ip,1] < xval:
                rmpolygon.append([gen.fillValueF, gen.fillValueF])
            else:
                rmpolygon.append(cutpolygon[ip,:])

    rmpolygon = np.array(rmpolygon)
    cutpolygon = rm_consecpt_polygon(rmpolygon)
    Npts = cutpolygon.shape[0]

    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)

    return Npts, cutpolygon

def cut_between_ypolygon(polygon, yval1, yval2, Nadd=20):
    """ Function to cut a polygon between 2 given value of the y-axis
      polygon: polygon to cut
      yval1: first value to use to cut the polygon
      yval2: first value to use to cut the polygon
      Nadd: additional points to add to draw the line (20, default)
    """
    fname = 'cut_betwen_ypolygon'

    N = polygon.shape[0]

    if yval1 > yval2:
        print errormsg
        print '  ' + fname + ': wrong between cut values !!'
        print '     it is expected yval1 < yval2'
        print '     values provided yval1: (', yval1, ')> yval2 (', yval2, ')'
        quit(-1)

    yvals = [yval1, yval2]

    ipt = None
    ept = None

    cuts = {}
    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
      type(gen.mamat.mask[1]):
        for ic in range(2):
            yval = yvals[ic]
            # There might be more than 1 cut ...
            icut = []
            ecut = []
            ipt = []
            ept = []
            Ncuts = 0
            # Assuming clockwise polygons
            for ip in range(N-1):
                if not polygon.mask[ip,0]:
                    eep = ip + 1
                    if eep == N: eep = 0
      
                    if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
                        icut.append(ip)
                        dx = polygon[eep,1] - polygon[ip,1]
                        dy = polygon[eep,0] - polygon[ip,0]
                        dd = yval - polygon[ip,0]
                        ipt.append([yval, polygon[ip,1]+dx*dd/dy])

                    if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                        ecut.append(ip)
                        dx = polygon[eep,1] - polygon[ip,1]
                        dy = polygon[eep,0] - polygon[ip,0]
                        dd = yval - polygon[ip,0]
                        ept.append([yval, polygon[ip,1]+dx*dd/dy])
                        Ncuts = Ncuts + 1

            # Looking for repeated
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            newNcuts = Ncuts
            for icp in range(newNcuts-1):
                for ic2 in range(icp+1,newNcuts):
                    if newipt[icp] == newipt[ic2]:
                        Ncuts = Ncuts-1
                        icut.pop(ic2)
                        ecut.pop(ic2)
                        ipt.pop(ic2)
                        ept.pop(ic2)
                        newNcuts = Ncuts + 0

            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
    else:
        for ic in range(2):
            yval = yvals[ic]
            # There might be more than 1 cut ...
            icut = []
            ecut = []
            ipt = []
            ept = []
            Ncuts = 0
            # Assuming clockwise polygons
            for ip in range(N-1):
                eep = ip + 1
                if eep == N: eep = 0
      
                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
                    icut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])

                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
                    ecut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
                    Ncuts = Ncuts + 1
            # Looking for repeated
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            newNcuts = Ncuts
            for icp in range(newNcuts-1):
                for ic2 in range(icp+1,newNcuts):
                    if newipt[icp] == newipt[ic2]:
                        Ncuts = Ncuts-1
                        icut.pop(ic2)
                        ecut.pop(ic2)
                        ipt.pop(ic2)
                        ept.pop(ic2)
                        newNcuts = Ncuts + 0

            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]

    Naddlines = {}
    for icc in range(2):
        cutv = cuts[icc]
        Ncuts = cutv[4]
        # Length of joining lines
        Nadds = []
        if Ncuts > 1:
            Naddc = (Nadd-Ncuts)/(Ncuts)
            if Naddc < 3:
                print errormsg
                print '  ' + fname + ': too few points for jioning lines !!'
                print '    increase Nadd at least to:', Ncuts*3+Ncuts
                quit(-1)
            for ic in range(Ncuts-1):
                Nadds.append(Naddc)

            Nadds.append(Nadd-Naddc*(Ncuts-1))
        else:
            Nadds.append(Nadd)

        # Total points cut polygon
        Ntotpts = 0
        Ncpts = []
        for ic in range(Ncuts):
            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1

            Ntotpts = Ntotpts + dpts
            Ncpts.append(ecut[ic] - icut[ic])

        Naddlines[icc] = [Nadds, Ntotpts, Ncpts]

    cutv1 = cuts[0]
    addv1 = Naddlines[0]
    Nadds1 = addv1[0]
    Ncuts1 = cutv1[4]

    cutv2 = cuts[1]
    addv2 = Naddlines[1]
    Nadds2 = addv2[0]
    Ncuts2 = cutv2[4]

    if Ncuts1 != Ncuts2: 
        print errormsg
        print '  ' + fname + ": different number of cuts !!"
        print '    yval1:', yval1, 'Ncuts=', Ncuts1
        print '    yval2:', yval2, 'Ncuts=', Ncuts2
        print '      I am not prepare to deal with it'
        quit(-1)
    #else:
    #    print '  ' + fname + ' _______'
    #    print '    yval1:', yval1, 'Ncuts=', Ncuts1
    #    print '    yval2:', yval2, 'Ncuts=', Ncuts2

    icut1 = cutv1[0]
    ecut1 = cutv1[1]
    ipt1 = cutv1[2]
    ept1 = cutv1[3]
    icut2 = cutv2[0]
    ecut2 = cutv2[1]
    ipt2 = cutv2[2]
    ept2 = cutv2[3]

    # Looking for pairs of cuts. Grouping for smallest x distance between initial    
    #   points of each cut
    cutpolygons = []
    for iic1 in range(Ncuts1):
        iip = 0
        cutpolygon = []
        ic1 = icut1[iic1]
        ec1 = ecut1[iic1]
        ip1 = ipt1[iic1]
        ep1 = ept1[iic1]

        ipx2s = []
        for ip in range(Ncuts2): 
            ip2 = ipt2[ip]
            ipx2s.append(ip2[1])
        dxps = ipx2s - ip1[1]
        dxps = np.where(dxps < 0., gen.fillValueF, dxps)
        ndxps = np.min(dxps)
        iic12 = gen.index_vec(dxps,ndxps)

        ic2 = icut2[iic12]
        ec2 = ecut2[iic12]
        ip2 = ipt2[iic12]
        ep2 = ept2[iic12]

        #print 'Lluis iic1', iic1, 'ic1', ic1, 'ec1', ec1, 'ipt1', ip1, 'ept1', ep1, 'Nadds1', Nadds1
        #print '    iic12', iic12, 'ic2', ic2, 'ec2', ec2, 'ipt2', ip2, 'ept2', ep2, 'Nadds2', Nadds2

        cutpolygon.append(ip1)
        for ip in range(ic1+1,ic2-1):
            cutpolygon.append(polygon[ip,:])
        iip = ic2-ic1
        # cutting line 1
        Nadd2 = Nadds1[iic1]
        cutlines = np.zeros((Nadd2,2), dtype=np.float)
        dx = (ep2[1] - ip2[1])/(Nadd2-2)
        dy = (ep2[0] - ip2[0])/(Nadd2-2)
        cutlines[0,:] = ip2
        for ip in range(1,Nadd2-1):
            cutlines[ip,:] = ip2 + np.array([dy*ip,dx*ip])
        cutlines[Nadd2-1,:] = ep2
        for ip in range(Nadd2): cutpolygon.append(cutlines[ip,:])
        iip = iip + Nadd2

        for ip in range(ec2,ec1):
            cutpolygon.append(polygon[ip,:])
        iip = iip + ec1-ec2
        # cutting line 2
        Nadd2 = Nadds2[iic12]
        cutlines = np.zeros((Nadd2,2), dtype=np.float)
        dx = (ep1[1] - ip1[1])/(Nadd2-2)
        dy = (ep1[0] - ip1[0])/(Nadd2-2)
        cutlines[0,:] = ip1
        for ip in range(1,Nadd2-1):
            cutlines[ip,:] = ip1 + np.array([dy*ip,dx*ip])
        cutlines[Nadd2-1,:] = ep1
        for ip in range(Nadd2-1,0,-1):
            cutpolygon.append(cutlines[ip,:])

        cutpolygon.append(ip1)

        cutpolygon.append([gen.fillValueF,gen.fillValueF])
        if len(cutpolygons) == 0: cutpolygons = cutpolygon
        else: cutpolygons = cutpolygons + cutpolygon 

    cutpolygons = np.array(cutpolygons)
    cutpolygons = rm_consecpt_polygon(cutpolygons)
    cutpolygons = ma.masked_equal(cutpolygons, gen.fillValueF)

    Npts = cutpolygons.shape[0]

    return Npts, cutpolygons

def cut_between_xpolygon(polygon, xval1, xval2, Nadd=20):
    """ Function to cut a polygon between 2 given value of the x-axis
      polygon: polygon to cut
      xval1: first value to use to cut the polygon
      xval2: first value to use to cut the polygon
      Nadd: additional points to add to draw the line (20, default)
    """
    fname = 'cut_betwen_xpolygon'

    N = polygon.shape[0]

    if xval1 > xval2:
        print errormsg
        print '  ' + fname + ': wrong between cut values !!'
        print '     it is expected xval1 < xval2'
        print '     values provided xval1: (', xval1, ')> xval2 (', xval2, ')'
        quit(-1)

    xvals = [xval1, xval2]

    ipt = None
    ept = None

    cuts = {}
    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
      type(gen.mamat.mask[1]):
        for ic in range(2):
            xval = xvals[ic]
            # There might be more than 1 cut ...
            icut = []
            ecut = []
            ipt = []
            ept = []
            Ncuts = 0
            # Assuming clockwise polygons
            for ip in range(N-1):
                if not polygon.mask[ip,0]:
                    eep = ip + 1
                    if eep == N: eep = 0
      
                    if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
                        icut.append(ip)
                        dx = polygon[eep,1] - polygon[ip,1]
                        dy = polygon[eep,0] - polygon[ip,0]
                        dd = xval - polygon[ip,1]
                        ipt.append([polygon[ip,0]+dy*dd/dx, xval])

                    if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
                        ecut.append(ip)
                        dx = polygon[eep,1] - polygon[ip,1]
                        dy = polygon[eep,0] - polygon[ip,0]
                        dd = xval - polygon[ip,1]
                        ept.append([polygon[ip,0]+dy*dd/dx, xval])
                        Ncuts = Ncuts + 1

            # Looking for repeated
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            newNcuts = Ncuts
            for icp in range(newNcuts-1):
                for ic2 in range(icp+1,newNcuts):
                    if newipt[icp] == newipt[ic2]:
                        Ncuts = Ncuts-1
                        icut.pop(ic2)
                        ecut.pop(ic2)
                        ipt.pop(ic2)
                        ept.pop(ic2)
                        newNcuts = Ncuts + 0

            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
    else:
        for ic in range(2):
            xval = xvals[ic]
            # There might be more than 1 cut ...
            icut = []
            ecut = []
            ipt = []
            ept = []
            Ncuts = 0
            # Assuming clockwise polygons
            for ip in range(N-1):
                eep = ip + 1
                if eep == N: eep = 0
      
                if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
                    icut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ipt.append([polygon[ip,0]+dy*dd/dx, xval])

                if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
                    ecut.append(ip)
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ept.append([polygon[ip,0]+dy*dd/dx, xval])
                    Ncuts = Ncuts + 1
            # Looking for repeated
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            newNcuts = Ncuts
            for icp in range(newNcuts-1):
                for ic2 in range(icp+1,newNcuts):
                    if newipt[icp] == newipt[ic2]:
                        Ncuts = Ncuts-1
                        icut.pop(ic2)
                        ecut.pop(ic2)
                        ipt.pop(ic2)
                        ept.pop(ic2)
                        newNcuts = Ncuts + 0

            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]

    for iic in range(1):
        cutvs = cuts[iic]
        icut = cutvs[0]
        ecut = cutvs[1]
        ipt = cutvs[2]
        ept = cutvs[3]
        Ncuts = cutvs[4]
        if Ncuts > 0:
            # Re-shifting cuts by closest heigth.
            yis = []
            yes = []
            for ic in range(Ncuts):
                yp = ipt[ic]
                yis.append(yp[0])
                yp = ept[ic]
                yes.append(yp[0])
            ys = yis + yes
            ys.sort()
            newicut = icut + []
            newecut = ecut + []
            newipt = ipt + []
            newept = ept + []
            icut = []
            ecut = []
            ipt = []
            ept = []
            for yv in ys:
                ic = yis.count(yv)
                if ic != 0: 
                    icc = yis.index(yv)
                    if len(icut) > len(ecut): 
                        ecut.append(newicut[icc])
                        ept.append(newipt[icc])
                    else: 
                        icut.append(newicut[icc])
                        ipt.append(newipt[icc])
                else:
                    icc = yes.index(yv)
                    if len(icut) > len(ecut): 
                        ecut.append(newecut[icc])
                        ept.append(newept[icc])
                    else: 
                        icut.append(newecut[icc])
                        ipt.append(newept[icc])

        cuts[iic] = [icut, ecut, ipt, ept, Ncuts]

    Naddlines = {}
    for icc in range(2):
        cutv = cuts[icc]
        Ncuts = cutv[4]
        if Ncuts == 0:
            print errormsg
            print '  ' + fname + ": no cuts for xval=", xvals[icc], '!!'
            quit(-1)
        #print '  icc:', icc, 'ic ec ipt ept _______'
        #for ic in range(Ncuts):
        #    print ic, ':', cutv[0][ic], cutv[1][ic], cutv[2][ic], cutv[3][ic]

        # Length of joining lines
        Nadds = []
        if Ncuts > 1:
            Naddc = (Nadd-Ncuts)/(Ncuts)
            if Naddc < 3:
                print errormsg
                print '  ' + fname + ': too few points for jioning lines !!'
                print '    increase Nadd at least to:', Ncuts*3+Ncuts
                quit(-1)
            for ic in range(Ncuts-1):
                Nadds.append(Naddc)

            Nadds.append(Nadd-Naddc*(Ncuts-1))
        else:
            Nadds.append(Nadd)

        Naddlines[icc] = Nadds

    # sides
    sides = {}
    for iic in range(2):
        cutvs = cuts[iic]
        icut = cutvs[0]
        ecut = cutvs[1]
        ipt = cutvs[2]
        ept = cutvs[3]
        Ncuts = cutvs[4]
        Nadds = Naddlines[iic]
        cutpolygon = []
        # left side
        if iic == 0:
            for ic in range(Ncuts-1):
                cutpolygon.append(ipt[ic])
                dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
                dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
                for ip in range(1,Nadds[ic]-1):
                    cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
                cutpolygon.append(ept[ic])
                for ip in range(ecut[ic]+1,icut[ic+1]): cutpolygon.append(polygon[ip,:])

            ic = Ncuts-1
            cutpolygon.append(ipt[ic])
            dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
            dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
            for ip in range(1,Nadds[ic]-1):
                cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
        # right side
        else:
            for ic in range(Ncuts-1):
                cutpolygon.append(ipt[ic])

                # line
                dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
                dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
                for ip in range(1,Nadds[ic]-1):
                    cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
                cutpolygon.append(ept[ic])
                for ip in range(ecut[ic],icut[ic+1]): cutpolygon.append(polygon[ip,:])

            ic = Ncuts-1
            cutpolygon.append(ipt[ic])
            # line
            dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
            dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
            for ip in range(1,Nadds[ic]-1):
                cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
            cutpolygon.append(ept[ic])
        sides[iic] = cutpolygon
        
    # joining sides by e1[Ncuts1-1] --> i2[0]; e2[Ncuts2-1] --> i1[0]
    cutv1 = cuts[0]
    Ncuts1 = cutv1[4]
    ec1 = cutv1[1][np.max([0,Ncuts1-1])]
    ic1 = cutv1[0][0]
    ept1 = cutv1[3][np.max([0,Ncuts1-1])]
    ipt1 = cutv1[2][0]

    cutv2 = cuts[1]
    Ncuts2 = cutv2[4]
    ec2 = cutv2[1][np.max([0,Ncuts2-1])]
    ic2 = cutv2[0][0]
    ept2 = cutv2[3][np.max([0,Ncuts2-1])]
    ipt2 = cutv2[2][0]

    finalcutpolygon = sides[0]
    for ip in range(ec1+1,ic2): finalcutpolygon.append(polygon[ip,:])
    finalcutpolygon = finalcutpolygon + sides[1]
    for ip in range(ec2+1,ic1): finalcutpolygon.append(polygon[ip,:])
    finalcutpolygon.append(ipt1)

    finalcutpolygon = np.array(finalcutpolygon)

    finalcutpolygon = rm_consecpt_polygon(finalcutpolygon)
    finalcutpolygon = ma.masked_equal(finalcutpolygon, gen.fillValueF)

    Npts = finalcutpolygon.shape[0]

    return Npts, finalcutpolygon

def pile_polygons(polyns, polygons):
    """ Function to pile polygons one over the following one
      polyns: ordered list of polygons. First over all. last below all
      polygons: dictionary with the polygons
    >>> pns = ['sqra', 'sqrb']
    >>> polya = np.array([[-0.5, -0.75], [0.5, -0.75], [0.5, 0.75], [-0.5, 0.75]])
    >>> polyb = np.array([[-0.75, -0.5], [0.75, -0.5], [0.75, 0.5], [-0.75, 0.5]])
    >>> plgs = {'sqra': polya, 'sqrb': polyb}
    >>> pile_polygons(pns, plgs)
    #  sqrb : 
    [[-0.75 -0.5]
     [-0.5 -0.5]
     [-- --]
     [0.5 -0.5]
     [0.75 -0.5]
     [0.75 0.5]
     [0.5 0.5]
     [-- --]
     [-0.5 0.5]
     [-0.75 0.5]
     [-0.75 -0.5]]
    #  sqra : 
     [[-0.5  -0.75]
      [ 0.5  -0.75]
      [ 0.5   0.75]
      [-0.5   0.75]
      [-0.5  -0.75]]
    """
    fname = 'pile_polygons'
    pilepolygons = dict(polygons)
    Npolys = len(polyns)

    for ipolyp in range(Npolys-2,-1,-1):
        polyn = polyns[ipolyp]
        poly = pilepolygons[polyn]
        Npts = poly.shape[0]
        for ipolyi in range(ipolyp+1,Npolys,1):
            ipolyn = polyns[ipolyi]
            #print '  Lluis ' + polyn + ' above ' + ipolyn
            ipoly = pilepolygons[ipolyn]
            iNpts = ipoly.shape[0]
            newipoly = []

            Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys(     \
              nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,      \
              polyb=poly)
            # We're in C-mode !
            inti = inti-1
            intp = intp-1

            # Re-constructing below polygon looking for respective crossings
            linti = list(inti)
            for ip in range(iNpts):
                iip1 = ip+1
                if ip == iNpts-1: iip1 = 0
                #print ip, ipoly[ip,:], ':', ipoly[iip1,:]
                Nc = linti.count(ip)
                ldists = []
                ddists = {}
                if Nc > 0:
                    iic = gen.multi_index_vec(inti,ip)
                    mindist = 1000000.
                    # Sorting from distance respect the vertex ip
                    for ic in range(Nc):
                        ddists[iic[ic]] = dist_points(ipoly[ip,:], pts[iic[ic],:])
                        ldists.append(ddists[iic[ic]])
                        #print '  ', ic, ';', iic[ic], '=', pts[iic[ic],:], ddists[iic[ic]]
                    ldists.sort()
                    #print '    ldists', ldists
                    newipoly.append(ipoly[ip,:])
                    for ic in range(Nc):
                        iic = gen.dictionary_key(ddists,ldists[ic])
                        newipoly.append(pts[iic,:])
                        #print '  ', ic, '|', iic, ';', pts[iic,:]
                        if ic < Nc-1: newipoly.append([gen.fillValueF, gen.fillValueF])
                    newipoly.append(ipoly[iip1,:])

                else:
                    newipoly.append(ipoly[ip,:])
    
            newipoly = np.array(newipoly)
            pilepolygons[polyns[Npolys-1]] = rm_consecpt_polygon(newipoly)
    
        for polyn in polyns:
            poly = pilepolygons[polyn]
            poly = ma.masked_equal(poly, gen.fillValueF)
            pilepolygons[polyn] = poly

    return pilepolygons

####### ###### ##### #### ### ## #
# Shapes/objects

def surface_sphere(radii,Npts):
    """ Function to provide an sphere as matrix of x,y,z coordinates
      radii: radii of the sphere
      Npts: number of points to discretisize longitues (half for latitudes)
    """
    fname = 'surface_sphere'

    sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float)
    spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float)
    for ia in range(Npts):
        alpha = ia*2*np.pi/(Npts-1)
        for ib in range(Npts/2):
            beta = ib*np.pi/(2.*(Npts/2-1))
            sphereup[:,ib,ia] = position_sphere(radii, alpha, beta)
        for ib in range(Npts/2):
            beta = -ib*np.pi/(2.*(Npts/2-1))
            spheredown[:,ib,ia] = position_sphere(radii, alpha, beta)

    return sphereup, spheredown

def ellipse_polar(c, a, b, Nang=100):
    """ Function to determine an ellipse from its center and polar coordinates
        FROM: https://en.wikipedia.org/wiki/Ellipse
      c= coordinates of the center
      a= distance major axis
      b= distance minor axis
      Nang= number of angles to use
    """
    fname = 'ellipse_polar'

    if np.mod(Nang,2) == 0: Nang=Nang+1
  
    dtheta = 2*np.pi/(Nang-1)

    ellipse = np.zeros((Nang,2), dtype=np.float)
    for ia in range(Nang):
        theta = dtheta*ia
        rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 )
        x = rad*np.cos(theta)
        y = rad*np.sin(theta)
        ellipse[ia,:] = [y+c[0],x+c[1]]

    return ellipse

def hyperbola_polar(a, b, Nang=100):
    """ Fcuntion to determine an hyperbola in polar coordinates
        FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates
          x^2/a^2 - y^2/b^2 = 1
      a= x-parameter
      y= y-parameter
      Nang= number of angles to use
      DOES NOT WORK!!!!
    """
    fname = 'hyperbola_polar'

    dtheta = 2.*np.pi/(Nang-1)

    # Positive branch
    hyperbola_p = np.zeros((Nang,2), dtype=np.float)
    for ia in range(Nang):
        theta = dtheta*ia
        x = a*np.cosh(theta)
        y = b*np.sinh(theta)
        hyperbola_p[ia,:] = [y,x]

    # Negative branch
    hyperbola_n = np.zeros((Nang,2), dtype=np.float)
    for ia in range(Nang):
        theta = dtheta*ia
        x = -a*np.cosh(theta)
        y = b*np.sinh(theta)
        hyperbola_n[ia,:] = [y,x]

    return hyperbola_p, hyperbola_n

def circ_sec(ptA, ptB, radii, arc='short', pos='left', Nang=100):
    """ Function union of point A and B by a section of a circle
      ptA= coordinates od the point A [yA, xA]
      ptB= coordinates od the point B [yB, xB]
      radii= radi of the circle to use to unite the points
      arc= which arc to be used ('short', default)
        'short': shortest angle between points
        'long': largest angle between points
      pos= orientation of the arc following clockwise union of points ('left', default)
        'left': to the left of union
        'right': to the right of union
      Nang= amount of angles to use
    """
    fname = 'circ_sec'
    availarc = ['short', 'long']
    availpos = ['left', 'right']

    distAB = dist_points(ptA,ptB)

    if distAB > radii:
        print errormsg
        print '  ' + fname + ': radii=', radii, " too small for the distance " +     \
          "between points !!"
        print '    distance between points:', distAB
        quit(-1)

    # Coordinate increments
    dAB = np.abs(ptA-ptB)

    # angle of the circular section joining points
    alpha = 2.*np.arcsin((distAB/2.)/radii)

    # center along coincident bisection of the union
    xcc = -radii
    ycc = 0.

    # Getting the arc of the circle at the x-axis
    if arc == 'short':
        dalpha = alpha/(Nang-1)
    elif arc == 'long':
        dalpha = (2.*np.pi - alpha)/(Nang-1)
    else:
        print errormsg
        print '  ' + fname + ": arc '" + arc + "' not ready !!" 
        print '    available ones:', availarc
        quit(-1)
    if pos == 'left': sign=-1.
    elif pos == 'right': sign=1.
    else:
        print errormsg
        print '  ' + fname + ": position '" + pos + "' not ready !!" 
        print '     available ones:', availpos
        quit(-1)

    circ_sec = np.zeros((Nang,2), dtype=np.float)
    for ia in range(Nang):
        alpha = sign*dalpha*ia
        x = radii*np.cos(alpha)
        y = radii*np.sin(alpha)

        circ_sec[ia,:] = [y+ycc,x+xcc]

    # Angle of the points
    theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])

    # rotating angle of the circ
    if pos == 'left': 
        rotangle = theta + np.pi/2. - alpha/2.
    elif pos == 'right':
        rotangle = theta + 3.*np.pi/2. - alpha/2.
    else:
        print errormsg
        print '  ' + fname + ": position '" + pos + "' not ready !!" 
        print '     available ones:', availpos
        quit(-1)

    #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi
  
    # rotating the arc along the x-axis
    rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle)

    # Moving arc to the ptA
    circ_sec = rotcirc_sec + ptA

    return circ_sec

def p_square(face, N=5):
    """ Function to get a polygon square
      face: length of the face of the square
      N: number of points of the polygon
    """
    fname = 'p_square'

    square = np.zeros((N,2), dtype=np.float)

    f2 = face/2.
    N4 = N/4
    df = face/(N4)
    # SW-NW
    for ip in range(N4):
        square[ip,:] = [-f2+ip*df,-f2]
    # NW-NE
    for ip in range(N4):
        square[ip+N4,:] = [f2,-f2+ip*df]
    # NE-SE
    for ip in range(N4):
        square[ip+2*N4,:] = [f2-ip*df,f2]
    N42 = N-3*N4-1
    df = face/(N42)
    # SE-SW
    for ip in range(N42):
        square[ip+3*N4,:] = [-f2,f2-ip*df]
    square[N-1,:] = [-f2,-f2]

    return square


def p_prism(base, height, N=5):
    """ Function to get a polygon prism
      base: length of the base of the prism
      height: length of the height of the prism
      N: number of points of the polygon
    """
    fname = 'p_prism'

    prism = np.zeros((N,2), dtype=np.float)

    b2 = base/2.
    h2 = height/2.
    N4 = N/4
    dh = height/(N4)
    db = base/(N4)

    # SW-NW
    for ip in range(N4):
        prism[ip,:] = [-h2+ip*dh,-b2]
    # NW-NE
    for ip in range(N4):
        prism[ip+N4,:] = [h2,-b2+ip*db]
    # NE-SE
    for ip in range(N4):
        prism[ip+2*N4,:] = [h2-ip*dh,b2]
    N42 = N-3*N4-1
    db = base/(N42)
    # SE-SW
    for ip in range(N42):
        prism[ip+3*N4,:] = [-h2,b2-ip*db]
    prism[N-1,:] = [-h2,-b2]

    return prism

def p_circle(radii, N=50):
    """ Function to get a polygon of a circle
      radii: length of the radii of the circle
      N: number of points of the polygon
    """
    fname = 'p_circle'

    circle = np.zeros((N,2), dtype=np.float)

    dangle = 2.*np.pi/(N-1)

    for ia in range(N):
        circle[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]

    circle[N-1,:] = [0., radii]

    return circle

def p_triangle(p1, p2, p3, N=4):
    """ Function to provide the polygon of a triangle from its 3 vertices
      p1: vertex 1 [y,x]
      p2: vertex 2 [y,x]
      p3: vertex 3 [y,x]
      N: number of vertices of the triangle
    """
    fname = 'p_triangle'

    triangle = np.zeros((N,2), dtype=np.float)

    N3 = N / 3
    # 1-2
    dx = (p2[1]-p1[1])/N3
    dy = (p2[0]-p1[0])/N3
    for ip in range(N3):
        triangle[ip,:] = [p1[0]+ip*dy,p1[1]+ip*dx]
    # 2-3
    dx = (p3[1]-p2[1])/N3
    dy = (p3[0]-p2[0])/N3
    for ip in range(N3):
        triangle[ip+N3,:] = [p2[0]+ip*dy,p2[1]+ip*dx]
    # 3-1
    N32 = N - 2*N/3
    dx = (p1[1]-p3[1])/N32
    dy = (p1[0]-p3[0])/N32
    for ip in range(N32):
        triangle[ip+2*N3,:] = [p3[0]+ip*dy,p3[1]+ip*dx]

    triangle[N-1,:] = p1

    return triangle

def p_spiral(loops, eradii, N=1000):
    """ Function to provide a polygon of an Archimedean spiral
        FROM: https://en.wikipedia.org/wiki/Spiral
      loops: number of loops of the spiral
      eradii: length of the radii of the final spiral
      N: number of points of the polygon
    """
    fname = 'p_spiral'

    spiral = np.zeros((N,2), dtype=np.float)

    dangle = 2.*np.pi*loops/(N-1)
    dr = eradii*1./(N-1)

    for ia in range(N):
        radii = dr*ia
        spiral[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]

    return spiral

def p_reg_polygon(Nv, lf, N=50):
    """ Function to provide a regular polygon of Nv vertices
      Nv: number of vertices
      lf: length of the face
      N: number of points
    """
    fname = 'p_reg_polygon'

    reg_polygon = np.zeros((N,2), dtype=np.float)

    # Number of points per vertex
    Np = N/Nv
    # Angle incremental between vertices
    da = 2.*np.pi/Nv
    # Radii of the circle according to lf
    radii = lf*Nv/(2*np.pi)

    iip = 0
    for iv in range(Nv-1):
        # Characteristics between vertices iv and iv+1
        av1 = da*iv
        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
        av2 = da*(iv+1)
        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
        dx = (v2[1]-v1[1])/Np
        dy = (v2[0]-v1[0])/Np
        for ip in range(Np):
            reg_polygon[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]

    # Characteristics between vertices Nv and 1

    # Number of points per vertex
    Np2 = N - Np*(Nv-1)

    av1 = da*Nv
    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
    av2 = 0.
    v2 = [radii*np.sin(av2), radii*np.cos(av2)]
    dx = (v2[1]-v1[1])/Np2
    dy = (v2[0]-v1[0])/Np2
    for ip in range(Np2):
        reg_polygon[ip+(Nv-1)*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]

    return reg_polygon

def p_reg_star(Nv, lf, freq, vs=0, N=50):
    """ Function to provide a regular star of Nv vertices
      Nv: number of vertices
      lf: length of the face of the regular polygon
      freq: frequency of union of vertices ('0', for just centered to zero arms)
      vs: vertex from which start (0 being first [0,lf])
      N: number of points
    """
    fname = 'p_reg_star'

    reg_star = np.zeros((N,2), dtype=np.float)

    # Number of arms of the star
    if freq != 0 and np.mod(Nv,freq) == 0: 
        Na = Nv/freq + 1
    else:
        Na = Nv

    # Number of points per arm
    Np = N/Na
    # Angle incremental between vertices
    da = 2.*np.pi/Nv
    # Radii of the circle according to lf
    radii = lf*Nv/(2*np.pi)

    iip = 0
    av1 = vs*da
    for iv in range(Na-1):
        # Characteristics between vertices iv and iv+1
        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
        if freq != 0:
            av2 = av1 + da*freq
            v2 = [radii*np.sin(av2), radii*np.cos(av2)]
        else:
            v2 = [0., 0.]
            av2 = av1 + da
        dx = (v2[1]-v1[1])/(Np-1)
        dy = (v2[0]-v1[0])/(Np-1)
        for ip in range(Np):
            reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
        if av2 > 2.*np.pi: av1 = av2 - 2.*np.pi
        else: av1 = av2 + 0.

    iv = Na-1
    # Characteristics between vertices Na and 1
    Np2 = N-Np*iv
    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
    if freq != 0:
        av2 = vs*da
        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
    else:
        v2 = [0., 0.]
    dx = (v2[1]-v1[1])/(Np2-1)
    dy = (v2[0]-v1[0])/(Np2-1)
    for ip in range(Np2):
        reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]

    return reg_star

def p_sinusoide(length=10., amp=5., lamb=3., ival=0., func='sin', N=100):
    """ Function to get coordinates of a sinusoidal curve
      length: length of the line (default 10.)
      amp: amplitude of the peaks (default 5.)
      lamb: wave longitude (defalult 3.)
      ival: initial angle (default 0. in degree)
      func: function to use: (default sinus)
        'sin': sinus
        'cos': cosinus
      N: number of points (default 100)
    """
    fname = 'p_sinusoide'
    availfunc = ['sin', 'cos']

    dx = length/(N-1)
    ia = ival*np.pi/180.
    da = 2*np.pi*dx/lamb

    sinusoide = np.zeros((N,2), dtype=np.float)
    if func == 'sin':
        for ix in range(N):
            sinusoide[ix,:] = [amp*np.sin(ia+da*ix),dx*ix]
    elif func == 'cos':
        for ix in range(N):
            sinusoide[ix,:] = [amp*np.cos(ia+da*ix),dx*ix]
    else:
        print errormsg
        print '  ' + fname + ": function '" + func + "' not ready !!"
        print '    available ones:', availfunc
        quit(-1)

    sinusoidesecs = ['sinusoide']
    sinusoidedic = {'sinusoide': [sinusoide, '-', '#000000', 1.]}

    return sinusoide, sinusoidesecs, sinusoidedic

def p_doubleArrow(length=5., angle=45., width=1., alength=0.10, N=50):
    """ Function to provide an arrow with double lines
      length: length of the arrow (5. default)
      angle: angle of the head of the arrow (45., default)
      width: separation between the two lines (2., default)
      alength: length of the head (as percentage in excess of width, 0.1 default)
      N: number of points (50, default)
    """
    function = 'p_doubleArrow'

    doubleArrow = np.zeros((50,2), dtype=np.float)
    N4 = int((N-3)/4)

    doublearrowdic = {}
    ddy = width*np.tan(angle*np.pi/180.)/2.
    # Arms
    dx = (length-ddy)/(N4-1)
    for ix in range(N4):
        doubleArrow[ix,:] = [dx*ix,-width/2.]
    doublearrowdic['leftarm'] = [doubleArrow[0:N4,:], '-', '#000000', 2.]
    doubleArrow[N4,:] = [gen.fillValueF,gen.fillValueF]
    for ix in range(N4):
        doubleArrow[N4+1+ix,:] = [dx*ix,width/2.]
    doublearrowdic['rightarm'] = [doubleArrow[N4+1:2*N4+1,:], '-', '#000000', 2.]
    doubleArrow[2*N4+1,:] = [gen.fillValueF,gen.fillValueF]

    # Head
    N42 = int((N-2 - 2*N4)/2)
    dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N42-1)
    dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N42-1)
    for ix in range(N42):
        doubleArrow[2*N4+2+ix,:] = [length-dy*ix,-dx*ix]
    doublearrowdic['lefthead'] = [doubleArrow[2*N4:2*N4+N42,:], '-', '#000000', 2.]
    doubleArrow[2*N4+2+N42,:] = [gen.fillValueF,gen.fillValueF]

    N43 = N-3 - 2*N4 - N42 + 1
    dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N43-1)
    dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N43-1)
    for ix in range(N43):
        doubleArrow[2*N4+N42+2+ix,:] = [length-dy*ix,dx*ix]
    doublearrowdic['rightthead'] = [doubleArrow[2*N4+N42+2:51,:], '-', '#000000', 2.]

    doubleArrow = ma.masked_equal(doubleArrow, gen.fillValueF)
    doublearrowsecs = ['leftarm', 'rightarm', 'lefthead', 'righthead']

    return doubleArrow, doublearrowsecs, doublearrowdic 

def p_angle_triangle(pi=np.array([0.,0.]), angle1=60., length1=1., angle2=60., N=100):
    """ Function to draw a triangle by an initial point and two consecutive angles 
        and the first length of face. The third angle and 2 and 3rd face will be 
        computed accordingly the provided values:
           length1 / sin(angle1) = length2 / sin(angle2) = length3 / sin(angle3)
           angle1 + angle2 + angle3 = 180.
      pi: initial point ([0., 0.], default)
      angle1: first angle from pi clockwise (60., default)
      length1: length of face from pi by angle1 (1., default)
      angle2: second angle from second point (60., default)
      length2: length of face from p2 by angle2 (1., default)
      N: number of points (100, default)
    """
    fname = 'p_angle_triangle'

    angle3 = 180. - angle1 - angle2
    length2 = np.sin(angle2*np.pi/180.)*length1/np.sin(angle1*np.pi/180.)
    length3 = np.sin(angle3*np.pi/180.)*length1/np.sin(angle1*np.pi/180.)

    triangle = np.zeros((N,2), dtype=np.float)

    N3 = int(N/3)
    # first face
    ix = pi[1]
    iy = pi[0]
    dx = length1*np.cos(angle1*np.pi/180.)/(N3-1)
    dy = length1*np.sin(angle1*np.pi/180.)/(N3-1)
    for ip in range(N3):
        triangle[ip,:] = [iy+dy*ip, ix+dx*ip]

    # second face
    ia = -90. - (90.-angle1)
    ix = triangle[N3-1,1]
    iy = triangle[N3-1,0]
    dx = length2*np.cos((ia+angle2)*np.pi/180.)/(N3-1)
    dy = length2*np.sin((ia+angle2)*np.pi/180.)/(N3-1)
    for ip in range(N3):
        triangle[N3+ip,:] = [iy+dy*ip, ix+dx*ip]

    # third face
    N32 = N - 2*N3
    ia = -180. - (90.-angle2)
    ix = triangle[2*N3-1,1]
    iy = triangle[2*N3-1,0]
    angle3 = np.arctan2(pi[0]-iy, pi[1]-ix)
    dx = (pi[1]-ix)/(N32-1)
    dy = (pi[0]-iy)/(N32-1)
    for ip in range(N32):
        triangle[2*N3+ip,:] = [iy+dy*ip, ix+dx*ip]

    return triangle

def p_cross_width(larm=5., width=1., Narms=4, N=200):
    """ Function to draw a cross with arms with a given width and an angle
      larm: legnth of the arms (5., default)
      width: width of the arms (1., default)
      Narms: Number of arms (4, default)
      N: number of points to us (200, default)
    """
    fname = 'p_cross_width'

    Narm = int((N-Narms)/Narms)

    larm2 = larm/2.
    width2 = width/2.

    cross = np.ones((N,2), dtype=np.float)*gen.fillValueF
    da = np.pi/Narms

    N1 = int(Narm*3./8.)
    N2 = int((Narm - 2*N1)/2.)
    N21 = Narm - 2*N1 - N2

    if N2 < 3:
        print errormsg
        print '  ' + fname + ": too few points for ", Narms, " arms !!"
        print "    increase number 'N' at least up to '", 25*Narms
        quit(-1)

    crosssecs = []
    crossdic = {}
    Npot = int(np.log10(Narms))+1

    iip = 0
    for iarm in range(Narms-1):

        a = da*iarm
        iip0 = iip

        # bottom coordinate
        bx = larm*np.cos(a+np.pi)
        by = larm*np.sin(a+np.pi)

        # upper coordinate
        ux = larm*np.cos(a)
        uy = larm*np.sin(a)

        rela = a+np.pi*3./2.
        # SW-NW
        ix = bx + width2*np.cos(rela)
        iy = by + width2*np.sin(rela)
        ex = ux + width2*np.cos(rela)
        ey = uy + width2*np.sin(rela)
        dx = (ex-ix)/(N1-1)
        dy = (ey-iy)/(N1-1)
        for ip in range(N1):
            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
        iip = iip + N1

        # NW-NE
        ix = ex + 0.
        iy = ey + 0.
        ex = ux - width2*np.cos(rela)
        ey = uy - width2*np.sin(rela)
        dx = (ex-ix)/(N2-1)
        dy = (ey-iy)/(N2-1)
        for ip in range(N2):
            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
        iip = iip + N2

        # NW-SW
        ix = ex + 0.
        iy = ey + 0.
        ex = bx - width2*np.cos(rela)
        ey = by - width2*np.sin(rela)
        dx = (ex-ix)/(N1-1)
        dy = (ey-iy)/(N1-1)
        for ip in range(N1):
            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
        iip = iip + N1

        # SW-SE
        ix = ex + 0.
        iy = ey + 0.
        ex = bx + width2*np.cos(rela)
        ey = by + width2*np.sin(rela)
        dx = (ex-ix)/(N21-1)
        dy = (ey-iy)/(N21-1)
        for ip in range(N21):
            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
        iip = iip + N21 + 1

        iarmS = str(iarm).zfill(Npot)
        crosssecs.append(iarmS)
        crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.']

    iip0 = iip

    Narm = N - Narm*(Narms-1) - Narms

    N1 = int(Narm*3./8.)
    N2 = int((Narm - 2*N1)/2.)
    N21 = Narm - 2*N1 - N2

    iarm = Narms-1
    a = da*iarm

    # bottom coordinate
    bx = larm*np.cos(a+np.pi)
    by = larm*np.sin(a+np.pi)

    # upper coordinate
    ux = larm*np.cos(a)
    uy = larm*np.sin(a)

    rela = a+np.pi*3./2.
    # SW-NW
    ix = bx + width2*np.cos(rela)
    iy = by + width2*np.sin(rela)
    ex = ux + width2*np.cos(rela)
    ey = uy + width2*np.sin(rela)
    dx = (ex-ix)/(N1-1)
    dy = (ey-iy)/(N1-1)
    for ip in range(N1):
      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    iip = iip + N1

    # NW-NE
    ix = ex + 0.
    iy = ey + 0.
    ex = ux - width2*np.cos(rela)
    ey = uy - width2*np.sin(rela)
    dx = (ex-ix)/(N2-1)
    dy = (ey-iy)/(N2-1)
    for ip in range(N2):
      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    iip = iip + N2

    # NW-SW
    ix = ex + 0.
    iy = ey + 0.
    ex = bx - width2*np.cos(rela)
    ey = by - width2*np.sin(rela)
    dx = (ex-ix)/(N1-1)
    dy = (ey-iy)/(N1-1)
    for ip in range(N1):
      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    iip = iip + N1

    # SW-SE
    ix = ex + 0.
    iy = ey + 0.
    ex = bx + width2*np.cos(rela)
    ey = by + width2*np.sin(rela)
    dx = (ex-ix)/(N21-1)
    dy = (ey-iy)/(N21-1)
    for ip in range(N21):
      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    iip = iip + N21

    iarmS = str(iarm).zfill(Npot)
    crosssecs.append(iarmS)
    crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.']

    cross = ma.masked_equal(cross, gen.fillValueF)

    return cross, crosssecs, crossdic

####### ####### ##### #### ### ## #
# Plotting

def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10,                   \
  drwsfc=[True,True], colsfc=['#AAAAAA','#646464'],                                  \
  drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.],          \
  drwzline = True, linez=['-.','g',2.], drwxcline=[True,True],                       \
  linexc=[['-','#646400',1.],['--','#646400',1.]],                                   \
  drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]],           \
  drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]):
    """ Function to plot an sphere and determine which standard lines will be also 
        drawn
      iazm: azimut of the camera form the sphere
      iele: elevation of the camera form the sphere
      dist: distance of the camera form the sphere
      Npts: Resolution for the sphere
      radii: radius of the sphere
      drwsfc: whether 'up' and 'down' portions of the sphere should be drawn
      colsfc: colors of the surface of the sphere portions ['up', 'down']
      drwxline: whether x-axis line should be drawn
      linex: properties of the x-axis line ['type', 'color', 'wdith']
      drwyline: whether y-axis line should be drawn
      liney: properties of the y-axis line ['type', 'color', 'wdith']
      drwzline: whether z-axis line should be drawn
      linez: properties of the z-axis line ['type', 'color', 'wdith']
      drwequator: whether 'front' and 'back' portions of the Equator should be drawn
      lineeq: properties of the lines 'front' and 'back' of the Equator
      drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn
      linegw: properties of the lines 'front' and 'back' Greenwhich
      drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn
      linexc: properties of the lines 'front' and 'back' for the 90 line
    """
    fname = 'plot_sphere'

    iazmrad = iazm*np.pi/180.
    ielerad = iele*np.pi/180.

    # 3D surface Sphere
    sfcsphereu, sfcsphered = surface_sphere(radii,Npts)
    
    # greenwhich
    if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.:
        ia=np.pi-ielerad
    else:
        ia=0.-ielerad
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    beta = np.arange(ia,ea+da,da)[0:Npts]
    alpha = np.zeros((Npts), dtype=np.float)
    greenwhichc = spheric_line(radii,alpha,beta)
    ia=ea+0.
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    beta = np.arange(ia,ea+da,da)[0:Npts]
    greenwhichd = spheric_line(radii,alpha,beta)

    # Equator
    ia=np.pi-iazmrad/2.
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    alpha = np.arange(ia,ea+da,da)[0:Npts]
    beta = np.zeros((Npts), dtype=np.float)
    equatorc = spheric_line(radii,alpha,beta)
    ia=ea+0.
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    alpha = np.arange(ia,ea+da,da)[0:Npts]
    equatord = spheric_line(radii,alpha,beta)

    # 90 line
    if iazmrad > np.pi and iazmrad < 2.*np.pi:
        ia=3.*np.pi/2. + ielerad
    else:
        ia=np.pi/2. - ielerad
    if ielerad < 0.:
        ia = ia + np.pi
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    beta = np.arange(ia,ea+da,da)[0:Npts]
    alpha = np.ones((Npts), dtype=np.float)*np.pi/2.
    xclinec = spheric_line(radii,alpha,beta)
    ia=ea+0.
    ea=ia+np.pi
    da = (ea-ia)/(Npts-1)
    beta = np.arange(ia,ea+da,da)[0:Npts]
    xclined = spheric_line(radii,alpha,beta)

    # x line
    xline = np.zeros((2,3), dtype=np.float)
    xline[0,:] = position_sphere(radii, 0., 0.)
    xline[1,:] = position_sphere(radii, np.pi, 0.)

    # y line
    yline = np.zeros((2,3), dtype=np.float)
    yline[0,:] = position_sphere(radii, np.pi/2., 0.)
    yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.)

    # z line
    zline = np.zeros((2,3), dtype=np.float)
    zline[0,:] = position_sphere(radii, 0., np.pi/2.)
    zline[1,:] = position_sphere(radii, 0., -np.pi/2.)

    fig = plt.figure()
    ax = fig.gca(projection='3d')

    # Sphere surface
    if drwsfc[0]:
        ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:],     \
          color=colsfc[0])
    if drwsfc[1]:
        ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:],     \
          color=colsfc[1])

    # greenwhich
    linev = linegw[0]
    if drwgreeenwhich[0]:
        ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0],      \
          color=linev[1], linewidth=linev[2],  label='Greenwich')
    linev = linegw[1]
    if drwgreeenwhich[1]:
        ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0],      \
          color=linev[1], linewidth=linev[2])

    # Equator
    linev = lineeq[0]
    if drwequator[0]:
        ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0],               \
          color=linev[1], linewidth=linev[2], label='Equator')
    linev = lineeq[1]
    if drwequator[1]:
        ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0],               \
          color=linev[1], linewidth=linev[2])

    # 90line
    linev = linexc[0]
    if drwxcline[0]:
        ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1],  \
          linewidth=linev[2], label='90-line')
    linev = linexc[1]
    if drwxcline[1]:
        ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1],  \
          linewidth=linev[2])

    # x line
    linev = linex
    if drwxline:
        ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]],                    \
          [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
          label='xline')

    # y line
    linev = liney
    if drwyline:
        ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]],                    \
          [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
          label='yline')

    # z line
    linev = linez
    if drwzline:
        ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]],                    \
          [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
          label='zline')

    plt.legend()

    return fig, ax

def draw_secs(objdic):
    """ Function to draw an object according to its dictionary
      objdic: dictionary with the parts to draw [polygon, ltype, lcol, lw]
    """
    fname = 'draw_secs'

    for secn in objdic.keys():
        secv = objdic[secn]
        poly = secv[0]
        lt = secv[1]
        lc = secv[2]
        lw = secv[3]

        plt.plot(poly[:,1], poly[:,0], lt, color=lc, linewidth=lw)

    return

def paint_filled(objdic, fillsecs):
    """ Function to draw an object filling given sections
      objdic: dictionary of the object
      filesecs: list of sections to be filled
    """
    fname = 'paint_filled'

    Nsecs = len(fillsecs)

    for secn in fillsecs:
        secvals=objdic[secn]
        pvals = secvals[0]
        fillsecs = []
        Nvals = pvals.shape[0]
        # re-sectionning to plot without masked values
        for ip in range(Nvals-1):
            if type(pvals[ip][0]) == type(gen.mamat[1]): fillsecs.append(ip) 
        Nsecs = len(fillsecs)
        iisc = 0
        for isc in range(Nsecs):
            plt.fill(pvals[iisc:fillsecs[isc],1], pvals[iisc:fillsecs[isc],0],       \
              color=secvals[2])
            iisc = fillsecs[isc]+1
        plt.fill(pvals[iisc:Nvals-1,1], pvals[iisc:Nvals-1,0], color=secvals[2])

    return