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

# 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:
# 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
# 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
# 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
# circ_sec: Function union of point A and B by a section of a circle
# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
# p_angle_triangle: Function to draw a triangle by an initial point and two 
#   consecutive angles and the first length of face. The third angle and 2 and 3rd 
#   face will be computed accordingly the provided values
# p_doubleArrow: Function to provide an arrow with double lines
# p_circle: Function to get a polygon of a circle
# p_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
# band_lighthouse: Function to plot a lighthouse with spiral bands
# red_buoy1: Function to draw a red mark buoy using buoy1
# safewater_buoy1: Function to draw a safe water 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 cut_ypolygon(polygon, yval, keep='bottom', 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 ('bottom', default)
         'bottom': 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 = ['bottom', '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

    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,0] <= yval and polygon[eep,0] >= yval:
                    icut = ip
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ipt = [yval, polygon[ip,1]+dx*dd/dy]

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

            if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
                ecut = ip
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = yval - polygon[ip,0]
                ept = [yval, polygon[ip,1]+dx*dd/dy]

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

    if keep == 'bottom':
        Npts = icut + (N-ecut) + Nadd
        cutpolygon = np.zeros((Npts,2), dtype=np.float)
        if type(polygon) == type(gen.mamat):
            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:]
        else:
            cutpolygon[0:icut+1,:] = polygon[0:icut+1,:]
        iip = icut+1
    else:
        Npts = ecut - icut + Nadd-1
        cutpolygon = np.zeros((Npts,2), dtype=np.float)
        cutpolygon[0,:] = ipt
        cutpolygon[1:ecut-icut,:] = polygon[icut+1:ecut,:]
        iip = ecut-icut-1

    # cutting line
    cutline = np.zeros((Nadd,2), dtype=np.float)
    dx = (ept[1] - ipt[1])/(Nadd-2)
    dy = (ept[0] - ipt[0])/(Nadd-2)
    cutline[0,:] = ipt
    for ip in range(1,Nadd-1):
        cutline[ip,:] = ipt + np.array([dy*ip,dx*ip])
    cutline[Nadd-1,:] = ept
    if keep == 'bottom':
        cutpolygon[iip:iip+Nadd,:] = cutline
        if type(polygon) == type(gen.mamat):
            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
        else:
            cutpolygon[iip+Nadd:Npts,:] = polygon[ecut+1:N,:]
    else:
        cutpolygon[iip:iip+Nadd,:] = cutline[::-1,:]

    rmpolygon = []
    if keep == 'bottom':
        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

    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:
                    icut = ip
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ipt = [polygon[ip,0]+dy*dd/dx, xval]

                if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
                    ecut = ip
                    dx = polygon[eep,1] - polygon[ip,1]
                    dy = polygon[eep,0] - polygon[ip,0]
                    dd = xval - polygon[ip,1]
                    ept = [polygon[ip,0]+dy*dd/dx, 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:
                icut = ip
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = xval - polygon[ip,1]
                ipt = [polygon[ip,0]+dy*dd/dx, xval]

            if polygon[ip,1] >= xval and polygon[eep,1] <= xval:
                ecut = ip
                dx = polygon[eep,1] - polygon[ip,1]
                dy = polygon[eep,0] - polygon[ip,0]
                dd = xval - polygon[ip,1]
                ept = [polygon[ip,0]+dy*dd/dx, xval]

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

    if keep == 'left':
        Npts = icut + (N-ecut) + Nadd + 1
        cutpolygon = np.zeros((Npts,2), dtype=np.float)
        cutpolygon[0,:] = ept
        cutpolygon[1:N-ecut,:] = polygon[ecut+1:N,:]
        iip = N-ecut
        cutpolygon[iip:iip+icut+1,:] = polygon[0:icut+1,:]
        iip = iip + icut+1
    else:
        Npts = ecut - icut + Nadd-1
        cutpolygon = np.zeros((Npts,2), dtype=np.float)
        cutpolygon[0,:] = ipt
        cutpolygon[1:ecut-icut,:] = polygon[icut+1:ecut,:]
        iip = ecut-icut-1

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

    rmpolygon = []
    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)
    cutpolygon = np.array(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]

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

    for ic in range(2):
        yval = yvals[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,0] <= yval and polygon[eep,0] >= yval:
                        icut[ic] = ip
                        dx[ic] = polygon[eep,1] - polygon[ip,1]
                        dy[ic] = polygon[eep,0] - polygon[ip,0]
                        dd = yval - polygon[ip,0]
                        ipt[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]

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

                if polygon[ip,0] >= yval and polygon[eep,0] <= yval:
                    ecut[ic] = ip
                    dx[ic] = polygon[eep,1] - polygon[ip,1]
                    dy[ic] = polygon[eep,0] - polygon[ip,0]
                    dd = yval - polygon[ip,0]
                    ept[ic,:] = [yval, polygon[ip,1]+dx[ic]*dd/dy[ic]]

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

    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):
        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
        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

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]

    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:
                        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:
                        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:
                    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] >= 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]
                    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):
        dx = (ept[ic,1] - ipt[ic,1])/(Nadd2-2)
        dy = (ept[ic,0] - ipt[ic,0])/(Nadd2-2)
        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_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

# 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-1)
    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-1)
    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-1)
    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-1)
    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
    dx = -width/(Np-1)
    dy = 0.
    for ip in range(Np):
        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
    NNp = NNp + Np
    iip = NNp

    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='bottom')
    Ncut, halfup = cut_ypolygon(buoy1v, yval=Height/2., 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], '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='bottom')
    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='bottom')
    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='bottom')
    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', 'quarter1', 'quarter2', 'quarter3', 'quarter4']
    buoydic = {'buoy': [buoy[0:N2,:],'-','r',1.5],                                   \
      'sign': [buoy[N2+1:N2+N3+1,:],'-','r',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]
        plt.fill(pvals[:,1], pvals[:,0], color=secvals[2])

    return