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

# Python tools to manage netCDF files.
# L. Fita, CIMA. Mrch 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 sci

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
# 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]
# 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
# 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_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
# buoy1: Function to draw a buoy as superposition of prism and section of ball
# band_lighthouse: Function to plot a lighthouse with spiral bands
# 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
# green_buoy1: Function to draw a green mark buoy using buoy1
# isolateddanger_buoy1: Function to draw an isolated danger buoy using buoy1
# 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_sinusiode: 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
# prefchannelport[A/B]_buoy1: Function to draw a preferred channel port system 
#   [A/B] buoy using buoy1
# prefchannelstarboard[A/B]_buoy1: Function to draw a preferred channel starboard 
#   system [A/B] buoy using buoy1
# red_buoy1: Function to draw a red mark buoy using buoy1
# safewater_buoy1: Function to draw a safe water mark buoy using buoy1
# special_buoy1: Function to draw an special mark buoy using buoy1
# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates
# z_boat: Function to define an schematic boat from the z-plane
# zsailing_boat: Function to define an schematic sailing boat from the z-plane with sails
# zisland1: Function to draw an island from z-axis as the union of a series of points by
#   circular segments

## Plotting
# 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
# [north/east/south/west_buoy1: Function to draw a [North/East/South/West] danger buoy using buoy1

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: 
        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

    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)

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

    return rotvecs

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 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., 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)
      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, 'short', 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., 'short', 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]-1):
        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 = 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]
        if keep == 'left':
            if ic == 0:
                add_secpolygon_list(cutpolygon,0,icut[ic],polygon)
                iipc = icut[ic]
            else:
                add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
                #cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:]
                iipc = iipc + dcpt -1
        else:
            cutpolygon.append(ipt[ic])
            add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
            #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]
        print 'Lluis ipt0', ipt[ic], 'cutline[Nadds-ic,:]', cutline[Nadds[ic]-1,:], 'ept', ept[ic]
        if keep == 'left':
            if ic == 0: 
                for ip in range(Nadds[ic]): cutpolygon.append(cutline[ip,:])
            else: 
                for ip in range(Nadds[ic]-1,0,-1): cutpolygon.append(cutline[ip,:])
            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
                add_secpolygon_list(cutpolygon,ecut[ic]-2,N,polygon)
                #cutpolygon[iipc:iipc+N-ecut[ic]+2,:] = polygon[ecut[ic]-2:N,:]
                iipc = iipc + N-ecut[ic]+2
                cutpolygon.append(polygon[0,:])
                iipc = iipc + 1
        else:
            for ip in range(Nadds[ic]-1,0,-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,:])
    Npts = len(rmpolygon)

    rmpolygon = np.array(rmpolygon)
    cutpolygon = rm_consecpt_polygon(rmpolygon)

    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 = 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,:])

            print 'Lluis ', Ncuts, 'ic', Ncuts-1
            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 = ma.masked_equal(finalcutpolygon, gen.fillValueF)

    Npts = finalcutpolygon.shape[0]

    return Npts, finalcutpolygon

def cut_between_xpolygon_old(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]

    ipt = None
    ept = None

    dx = np.zeros((2), dtype=np.float)
    dy = np.zeros((2), dtype=np.float)
    icut = np.zeros((2), dtype=int)
    ecut = np.zeros((2), dtype=int)
    ipt = np.zeros((2,2), dtype=np.float)
    ept = np.zeros((2,2), dtype=np.float)

    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]

    for ic in range(2):
        xval = xvals[ic]
        if type(polygon) == type(gen.mamat):
            # Assuming clockwise polygons
            for ip in range(N-1):
                if not polygon.mask[ip,0]:
                    eep = ip + 1
                    if eep == N: eep = 0
      
                    if (polygon[ip,1] <= xval and polygon[eep,1] > xval) or          \
                      (polygon[ip,1] < xval and polygon[eep,1] >= xval):
                        icut[ic] = ip
                        dx[ic] = polygon[eep,1] - polygon[ip,1]
                        dy[ic] = polygon[eep,0] - polygon[ip,0]
                        dd = xval - polygon[ip,1]
                        ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]

                    if (polygon[ip,1] >= yval and polygon[eep,1] < xval) or          \
                      (polygon[ip,1] > yval and polygon[eep,1] <= xval):
                        ecut[ic] = ip
                        dx[ic] = polygon[eep,1] - polygon[ip,1]
                        dy[ic] = polygon[eep,0] - polygon[ip,0]
                        dd = xval - polygon[ip,1]
                        ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]
        else:
            # Assuming clockwise polygons
            for ip in range(N-1):
                eep = ip + 1
                if eep == N: eep = 0
      
                if (polygon[ip,1] <= xval and polygon[eep,1] > xval) or              \
                  (polygon[ip,1] < xval and polygon[eep,1] >= xval):
                    icut[ic] = ip
                    dx[ic] = polygon[eep,1] - polygon[ip,1]
                    dy[ic] = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    print 'Lluis ip', ip, 'poly:', polygon[ip,:], 'xval:', xval, 'ip+1', polygon[eep,:]
                    print '  dx:', dx, 'dy:', dy, 'dd', dd
                    ipt[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]

                if (polygon[ip,1] >= xval and polygon[eep,1] < xval) or              \
                  (polygon[ip,1] > xval and polygon[eep,1] <= xval):
                    ecut[ic] = ip
                    dx[ic] = polygon[eep,1] - polygon[ip,1]
                    dy[ic] = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    if dx[ic] == 0.:
                        ept[ic,:] = [polygon[eep,0], xval]
                    else:
                        ept[ic,:] = [polygon[ip,0]+dy[ic]*dd/dx[ic], xval]

        if ipt is None or ept is None:
            print errormsg
            print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'

    Npts = icut[1] - icut[0] + Nadd + ecut[0] - ecut[1]
    cutpolygon = np.zeros((Npts,2), dtype=np.float)
    cutpolygon[0,:] = ipt[0,:]
    cutpolygon[1:icut[1]-icut[0]+1,:] = polygon[icut[0]+1:icut[1]+1,:]
    iip = icut[1]-icut[0]

    # cutting lines
    Nadd2 = int(Nadd/2)
    cutlines = np.zeros((2,Nadd2,2), dtype=np.float)

    for ic in range(2):
        print ic, 'Lluis ipt:', ipt[ic,:], 'ept:', ept[ic,:]
        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
        print '    dx:', dx, 'dy', dy
        cutlines[ic,0,:] = ipt[ic,:]
        for ip in range(1,Nadd2-1):
            cutlines[ic,ip,:] = ipt[ic,:] + np.array([dy*ip,dx*ip])
        cutlines[ic,Nadd2-1,:] = ept[ic,:]

    cutpolygon[iip:iip+Nadd2,:] = cutlines[1,:,:]
    iip = iip + Nadd2
    cutpolygon[iip:iip+(ecut[0]-ecut[1]),:] = polygon[ecut[1]+1:ecut[0]+1,:]
    iip = iip + ecut[0]-ecut[1]
    cutpolygon[iip:iip+Nadd2,:] = cutlines[0,::-1,:]

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

    return Npts, cutpolygon

####### ###### ##### #### ### ## #
# 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_sinusiode(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_sinusiode'
    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

# Combined objects
##

# FROM: http://www.photographers1.com/Sailing/NauticalTerms&Nomenclature.html
def zboat(length=10., beam=1., lbeam=0.4, sternbp=0.5):
    """ Function to define an schematic boat from the z-plane
      length: length of the boat (without stern, default 10)
      beam: beam of the boat (default 1)
      lbeam: length at beam (as percentage of length, default 0.4)
      sternbp: beam at stern (as percentage of beam, default 0.5)
    """
    fname = 'zboat'

    bow = np.array([length, 0.])
    maxportside = np.array([length*lbeam, -beam])
    maxstarboardside = np.array([length*lbeam, beam])
    portside = np.array([0., -beam*sternbp])
    starboardside = np.array([0., beam*sternbp])

    # forward section
    fportside = circ_sec(maxportside, bow, length*2)
    fstarboardside = circ_sec(bow, maxstarboardside, length*2)
    # aft section
    aportside = circ_sec(portside, maxportside, length*2)
    astarboardside = circ_sec(maxstarboardside, starboardside, length*2)
    # stern
    stern = circ_sec(starboardside, portside, length*2)

    dpts = stern.shape[0]
    boat = np.zeros((dpts*5,2), dtype=np.float)

    boat[0:dpts,:] = aportside
    boat[dpts:2*dpts,:] = fportside
    boat[2*dpts:3*dpts,:] = fstarboardside
    boat[3*dpts:4*dpts,:] = astarboardside
    boat[4*dpts:5*dpts,:] = stern

    fname = 'boat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) + '_lb' +  \
      str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) + '.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
          ' %\n')
        for ip in range(dpts*5):
            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"


    # Center line extending [fcl] percentage from length on aft and stern
    fcl = 0.15
    centerline = np.zeros((dpts,2), dtype=np.float)
    dl = length*(1.+fcl*2.)/(dpts-1)
    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)

    # correct order of sections
    boatsecs = ['aportside', 'fportside', 'fstarboardside', 'astarboardside',        \
      'stern', 'centerline']

    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
    dicboat = {'fportside': [fportside, '-', '#8A5900', 2.],                         \
      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
      'stern': [stern, '-', '#8A5900', 2.],                                          \
      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
      'centerline': [centerline, '-.', '#AA6464', 1.5]}
    
    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +'.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' +str(sternbp)+'\n')
        for ip in range(dpts*5):
            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"
  
    return boat, boatsecs, dicboat

def zsailing_boat(length=10., beam=1., lbeam=0.4, sternbp=0.5, lmast=0.6, wmast=0.1, \
  hsd=5., msd=5., lheads=0.38, lmains=0.55):
    """ Function to define an schematic sailing boat from the z-plane with sails
      length: length of the boat (without stern, default 10)
      beam: beam of the boat (default 1)
      lbeam: length at beam (as percentage of length, default 0.4)
      sternbp: beam at stern (as percentage of beam, default 0.5)
      lmast: position of the mast (as percentage of length, default 0.6)
      wmast: width of the mast (default 0.1)
      hsd: head sail direction respect to center line (default 5., -999.99 for upwind)
      msd: main sail direction respect to center line (default 5., -999.99 for upwind)
      lheads: length of head sail (as percentage of legnth, defaul 0.38)
      lmains: length of main sail (as percentage of legnth, defaul 0.55)
    """
    fname = 'zsailing_boat'

    bow = np.array([length, 0.])
    maxportside = np.array([length*lbeam, -beam])
    maxstarboardside = np.array([length*lbeam, beam])
    portside = np.array([0., -beam*sternbp])
    starboardside = np.array([0., beam*sternbp])

    aportside = circ_sec(portside, maxportside, length*2)
    fportside = circ_sec(maxportside, bow, length*2)
    fstarboardside = circ_sec(bow, maxstarboardside, length*2)
    astarboardside = circ_sec(maxstarboardside, starboardside, length*2)
    stern = circ_sec(starboardside, portside, length*2)
    dpts = fportside.shape[0]

    # correct order of sections
    sailingboatsecs = ['aportside', 'fportside', 'fstarboardside', 'astarboardside', \
      'stern', 'mast', 'hsail', 'msail', 'centerline']

    # forward section

    # aft section
    # stern
    # mast
    mast = p_circle(wmast,N=dpts)
    mast = mast + [length*lmast, 0.]
    # head sails
    lsail = lheads*length
    if hsd != -999.99:
        sailsa = np.pi/2. - np.pi*hsd/180.
        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
        endsail[0] = length - endsail[0]
        if bow[1] > endsail[1]:
            hsail = circ_sec(endsail, bow, lsail*2.15)
        else:
            hsail = circ_sec(bow, endsail, lsail*2.15)
    else:
        hsail0 = p_sinusiode(length=lsail, amp=0.2, lamb=0.75, N=dpts)
        hsail = np.zeros((dpts,2), dtype=np.float)
        hsail[:,0] = hsail0[:,1]
        hsail[:,1] = hsail0[:,0]
        hsail = bow - hsail

    # main sails
    lsail = lmains*length
    if msd != -999.99:
        sailsa = np.pi/2. - np.pi*msd/180.
        begsail = np.array([length*lmast, 0.])
        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
        endsail[0] = length*lmast - endsail[0]
        if endsail[1] > begsail[1]:
            msail = circ_sec(begsail, endsail, lsail*2.15)
        else:
            msail = circ_sec(endsail, begsail, lsail*2.15)
    else:
        msail0 = p_sinusiode(length=lsail, amp=0.25, lamb=1., N=dpts)
        msail = np.zeros((dpts,2), dtype=np.float)
        msail[:,0] = msail0[:,1]
        msail[:,1] = msail0[:,0]
        msail = [length*lmast,0] - msail

    sailingboat = np.zeros((dpts*8+4,2), dtype=np.float)

    sailingboat[0:dpts,:] = aportside
    sailingboat[dpts:2*dpts,:] = fportside
    sailingboat[2*dpts:3*dpts,:] = fstarboardside
    sailingboat[3*dpts:4*dpts,:] = astarboardside
    sailingboat[4*dpts:5*dpts,:] = stern
    sailingboat[5*dpts,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[5*dpts+1:6*dpts+1,:] = mast
    sailingboat[6*dpts+1,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[6*dpts+2:7*dpts+2,:] = hsail
    sailingboat[7*dpts+2,:] = [gen.fillValueF, gen.fillValueF]
    sailingboat[7*dpts+3:8*dpts+3,:] = msail
    sailingboat[8*dpts+3,:] = [gen.fillValueF, gen.fillValueF]

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

    # Center line extending [fcl] percentage from length on aft and stern
    fcl = 0.15
    centerline = np.zeros((dpts,2), dtype=np.float)
    dl = length*(1.+fcl*2.)/(dpts-1)
    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)

    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
    dicsailingboat = {'fportside': [fportside, '-', '#8A5900', 2.],                  \
      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
      'stern': [stern, '-', '#8A5900', 2.],                                          \
      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
      'mast': [mast, '-', '#8A5900', 2.], 'hsail': [hsail, '-', '#AAAAAA', 1.],      \
      'msail': [msail, '-', '#AAAAAA', 1.],                                          \
      'centerline': [centerline, '-.', '#AA6464', 1.5]}
    
    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +                \
      '_lm' + str(int(lmast*100.)) + '_wm' + str(int(wmast)) +                       \
      '_hsd' + str(int(hsd)) + '_hs' + str(int(lheads*100.)) +                       \
      '_ms' + str(int(lheads*100.)) + '_msd' + str(int(msd)) +'.dat'
    if not os.path.isfile(fname):
        print infmsg
        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
        of = open(fname, 'w')
        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
          ' % mast position: '+ str(lmast) + ' % mast width: ' + str(wmast) + ' ' +  \
          ' head sail direction:' + str(hsd) + ' head sail length: ' + str(lheads) + \
          ' %' + ' main sail length' + str(lmains) + ' main sail direction:' +       \
          str(msd) +'\n')
        for ip in range(dpts*5):
            of.write(str(sailingboat[ip,0]) + ' ' + str(sailingboat[ip,1]) + '\n')
        
        of.close()
        print fname + ": Successfull written '" + fname + "' !!"
 
    return sailingboat, sailingboatsecs, dicsailingboat

def zisland1(mainpts= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],\
  [2.8, -0.1], [0.1,-0.6]], dtype=np.float), radfrac=3., N=200):
    """ Function to draw an island from z-axis as the union of a series of points by
        circular segments
      mainpts: 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)
          * default= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],
            [2.8, -0.1], [0.1,-0.6]], dtype=np.float)
      radfrac: multiplicative factor of the distance between consecutive points to 
        draw the circular segment (3., default)
      N: number of points (200, default)
    """
    fname = 'zisland1'

    island1 = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # Coastline
    island1 = join_circ_sec_rand(mainpts, arc='short', pos='left')

    islandsecs = ['coastline']
    islanddic = {'coastline': [island1, '-', '#161616', 2.]}

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

    return island1, islandsecs, islanddic

def buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=300):
    """ Function to draw a buoy as superposition of prism and section of ball
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'buoy1'

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

    N3 = int(N/3/5)
    NNp = 0
    iip = 0
    # left lateral
    ix = -width/2.
    Np = N3
    iy = 0.
    dx = 0.
    dy = height/(Np)
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # left upper
    ix = -width/2.
    iy = height
    dx = (width/2.-bradii*bfrac)/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # ball
    p1 = np.array([height, -bradii*bfrac])
    p2 = np.array([height, bradii*bfrac])
    Np = int(2*N/3)
    buoy[iip:iip+Np,:] = circ_sec(p1, p2, 2.*bradii, 'long', 'left', Np)
    NNp = NNp + Np
    iip = NNp

    # right upper
    ix = bradii*bfrac
    iy = height
    Np = N3
    dx = (width/2.-bradii*bfrac)/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # right lateral
    ix = width/2.
    iy = height
    dx = 0.
    dy = -height/(Np)
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    # Base
    ix = width/2.
    iy = 0.
    Np = N - int(2*N/3) - 4*N3 - 1
    dx = -width/(Np)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    buoy[N-1,:] = buoy[0,:]

    buoysecs = ['base']
    buoydic = {'base': [buoy, '-', 'k', 1.5]}

    return buoy, buoysecs, buoydic

def band_lighthouse(height=10., width=2., hlight=3., bands=3, N=300):
    """ Function to plot a lighthouse with spiral bands
      height: height of the tower (10., default)
      width: width of the tower (2., default)
      hlight: height of the light (3., default)
      bands: number of spiral bands (3, default)
      N: number of points (300, default)
    """
    fname = 'band_lighthouse'

    lighthouse = np.ones((N,2), dtype=np.float)*gen.fillValueF
    lighthousesecs = []
    lighthousedic = {}

    # base Tower
    Nsec = int(0.30*N/7)
    p1=np.array([0., width/2.])
    p2=np.array([0., -width/2.])
    iip = 0
    lighthouse[0:Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec

    # left side
    ix=-width/2.
    iy=0.
    dx = 0.
    dy = height/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec

    # Top Tower
    p1=np.array([height, width/2.])
    p2=np.array([height, -width/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec

    # right side
    ix=width/2.
    iy=height
    dx = 0.
    dy = -height/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec + 1

    Ntower = iip-1
    lighthousesecs.append('tower')
    lighthousedic['tower'] = [lighthouse[0:iip-1], '-', 'k', 1.5]

    # Left light
    p1 = np.array([height, -width*0.8/2.])
    p2 = np.array([height+hlight, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*hlight, Nang=Nsec)
    iip = iip + Nsec
    
    # Top Light
    p1=np.array([height+hlight, width*0.8/2.])
    p2=np.array([height+hlight, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec + 1

    # Right light
    p1 = np.array([height+hlight, width*0.8/2.])
    p2 = np.array([height, width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*hlight, Nang=Nsec)
    iip = iip + Nsec

    # Base Light
    p1=np.array([height, width*0.8/2.])
    p2=np.array([height, -width*0.8/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
    iip = iip + Nsec + 1
    lighthousesecs.append('light')
    lighthousedic['light'] = [lighthouse[Ntower+1:iip-1], '-', '#EEEE00', 1.5]

    # Spiral bands
    hb = height/(2.*bands)
    Nsec2 = (N - Nsec*8 - 3)/bands
    for ib in range(bands-1):
        iband = iip
        Nsec = Nsec2/4
        bandS = 'band' + str(ib).zfill(2)
        # hband
        ix = -width/2.
        iy = hb*ib*2
        dx = 0.
        dy = hb/(Nsec-1)
        for ip in range(Nsec):
            lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
        iip = iip + Nsec
        # uband
        p1 = np.array([hb*(ib*2+1), -width/2.])
        p2 = np.array([hb*(ib*2+2), width/2.])
        lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='right', Nang=Nsec)
        iip = iip + Nsec
        # dband
        ix = width/2.
        iy = hb*(ib*2+2)
        dx = 0.
        dy = -hb/(Nsec-1)
        for ip in range(Nsec):
            lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
        iip = iip + Nsec
        # dband
        p1 = np.array([hb*(ib*2+1), width/2.])
        p2 = np.array([hb*ib*2, -width/2.])
        lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)
        iip = iip + Nsec + 1
        lighthousesecs.append(bandS)
        lighthousedic[bandS] = [lighthouse[iband:iip-1], '-', '#6408AA', 2.]

    ib = bands-1
    Nsec3 = (N - iip - 1)
    Nsec = int(Nsec3/4)
    bandS = 'band' + str(ib).zfill(2)
    # hband
    iband = iip
    ix = -width/2.
    iy = hb*ib*2
    dx = 0.
    dy = hb/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec
    # uband
    p1 = np.array([hb*(ib*2+1), -width/2.])
    p2 = np.array([hb*(ib*2+2), width/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='right', Nang=Nsec)
    iip = iip + Nsec
    # dband
    ix = width/2.
    iy = hb*(2+ib*2)
    dx = 0.
    dy = -hb/(Nsec-1)
    for ip in range(Nsec):
        lighthouse[iip+ip,:] = [iy+dy*ip, ix+dx*ip]
    iip = iip + Nsec
    # dband
    Nsec = N - iip
    p1 = np.array([hb*(1+ib*2), width/2.])
    p2 = np.array([hb*ib*2, -width/2.])
    lighthouse[iip:iip+Nsec,:] = circ_sec(p1, p2, 3*width, pos='left', Nang=Nsec)        
    lighthousesecs.append(bandS)
    lighthousedic[bandS] = [lighthouse[iband:iip-1], '-', '#6408AA', 2.]

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

    return lighthouse, lighthousesecs, lighthousedic

def north_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a North danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'north_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    triu = p_angle_triangle(N=N32)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.1*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    trid = p_angle_triangle(N=N323)
    trib = trid*lsign + [0.,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+1.1*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = cut_ypolygon(buoy1v, yval=Height/2., keep='below')
    Ncut, halfup = cut_ypolygon(buoy1v, yval=Height/2., keep='above')

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

    buoysecs = ['buoy', 'sign1', 'sign2', 'half1', 'half2']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'half1': [halfup, '-', 'k', 1.],    \
      'half2': [halfdown, '-', '#FFFF00', 1.]}

    return buoy, buoysecs, buoydic

def east_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a East danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'east_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    triu = p_angle_triangle(N=N32)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.1*lsign,0.]

    # down
    N323 = N - N32 - N2 - 2

    trid = p_angle_triangle(N=N323)
    trid = mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+0.9*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, halfbtw = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, halfup = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign1', 'sign2', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5],                                     \
      'third1': [halfup, '-', 'k', 1.], 'third2': [halfbtw, '-', '#FFFF00', 1.],     \
      'third3': [halfdown, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def south_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a South danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'south_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    trid = p_angle_triangle(N=N32)
    trid = mirror_polygon(trid, 'x')
    trib = trid*lsign + [0.,-lsign/2.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+2.9*lsign,0.]

    # down
    N323 = N - N32 - N2 - 2
    trid = p_angle_triangle(N=N323)
    trid = mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+0.9*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = cut_ypolygon(buoy1v, yval=Height/2., keep='below')
    Ncut, halfup = cut_ypolygon(buoy1v, yval=Height/2., keep='above')

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

    buoysecs = ['buoy', 'sign1', 'sign2', 'half1', 'half2']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'half1': [halfup, '-', '#FFFF00', 1.], \
      'half2': [halfdown, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def west_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.7, N=300):
    """ Function to draw a West danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.7, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'east_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns

    # down
    N32 = int(N3/2) 
    trid = p_angle_triangle(N=N32)
    trid = mirror_polygon(trid, 'x')
    trib = trid*lsign + [lsign,-lsign/2.] 
    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+1.9*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    triu = p_angle_triangle(N=N323)
    trib = triu*lsign + [0.,-lsign/2.] 

    buoy[N2+N323+2:N,:] = trib + [bottsigns+1.*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, halfdown = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, halfbtw1 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, halfup = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign1', 'sign2', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'third1': [halfdown, '-', '#FFFF00', 1.], 'third2': [halfbtw1, '-', 'k', 1.],  \
      'third3': [halfup, '-', '#FFFF00', 1.],                                        \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5]}

    return buoy, buoysecs, buoydic

def safewater_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a safe water mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'safewater_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns
    
    Height = np.max(buoy1v[:,0])
    sign = p_circle(lsign, N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    ix = -width/2.
    Ncut, quarter1 = cut_xpolygon(buoy1v, xval=ix+width/4., keep='left')
    Ncut, quarter2 = cut_between_xpolygon(buoy1v, xval1=ix+width/4., xval2=ix+width/2.)
    Ncut, quarter3 = cut_between_xpolygon(buoy1v, xval1=ix+width/2., xval2=ix+3.*width/4.)
    Ncut, quarter4 = cut_xpolygon(buoy1v, xval=ix+3.*width/4., keep='right')

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

    buoysecs = ['buoy', 'sign', 'quarter1', 'quarter2', 'quarter3', 'quarter4']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'quarter1': [quarter1,'-','r',1.], \
      'quarter2': [quarter2,'-','#FFFFFF',1.], 'quarter3': [quarter3,'-','r',1.],    \
      'quarter4': [quarter4,'-','#FFFFFF',1.]}

    return buoy, buoysecs, buoydic

def red_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a red mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'red_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5]}

    return buoy, buoysecs, buoydic

def green_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, N=300):
    """ Function to draw a green mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'green_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelportA_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, \
  N=300):
    """ Function to draw a preferred channel port system A buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelportA_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    Ncut, third1 = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5], 'third1': [third1,'-','g',1.5],    \
      'third2': [third2,'-','r',1.5], 'third3': [third3,'-','g',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelportB_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.3, \
  N=300):
    """ Function to draw a preferred channel port system B buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelportB_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    triu = p_angle_triangle(N=N3)
    sign = triu*lsign
    buoy[N2+1:N2+2+N3,:] = sign + [Height+0.2*lsign,-lsign/2.]

    # painting it
    Ncut, third1 = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign', 'third1', 'third2', 'third3']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'third1': [third1,'-','r',1.5],    \
      'third2': [third2,'-','g',1.5], 'third3': [third3,'-','r',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelstarboardA_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8,        \
  hsigns=0.3, N=300):
    """ Function to draw a preferred channel starboard system A buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelstarboardA_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    # painting it
    Ncut, third1 = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',1.5], 'third1': [third1,'-','r',1.5],    \
      'third2': [third2,'-','g',1.5], 'third3': [third3,'-','r',1.5]}

    return buoy, buoysecs, buoydic

def prefchannelstarboardB_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8,        \
  hsigns=0.3, N=300):
    """ Function to draw a preferred channel starboard system B buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.3, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'prefchannelstarboardB_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 1
    lsign = height*hsigns*2.
    
    Height = np.max(buoy1v[:,0])
    sign = p_prism(lsign, lsign*2, N=N3)
    buoy[N2+1:N2+2+N3,:] = sign + [Height+1.2*lsign,0.]

    # painting it
    # painting it
    Ncut, third1 = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','g',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','g',1.5], 'third1': [third1,'-','g',1.5],    \
      'third2': [third2,'-','r',1.5], 'third3': [third3,'-','g',1.5]}

    return buoy, buoysecs, buoydic

def isolateddanger_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5,   \
  N=300):
    """ Function to draw an isolated danger buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'isolateddanger_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    # signs
    N3 = N - N2 - 2
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    N32 = int(N3/2) 
    circle = p_circle(lsign/2., N=N32)
    trib = circle + [0.,0.] 

    buoy[N2+1:N2+1+N32,:] = trib + [bottsigns+3.2*lsign,0.]

    # up
    N323 = N - N32 - N2 - 2
    trid = p_circle(lsign/2., N=N32)
    trib = circle + [0.,0.] 
    buoy[N2+N32+2:N,:] = trib + [bottsigns+2.*lsign,0.]

    # painting it
    Height = np.max(buoy1v[:,0])

    Ncut, third1 = cut_ypolygon(buoy1v, yval=Height/3., keep='below')
    Ncut, third2 = cut_between_ypolygon(buoy1v, yval1=Height/3., yval2=Height*2./3.)
    Ncut, third3 = cut_ypolygon(buoy1v, yval=Height*2./3., keep='above')

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

    buoysecs = ['buoy', 'sign1', 'sign2', 'halfk', 'halfy']
    buoydic = {'buoy': [buoy[0:N2,:],'-','k',1.5],                                   \
      'sign1': [buoy[N2+1:N2+N32+1,:],'-','k',1.5],                                  \
      'sign2': [buoy[N2+N32+2:N,:],'-','k',1.5], 'third1': [third1, '-', 'k', 1.],   \
      'third2': [third2, '-', 'r', 1.], 'third3': [third3, '-', 'k', 1.]}

    return buoy, buoysecs, buoydic

def special_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5, N=300):
    """ Function to draw an special mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'special_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    Height = np.max(buoy1v[:,0])

    # sign
    N3 = N - N2 - 1
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    cross, crosssecs, crossdic = p_cross_width(lsign, width=0.3*lsign, Narms=2, N=N3)
    cross = rotate_polygon_2D(cross, 40.05)
    buoy[N2+1:N,:] = cross + [Height+1.1*lsign,0.]

    # painting it
    buoy = ma.masked_equal(buoy, gen.fillValueF)

    buoysecs = ['buoy', 'sign']
    buoydic = {'buoy': [buoy[0:N2,:],'-','#FFFF00',1.5],                             \
      'sign': [buoy[N2+1:N,:],'-','#FFFF00',1.5]}

    return buoy, buoysecs, buoydic

def emergency_buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, hsigns=0.5, N=300):
    """ Function to draw an eergency mark buoy using buoy1
      height: height of the prism (5., default) 
      width: width of the prism (10., default)
      bradii: radii of the ball (1.75, default)
      bfrac: fraction of the ball above the prism (0.8, default)
      hisgns: height of the signs [as reg. triangle] as percentage of the height 
        (0.5, default)
      N: total number of points of the buoy (300, default)
    """
    fname = 'emergency_buoy1'

    buoy = np.ones((N,2), dtype=np.float)*gen.fillValueF

    # buoy
    N2 = int(N/2)
    buoy1v, buoy1vsecs, buoy1vdic = buoy1(height=5., width=10., bradii=1.75,         \
      bfrac=0.8, N=N2)
    buoy[0:N2,:] = buoy1v

    Height = np.max(buoy1v[:,0])

    # sign
    N3 = N - N2 - 1
    
    bottsigns = 2.*bradii+height
    lsign = height*hsigns
    # up
    cross, crosssecs, crossdic = p_cross_width(lsign, width=0.3*lsign, Narms=2, N=N3)
    buoy[N2+1:N,:] = cross + [Height+1.1*lsign,0.]

    # painting it
    ix = -width/2.
    Ncut, fifth1 = cut_xpolygon(buoy1v, xval=ix+width/5., keep='left')
    Ncut, fifth2 = cut_between_xpolygon(buoy1v,xval1=ix+width/5.,xval2=ix+width*2./5.)
    Ncut, fifth3 = cut_between_xpolygon(buoy1v,xval1=ix+width*2./5.,xval2=ix+width*3./5.)
    Ncut, fifth4 = cut_between_xpolygon(buoy1v,xval1=ix+width*3./5.,xval2=ix+width*4./5.)
    Ncut, fifth5 = cut_xpolygon(buoy1v, xval=ix+width*4./5., keep='right')

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

    buoysecs = ['buoy', 'sign', 'fifth1', 'fifth2', 'fifth3', 'fifth4', 'fifth5']
    buoydic = {'buoy': [buoy[0:N2,:],'-','#FFFF00',1.5],                             \
      'sign': [buoy[N2+1:N,:],'-','#FFFF00',1.5],'fifth1':[fifth1,'-','#FFFF00',1.5],\
      'fifth2': [fifth2,'-','#0000FF',1.5],'fifth3': [fifth3,'-','#FFFF00',1.5],     \
      'fifth4': [fifth4,'-','#0000FF',1.5],'fifth5': [fifth5,'-','#FFFF00',1.5]}

    return buoy, buoysecs, buoydic

####### ####### ##### #### ### ## #
# 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 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]
        fillsec = []
        Nvals = pvals.shape[0]
        for ip in range(Nvals-1):
            if type(pvals[ip][0]) != type(gen.mamat[1]) and type(pvals[ip+1][0]) != type(gen.mamat[1]): 
                fillsec.append(pvals[ip,:])
#        plt.fill(pvals[:,1], pvals[:,0], color=secvals[2])
        fillsec = np.array(fillsec)
        plt.fill(fillsec[:,1], fillsec[:,0], color=secvals[2])

    return