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

Last change on this file since 2845 was 2817, checked in by lfita, 5 years ago

Fixing
File size: 83.8 KB
Line 
1# Python tools to manage netCDF files.
2# L. Fita, CIMA. March 2019
3# More information at: http://www.xn--llusfb-5va.cat/python/PyNCplot
4#
5# pyNCplot and its component geometry_tools.py comes with ABSOLUTELY NO WARRANTY.
6# This work is licendes under a Creative Commons
7#   Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0)
8#
9## Script for geometry calculations and operations as well as definition of different
10###    standard objects and shapes
11
12import numpy as np
13import matplotlib as mpl
14from mpl_toolkits.mplot3d import Axes3D
15import matplotlib.pyplot as plt
16import os
17import generic_tools as gen
18import numpy.ma as ma
19import module_ForSci as fsci
20
21errormsg = 'ERROR -- error -- ERROR -- error'
22infmsg = 'INFORMATION -- information -- INFORMATION -- information'
23
24####### Contents:
25# add_secpolygon_list: Function to add a range of points of a polygon into a list
26# angle_vectors2D: Angle between two vectors with sign
27# cut_between_[x/y]polygon: Function to cut a polygon between 2 given value of the [x/y]-axis
28# cut_[x/y]polygon: Function to cut a polygon from a given value of the [x/y]-axis
29# deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
30# displace_objdic_2D: Function to displace 2D plain the vertices of all polygons of an object
31# dist_points: Function to provide the distance between two points
32# join_circ_sec: Function to join aa series of points by circular segments
33# join_circ_sec_rand: Function to join aa series of points by circular segments with
34#   random perturbations
35# max_coords_poly: Function to provide the extremes of the coordinates of a polygon
36# mirror_polygon: Function to reflex a polygon for a given axis
37# mod_vec: Function to compute the module of a vector
38# position_sphere: Function to tranform fom a point in lon, lat deg coordinates to
39#   cartesian coordinates over an sphere
40# read_join_poly: Function to read an ASCII file with the combination of polygons
41# rm_consecpt_polygon: Function to remove consecutive same point of a polygon
42# rotate_2D: Function to rotate a vector by a certain angle in the plain
43# rotate_objdic_2D: Function to rotate 2D plain the vertices of all polygons of an object
44# rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon
45# rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a
46#   certain angle in the plain
47# rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y
48#   coordinates by a certain angle in the plain
49# spheric_line: Function to transform a series of locations in lon, lat coordinates
50#   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
51# val_consec_between: Function to provide if a given value is between two consecutive ones
52# write_join_poly: Function to write an ASCII file with the combination of polygons
53
54## Shapes/objects
55# circ_sec: Function union of point A and B by a section of a circle
56# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
57# p_angle_triangle: Function to draw a triangle by an initial point and two
58#   consecutive angles and the first length of face. The third angle and 2 and 3rd
59#   face will be computed accordingly the provided values
60# p_doubleArrow: Function to provide an arrow with double lines
61# p_circle: Function to get a polygon of a circle
62# p_cross_width: Function to draw a cross with arms with a given width and an angle
63# p_prism: Function to get a polygon prism
64# p_reg_polygon: Function to provide a regular polygon of Nv vertices
65# p_reg_star: Function to provide a regular star of Nv vertices
66# p_sinusoide: Function to get coordinates of a sinusoidal curve
67# p_square: Function to get a polygon square
68# p_spiral: Function to provide a polygon of an Archimedean spiral
69# p_triangle: Function to provide the polygon of a triangle from its 3 vertices
70# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates
71
72## Plotting
73# draw_secs: Function to draw an object according to its dictionary
74# paint_filled: Function to draw an object filling given sections
75# plot_sphere: Function to plot an sphere and determine which standard lines will be
76#   also drawn
77
78def deg_deci(angle):
79    """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
80      angle: list of [deg, minute, sec] to pass
81    >>> deg_deci([41., 58., 34.])
82    0.732621346072
83    """
84    fname = 'deg_deci'
85
86    deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600.
87
88    if angle[0] < 0.: deg = -deg*np.pi/180.
89    else: deg = deg*np.pi/180.
90
91    return deg
92
93def position_sphere(radii, alpha, beta):
94    """ Function to tranform fom a point in lon, lat deg coordinates to cartesian 
95          coordinates over an sphere
96      radii: radii of the sphere
97      alpha: longitude of the point
98      beta: latitude of the point
99    >>> position_sphere(10., 30., 45.)
100    (0.81031678432964027, -5.1903473778327376, 8.5090352453411846
101    """
102    fname = 'position_sphere'
103
104    xpt = radii*np.cos(beta)*np.cos(alpha)
105    ypt = radii*np.cos(beta)*np.sin(alpha)
106    zpt = radii*np.sin(beta)
107
108    return xpt, ypt, zpt
109
110def spheric_line(radii,lon,lat):
111    """ Function to transform a series of locations in lon, lat coordinates to x,y,z
112          over an 3D space
113      radii: radius of the sphere
114      lon: array of angles along longitudes
115      lat: array of angles along latitudes
116    """
117    fname = 'spheric_line'
118
119    Lint = lon.shape[0]
120    coords = np.zeros((Lint,3), dtype=np.float)
121
122    for iv in range(Lint):
123        coords[iv,:] = position_sphere(radii, lon[iv], lat[iv])
124
125    return coords
126
127def rotate_2D(vector, angle):
128    """ Function to rotate a vector by a certain angle in the plain
129      vector= vector to rotate [y, x]
130      angle= angle to rotate [rad]
131    >>> rotate_2D(np.array([1.,0.]), np.pi/4.)
132    [ 0.70710678 -0.70710678]
133    """
134    fname = 'rotate_2D'
135
136    rotmat = np.zeros((2,2), dtype=np.float)
137
138    rotmat[0,0] = np.cos(angle)
139    rotmat[0,1] = -np.sin(angle)
140    rotmat[1,0] = np.sin(angle)
141    rotmat[1,1] = np.cos(angle)
142
143    rotvector = np.zeros((2), dtype=np.float)
144
145    vecv = np.zeros((2), dtype=np.float)
146
147    # Unifying vector
148    modvec = vector[0]**2+vector[1]**2
149    if modvec != 0 and vector[0] != gen.fillValue: 
150        vecv[0] = vector[1]/modvec
151        vecv[1] = vector[0]/modvec
152
153        rotvec = np.matmul(rotmat, vecv)
154        rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec)
155
156        rotvector[0] = rotvec[1]*modvec
157        rotvector[1] = rotvec[0]*modvec
158    else:
159        rotvector = vector + 0.
160
161    return rotvector
162
163def rotate_polygon_2D(vectors, angle):
164    """ Function to rotate 2D plain the vertices of a polygon
165      line= matrix of vectors to rotate
166      angle= angle to rotate [rad]
167    >>> square = np.zeros((4,2), dtype=np.float)
168    >>> square[0,:] = [-0.5,-0.5]
169    >>> square[1,:] = [0.5,-0.5]
170    >>> square[2,:] = [0.5,0.5]
171    >>> square[3,:] = [-0.5,0.5]
172    >>> rotate_polygon_2D(square, np.pi/4.)
173    [[-0.70710678  0.        ]
174     [ 0.         -0.70710678]
175     [ 0.70710678  0.        ]
176     [ 0.          0.70710678]]
177    """
178    fname = 'rotate_polygon_2D'
179
180    rotvecs = np.zeros(vectors.shape, dtype=np.float)
181
182    mavec = False
183    if type(vectors) == type(gen.mamat):
184        mavec = True
185        vectors = ma.filled(vectors,gen.fillValueF)
186
187    Nvecs = vectors.shape[0]
188    for iv in range(Nvecs):
189        rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
190
191    if mavec: 
192        rotvecs = ma.masked_equal(rotvecs, gen.fillValueF)
193
194    return rotvecs
195
196def displace_objdic_2D(objdic, distance):
197    """ Function to displace 2D plain the vertices of all polygons of an object
198      objdic= dictionary with all the polygons of the object
199      distance= distance to displace [ydist, xdist]
200    """
201    fname = 'displace_objdic_2D'
202
203    disobjdic = dict(objdic)
204
205    for secn in objdic.keys():
206        objv = objdic[secn]
207        vectors = objv[0]
208        lt = objv[1]
209        lc = objv[2]
210        lw = objv[3]
211
212        disvecs = np.zeros(vectors.shape, dtype=np.float)   
213        disvecs = vectors + distance
214        disobjdic[secn] = [disvecs, lt, lc, lw]
215
216    return disobjdic
217
218def rotate_objdic_2D(objdic, angle):
219    """ Function to rotate 2D plain the vertices of all polygons of an object
220      objdic= dictionary with all the polygons of the object
221      angle= angle to rotate [rad]
222    """
223    fname = 'rotate_objdic_2D'
224
225    rotobjdic = dict(objdic)
226
227    for secn in objdic.keys():
228        objv = objdic[secn]
229        vectors = objv[0]
230        lt = objv[1]
231        lc = objv[2]
232        lw = objv[3]
233
234        rotvecs = np.zeros(vectors.shape, dtype=np.float)
235
236        Nvecs = vectors.shape[0]
237        for iv in range(Nvecs):
238            rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
239        rotobjdic[secn] = [rotvecs, lt, lc, lw]
240
241    return rotobjdic
242
243def rotate_line2D(line, angle):
244    """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain
245          angle in the plain
246      line= line to rotate as couple of points [[y0,x0], [y1,x1]]
247      angle= angle to rotate [rad]
248    >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.)
249    [[ 0.          0.        ]
250     [0.70710678  -0.70710678]]
251    """
252    fname = 'rotate_2D'
253
254    rotline = np.zeros((2,2), dtype=np.float)
255    rotline[0,:] = rotate_2D(line[0,:], angle)
256    rotline[1,:] = rotate_2D(line[1,:], angle)
257
258    return rotline
259
260def rotate_lines2D(lines, angle):
261    """ Function to rotate multiple lines given by mulitple pars of x,y coordinates 
262          by a certain angle in the plain
263      line= matrix of N couples of points [N, [y0,x0], [y1,x1]]
264      angle= angle to rotate [rad]
265    >>> square = np.zeros((4,2,2), dtype=np.float)
266    >>> square[0,0,:] = [-0.5,-0.5]
267    >>> square[0,1,:] = [0.5,-0.5]
268    >>> square[1,0,:] = [0.5,-0.5]
269    >>> square[1,1,:] = [0.5,0.5]
270    >>> square[2,0,:] = [0.5,0.5]
271    >>> square[2,1,:] = [-0.5,0.5]
272    >>> square[3,0,:] = [-0.5,0.5]
273    >>> square[3,1,:] = [-0.5,-0.5]
274    >>> rotate_lines2D(square, np.pi/4.)
275    [[[-0.70710678  0.        ]
276      [ 0.         -0.70710678]]
277
278     [[ 0.         -0.70710678]
279      [ 0.70710678  0.        ]]
280
281     [[ 0.70710678  0.        ]
282      [ 0.          0.70710678]]
283
284     [[ 0.          0.70710678]
285      [-0.70710678  0.        ]]]
286    """
287    fname = 'rotate_lines2D'
288
289    rotlines = np.zeros(lines.shape, dtype=np.float)
290
291    Nlines = lines.shape[0]
292    for il in range(Nlines):
293        line = np.zeros((2,2), dtype=np.float)
294        line[0,:] = lines[il,0,:]
295        line[1,:] = lines[il,1,:]
296
297        rotlines[il,:,:] = rotate_line2D(line, angle)
298
299    return rotlines
300
301def dist_points(ptA, ptB):
302    """ Function to provide the distance between two points
303      ptA: coordinates of the point A [yA, xA]
304      ptB: coordinates of the point B [yB, xB]
305    >>> dist_points([1.,1.], [-1.,-1.])
306    2.82842712475
307    """
308    fname = 'dist_points'
309
310    dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2)
311
312    return dist
313
314def mod_vec(vec):
315    """ Function to compute the module of a vector
316      vec: vector [y, x]
317    >>> mod_vec([1., 1.])
318    1.41421356237
319    """
320    fname = 'mod_vec'
321
322    v = np.array(vec, dtype=np.float)
323    vv = v*v
324    mod = np.sqrt(np.sum(vv[:]))
325
326    return mod
327
328def angle_vectors2D(veca, vecb):
329    """ Angle between two vectors with sign
330        FROM: https://stackoverflow.com/questions/5188561/signed-angle-between-two-3d-vectors-with-same-origin-within-the-same-plane
331      veca: angle A [ya, xa]
332      vecb: angle B [yb, xb]
333      NOTE: angle from A to B
334    >>> angle_vectors2D([1.,0.], [0.,1.])
335    1.57079632679
336    >>> angle_vectors2D([0.,1.], [1.,0.])
337    -1.57079632679
338    """
339    fname = 'angle_vectors2D'
340
341    v1 = np.array(veca, dtype=np.float)
342    v2 = np.array(vecb, dtype=np.float)
343
344    moda = mod_vec(v1)
345    modb = mod_vec(v2)
346    modab = mod_vec(v1*v2)
347
348    vc = np.cross(v1,v2)
349    theta = np.arcsin(vc/(moda*modb))
350
351    # Without sign
352    #alpha = np.arccos(modab/(moda*modb))
353
354    return theta
355
356def max_coords_poly(polygon):
357    """ Function to provide the extremes of the coordinates of a polygon
358      polygon: coordinates [Nvertexs, 2] of a polygon
359    >>> square = np.zeros((4,2), dtype=np.float)
360    >>> square[0,:] = [-0.5,-0.5]
361    >>> square[1,:] = [0.5,-0.5]
362    >>> square[2,:] = [0.5,0.5]
363    >>> square[3,:] = [-0.5,0.5]
364    >>> max_coords_poly(square)
365    [-0.5, 0.5], [-0.5, 0.5], [0.5, 0.5], 0.5
366    """
367    fname = 'max_coords_poly'
368
369    # x-coordinate min/max
370    nx = np.min(polygon[:,1])
371    xx = np.max(polygon[:,1])
372
373    # y-coordinate min/max
374    ny = np.min(polygon[:,0])
375    xy = np.max(polygon[:,0])
376
377    # x/y-coordinate maximum of absolute values
378    axx = np.max(np.abs(polygon[:,1]))
379    ayx = np.max(np.abs(polygon[:,0]))
380
381    # absolute maximum
382    xyx = np.max([axx, ayx])
383
384    return [nx, xx], [ny, xy], [ayx, axx], xyx
385
386def mirror_polygon(polygon,axis):
387    """ Function to reflex a polygon for a given axis
388      polygon: polygon to mirror
389      axis: axis at which mirror is located ('x' or 'y')
390    """
391    fname = 'mirror_polygon'
392
393    reflex = np.zeros(polygon.shape, dtype=np.float)
394
395    N = polygon.shape[0]
396    if axis == 'x':
397        for iv in range(N):
398            reflex[iv,:] = [-polygon[iv,0], polygon[iv,1]]
399    elif axis == 'y':
400        for iv in range(N):
401            reflex[iv,:] = [polygon[iv,0], -polygon[iv,1]]
402
403    return reflex
404
405def join_circ_sec(points, radfrac=3., arc='short', side='left', N=200):
406    """ Function to join aa series of points by circular segments
407      points: main points of the island (clockwise ordered, to be joined by circular
408        segments of radii as the radfrac factor of the distance between
409        consecutive points)
410      radfrac: multiplicative factor of the distance between consecutive points to
411        draw the circular segment (3., default)
412      arc: type of arc ('short', default)
413      pos: position of arc ('left', default)
414      N: number of points (200, default)
415    """
416    fname = 'join_circ_sec'
417
418    jcirc_sec = np.ones((N,2), dtype=np.float)
419
420    # main points
421    lpoints = list(points)
422    Npts = len(lpoints)
423    Np = int(N/(Npts+1))
424    for ip in range(Npts-1):
425        p1 = lpoints[ip]
426        p2 = lpoints[ip+1]
427        dps = dist_points(p1, p2)
428        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1,p2,dps*radfrac, arc, side, Np)
429
430    Np2 = N - (Npts-1)*Np
431    p1 = lpoints[Npts-1]
432    p2 = lpoints[0]
433    dps = dist_points(p1, p2)
434    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, side, Np2)
435
436    return jcirc_sec
437
438def join_circ_sec_rand(points, radfrac=3., Lrand=0.1, arc='short', pos='left', N=200):
439    """ Function to join aa series of points by circular segments with random
440        perturbations
441      points: main points of the island (clockwise ordered, to be joined by circular
442        segments of radii as the radfrac factor of the distance between
443        consecutive points)
444      radfrac: multiplicative factor of the distance between consecutive points to
445        draw the circular segment (3., default)
446      Lrand: maximum length of the random perturbation to be added perpendicularly to
447        the direction of the union line between points (0.1, default)
448      arc: type of arc ('short', default)
449      pos: position of arc ('left', default)
450      N: number of points (200, default)
451    """
452    import random
453    fname = 'join_circ_sec_rand'
454
455    jcirc_sec = np.ones((N,2), dtype=np.float)
456
457    # main points
458    lpoints = list(points)
459    Npts = len(lpoints)
460    Np = int(N/(Npts+1))
461    for ip in range(Npts-1):
462        p1 = lpoints[ip]
463        p2 = lpoints[ip+1]
464        dps = dist_points(p1, p2)
465        angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
466        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, arc, pos, Np)
467        drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
468        for iip in range(Np*ip,Np*(ip+1)):
469            jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.)
470
471    Np2 = N - (Npts-1)*Np
472    p1 = lpoints[Npts-1]
473    p2 = lpoints[0]
474    dps = dist_points(p1, p2)
475    angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
476    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., arc, pos, Np2)
477    drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
478    for iip in range(Np*(Npts-1),N):
479        jcirc_sec[iip,:] = jcirc_sec[iip,:] + drand*random.uniform(-1.,1.)
480
481    return jcirc_sec
482
483def write_join_poly(polys, flname='join_polygons.dat'):
484    """ Function to write an ASCII file with the combination of polygons
485      polys: dictionary with the names of the different polygons
486      flname: name of the ASCII file
487    """
488    fname = 'write_join_poly'
489
490    of = open(flname, 'w')
491
492    for polyn in polys.keys():
493        vertices = polys[polyn]
494        Npts = vertices.shape[0]
495        for ip in range(Npts):
496            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')
497
498    of.close()
499
500    return
501
502def read_join_poly(flname='join_polygons.dat'):
503    """ Function to read an ASCII file with the combination of polygons
504      flname: name of the ASCII file
505    """
506    fname = 'read_join_poly'
507
508    of = open(flname, 'r')
509
510    polys = {}
511    polyn = ''
512    poly = []
513    for line in of:
514        if len(line) > 1: 
515            linevals = line.replace('\n','').split(' ')
516            if polyn != linevals[0]:
517                if len(poly) > 1:
518                    polys[polyn] = np.array(poly)
519                polyn = linevals[0]
520                poly = []
521                poly.append([np.float(linevals[2]), np.float(linevals[1])])
522            else:
523                poly.append([np.float(linevals[2]), np.float(linevals[1])])
524
525    of.close()
526    polys[polyn] = np.array(poly)
527
528    return polys
529
530def val_consec_between(valA, valB, val):
531    """ Function to provide if a given value is between two consecutive ones
532      valA: first value
533      valB: second value
534      val: value to determine if it is between
535      >>> val_consec_between(0.5,1.5,0.8)
536      True
537      >>> val_consec_between(0.5,1.5.,-0.8)
538      False
539      >>> val_consec_between(0.5,1.5,0.5)
540      True
541      >>> val_consec_between(-1.58, -1.4, -1.5)
542      True
543      >>> val_consec_between(-1.48747753212, -1.57383530044, -1.5)
544      False
545    """
546    fname = 'val_consec_between'
547
548    btw = False
549    diffA = valA - val
550    diffB = valB - val
551    absdA = np.abs(diffA)
552    absdB = np.abs(diffB)
553    #if (diffA/absdA)* (diffB/absdB) < 0.: btw = True
554#    if valA < 0. and valB < 0. and val < 0.:
555#        if (valA >= val and valB < val) or (valA > val and valB <= val): btw =True
556#    else:
557#        if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True
558    if (valA <= val and valB > val) or (valA < val and valB >= val): btw =True
559
560    return btw
561
562def add_secpolygon_list(listv, iip, eep, polygon):
563    """ Function to add a range of points of a polygon into a list
564      listv: list into which add values of the polygon
565      iip: initial value of the range
566      eep: ending value of the range
567      polygon: array with the points of the polygon
568    """
569    fname = 'add_secpolygon_list'
570
571    if eep > iip:
572        for ip in range(iip,eep): listv.append(polygon[ip,:])
573    else:
574        for ip in range(iip,eep,-1): listv.append(polygon[ip,:])
575
576    return
577
578def rm_consecpt_polygon(polygon):
579    """ Function to remove consecutive same point of a polygon
580      poly: polygon
581    >>> poly = np.ones((5,2), dtype=np.float)
582    >>> poly[2,:] = [2., 1.]
583    rm_consecpt_polygon(poly)
584    [[ 1.  1.]
585     [ 2.  1.]
586     [ 1.  1.]]
587    """
588    fname = 'rm_consecpt_polygon'
589
590    newpolygon = []
591    prevpt = polygon[0,:]
592    newpolygon.append(prevpt)
593    for ip in range(1,polygon.shape[0]):
594        if polygon[ip,0] != prevpt[0] or polygon[ip,1] != prevpt[1]:
595            prevpt = polygon[ip,:]
596            newpolygon.append(prevpt)
597
598    newpolygon = np.array(newpolygon)
599
600    return newpolygon
601
602def cut_ypolygon(polygon, yval, keep='below', Nadd=20):
603    """ Function to cut a polygon from a given value of the y-axis
604      polygon: polygon to cut
605      yval: value to use to cut the polygon
606      keep: part to keep from the height ('below', default)
607         'below': below the height
608         'above': above the height
609      Nadd: additional points to add to draw the line (20, default)
610    """
611    fname = 'cut_ypolygon'
612
613    N = polygon.shape[0]
614    availkeeps = ['below', 'above']
615
616    if not gen.searchInlist(availkeeps, keep):
617        print errormsg
618        print '  ' + fname + ": wring keep '" + keep + "' value !!"
619        print '    available ones:', availkeeps
620        quit(-1)
621
622    ipt = None
623    ept = None
624
625    # There might be more than 1 cut...
626    Ncuts = 0
627    icut = []
628    ecut = []
629    ipt = []
630    ept = []
631
632    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
633      type(gen.mamat.mask[1]):
634        # Assuming clockwise polygons
635        for ip in range(N-1):
636            if not polygon.mask[ip,0]:
637                eep = ip + 1
638                if eep == N: eep = 0
639     
640                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
641                    icut.append(ip)
642                    dx = polygon[eep,1] - polygon[ip,1]
643                    dy = polygon[eep,0] - polygon[ip,0]
644                    dd = yval - polygon[ip,0]
645                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])
646
647                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
648                    ecut.append(ip)
649                    dx = polygon[eep,1] - polygon[ip,1]
650                    dy = polygon[eep,0] - polygon[ip,0]
651                    dd = yval - polygon[ip,0]
652                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
653                    Ncuts = Ncuts + 1
654    else:
655        # Assuming clockwise polygons
656        for ip in range(N-1):
657            eep = ip + 1
658            if eep == N: eep = 0
659     
660            if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
661                icut.append(ip)
662                dx = polygon[eep,1] - polygon[ip,1]
663                dy = polygon[eep,0] - polygon[ip,0]
664                dd = yval - polygon[ip,0]
665                ipt.append([yval, polygon[ip,1]+dx*dd/dy])
666
667            if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
668                ecut.append(ip)
669                dx = polygon[eep,1] - polygon[ip,1]
670                dy = polygon[eep,0] - polygon[ip,0]
671                dd = yval - polygon[ip,0]
672                ept.append([yval, polygon[ip,1]+dx*dd/dy])
673                Ncuts = Ncuts + 1
674
675    # Looking for repeated
676    newicut = icut + []
677    newecut = ecut + []
678    newipt = ipt + []
679    newept = ept + []
680    newNcuts = Ncuts
681    for ic in range(newNcuts-1):
682        for ic2 in range(ic+1,newNcuts):
683            if newipt[ic] == newipt[ic2]:
684                Ncuts = Ncuts-1
685                icut.pop(ic2)
686                ecut.pop(ic2)
687                ipt.pop(ic2)
688                ept.pop(ic2)
689                newNcuts = Ncuts + 0
690
691    if ipt is None or ept is None or Ncuts == 0:
692        print errormsg
693        print '  ' + fname + ': no cutting for polygon at y=', yval, '!!'
694    else:
695        print '  ' + fname + ': found ', Ncuts, ' Ncuts'
696        if Ncuts > 1 and keep == 'below':
697            # Re-shifting cuts by closest distance.
698            xis = []
699            xes = []
700            for ic in range(Ncuts):
701                xp = ipt[ic]
702                xis.append(xp[1])
703                xp = ept[ic]
704                xes.append(xp[1])
705            xs = xis + xes
706            xs.sort()
707            newicut = icut + []
708            newecut = ecut + []
709            newipt = ipt + []
710            newept = ept + []
711            icut = []
712            ecut = []
713            ipt = []
714            ept = []
715            for xv in xs:
716                ic = xis.count(xv)
717                if ic != 0: 
718                    icc = xis.index(xv)
719                    if len(icut) > len(ecut): 
720                        ecut.append(newicut[icc])
721                        ept.append(newipt[icc])
722                    else: 
723                        icut.append(newicut[icc])
724                        ipt.append(newipt[icc])
725                else:
726                    icc = xes.index(xv)
727                    if len(icut) > len(ecut): 
728                        ecut.append(newecut[icc])
729                        ept.append(newept[icc])
730                    else: 
731                        icut.append(newecut[icc])
732                        ipt.append(newept[icc])
733
734#            # Re-shifting cuts. 1st icut --> last ecut; 1st ecut as 1st icut;
735#            #    2nd icut --> last-1 ecut, ....
736#            newicut = icut + []
737#            newecut = ecut + []
738#            newipt = ipt + []
739#            newept = ept + []
740#            for ic in range(Ncuts-1):
741#                ecut[ic] = newecut[Ncuts-ic-1]
742#                ept[ic] = newept[Ncuts-ic-1]
743#                icut[ic+1] = newecut[ic]
744#                ipt[ic+1] = newept[ic]
745
746#            ecut[Ncuts-1] = newicut[Ncuts-1]
747#            ept[Ncuts-1] = newipt[Ncuts-1]
748
749##        print '    yval=', yval, 'cut, ip; ipt ep; ept ________'
750##        for ic in range(Ncuts):
751##            print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic]
752
753    # Length of joining lines
754    Nadds = []
755    if Ncuts > 1:
756        Naddc = (Nadd-Ncuts)/(Ncuts)
757        if Naddc < 3:
758            print errormsg
759            print '  ' + fname + ': too few points for jioning lines !!'
760            print '    increase Nadd at least to:', Ncuts*3+Ncuts
761            quit(-1)
762        for ic in range(Ncuts-1):
763            Nadds.append(Naddc)
764
765        Nadds.append(Nadd-Naddc*(Ncuts-1))
766    else:
767        Nadds.append(Nadd)
768
769    # Total points cut polygon
770    Ntotpts = 0
771    Ncpts = []
772    for ic in range(Ncuts):
773        if keep == 'below':
774            if ic == 0:
775                dpts = icut[ic] + Nadds[ic] + (N - ecut[ic])
776            else:
777                dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
778 
779            # Adding end of the polygon in 'left' keeps
780            if ic == Ncuts - 1: dpts = dpts + N-ecut[ic]
781        else:
782            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
783
784        Ntotpts = Ntotpts + dpts
785        Ncpts.append(ecut[ic] - icut[ic])
786
787    cutpolygon = np.ones((Ntotpts+Ncuts,2), dtype=np.float)*gen.fillValue
788
789    iipc = 0
790    for ic in range(Ncuts):
791        dcpt = Ncpts[ic]
792        if keep == 'below':
793            if ic == 0:
794                cutpolygon[0:icut[ic],:] = polygon[0:icut[ic],:]
795                iipc = icut[ic]
796            else:
797                cutpolygon[iipc:iipc+dcpt-1,:] = polygon[icut[ic]+1:ecut[ic],:]
798                iipc = iipc + dcpt -1
799        else:
800            cutpolygon[iipc,:] = ipt[ic]
801            cutpolygon[iipc:iipc+dcpt-1,:]=polygon[icut[ic]+1:ecut[ic],:]
802            iipc = iipc+dcpt-1
803
804        # cutting line
805        cutline = np.zeros((Nadds[ic],2), dtype=np.float)
806        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
807        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
808        cutline[0,:] = ipt[ic]
809        for ip in range(1,Nadds[ic]-1):
810            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
811        cutline[Nadds[ic]-1,:] = ept[ic]
812        if keep == 'below':
813            if ic == 0: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline
814            else: cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
815            iipc = iipc+Nadds[ic]
816            if ic == 0:
817                cutpolygon[iipc:iipc+N-ecut[ic]-1,:] = polygon[ecut[ic]+1:N,:]
818                iipc = iipc + N-ecut[ic]-1
819                cutpolygon[iipc,:] = polygon[0,:]
820                iipc = iipc + 1
821        else:
822            cutpolygon[iipc:iipc+Nadds[ic],:] = cutline[::-1,:]
823            iipc = iipc+Nadds[ic]
824        iipc = iipc + 1
825
826    rmpolygon = []
827    Npts = cutpolygon.shape[0]
828    if keep == 'below':
829        for ip in range(Npts):
830            if cutpolygon[ip,0] > yval:
831                rmpolygon.append([gen.fillValueF, gen.fillValueF])
832            else:
833                rmpolygon.append(cutpolygon[ip,:])
834    else:
835        for ip in range(Npts):
836            if cutpolygon[ip,0] < yval:
837                rmpolygon.append([gen.fillValueF, gen.fillValueF])
838            else:
839                rmpolygon.append(cutpolygon[ip,:])
840    Npts = len(rmpolygon)
841    cutpolygon = np.array(rmpolygon)
842    cutpolygon = rm_consecpt_polygon(cutpolygon)
843    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
844
845    return Npts, cutpolygon
846
847def cut_xpolygon(polygon, xval, keep='left', Nadd=20):
848    """ Function to cut a polygon from a given value of the x-axis
849      polygon: polygon to cut
850      yval: value to use to cut the polygon
851      keep: part to keep from the value ('left', default)
852         'left': left of the value
853         'right': right of the value
854      Nadd: additional points to add to draw the line (20, default)
855    """
856    fname = 'cut_xpolygon'
857
858    N = polygon.shape[0]
859    availkeeps = ['left', 'right']
860
861    if not gen.searchInlist(availkeeps, keep):
862        print errormsg
863        print '  ' + fname + ": wring keep '" + keep + "' value !!"
864        print '    available ones:', availkeeps
865        quit(-1)
866
867    ipt = None
868    ept = None
869
870    # There might be more than 1 cut ...
871    icut = []
872    ecut = []
873    ipt = []
874    ept = []
875    Ncuts = 0
876    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
877      type(gen.mamat.mask[1]):
878        # Assuming clockwise polygons
879        for ip in range(N-1):
880            if not polygon.mask[ip,1]:
881                eep = ip + 1
882                if eep == N: eep = 0
883     
884                if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
885                    icut.append(ip)
886                    dx = polygon[eep,1] - polygon[ip,1]
887                    dy = polygon[eep,0] - polygon[ip,0]
888                    dd = xval - polygon[ip,1]
889                    ipt.append([polygon[ip,0]+dy*dd/dx, xval])
890
891                if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
892                    ecut.append(ip)
893                    dx = polygon[eep,1] - polygon[ip,1]
894                    dy = polygon[eep,0] - polygon[ip,0]
895                    dd = xval - polygon[ip,1]
896                    ept.append([polygon[ip,0]+dy*dd/dx, xval])
897                    Ncuts = Ncuts + 1
898    else:
899        # Assuming clockwise polygons
900        for ip in range(N-1):
901            eep = ip + 1
902            if eep == N: eep = 0
903     
904            if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
905                icut.append(ip)
906                dx = polygon[eep,1] - polygon[ip,1]
907                dy = polygon[eep,0] - polygon[ip,0]
908                dd = xval - polygon[ip,1]
909                ipt.append([polygon[ip,0]+dy*dd/dx, xval])
910
911            if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
912                ecut.append(ip)
913                dx = polygon[eep,1] - polygon[ip,1]
914                dy = polygon[eep,0] - polygon[ip,0]
915                dd = xval - polygon[ip,1]
916                ept.append([polygon[ip,0]+dy*dd/dx, xval])
917                Ncuts = Ncuts + 1
918
919    # Looking for repeated
920    newicut = icut + []
921    newecut = ecut + []
922    newipt = ipt + []
923    newept = ept + []
924    newNcuts = Ncuts
925    for ic in range(newNcuts-1):
926        for ic2 in range(ic+1,newNcuts):
927            if newipt[ic] == newipt[ic2]:
928                Ncuts = Ncuts-1
929                icut.pop(ic2)
930                ecut.pop(ic2)
931                ipt.pop(ic2)
932                ept.pop(ic2)
933                newNcuts = Ncuts + 0
934
935    if ipt is None or ept is None or Ncuts == 0:
936        print errormsg
937        print '  ' + fname + ': no cutting for polygon at x=', xval, '!!'
938    else:
939        ##print '  ' + fname + ': found ', Ncuts, ' Ncuts'
940        if Ncuts >= 1 and keep == 'left':
941            # Re-shifting cuts by closest heigth.
942            yis = []
943            yes = []
944            for ic in range(Ncuts):
945                yp = ipt[ic]
946                yis.append(yp[0])
947                yp = ept[ic]
948                yes.append(yp[0])
949            ys = yis + yes
950            ys.sort()
951            newicut = icut + []
952            newecut = ecut + []
953            newipt = ipt + []
954            newept = ept + []
955            icut = []
956            ecut = []
957            ipt = []
958            ept = []
959            for yv in ys:
960                ic = yis.count(yv)
961                if ic != 0: 
962                    icc = yis.index(yv)
963                    if len(icut) > len(ecut): 
964                        ecut.append(newicut[icc])
965                        ept.append(newipt[icc])
966                    else: 
967                        icut.append(newicut[icc])
968                        ipt.append(newipt[icc])
969                else:
970                    icc = yes.index(yv)
971                    if len(icut) > len(ecut): 
972                        ecut.append(newecut[icc])
973                        ept.append(newept[icc])
974                    else: 
975                        icut.append(newecut[icc])
976                        ipt.append(newept[icc])
977        #print '    xval=', xval, 'cut, ip; ipt ep; ept ________'
978        #for ic in range(Ncuts):
979        #    print '      ', ic, icut[ic], ';', ipt[ic], ecut[ic], ';', ept[ic]
980
981    # Length of joining lines
982    Nadds = []   
983    if Ncuts > 1:
984        Naddc = (Nadd-Ncuts)/(Ncuts)
985        if Naddc < 3:
986            print errormsg
987            print '  ' + fname + ': too few points for jioning lines !!'
988            print '    increase Nadd at least to:', Ncuts*3+Ncuts
989            quit(-1)
990        for ic in range(Ncuts-1):
991            Nadds.append(Naddc)
992
993        Nadds.append(Nadd-Naddc*(Ncuts-1))
994    else:
995        Nadds.append(Nadd)
996
997   # Total points cut polygon
998    Ntotpts = 0
999    Ncpts = []
1000    for ic in range(Ncuts):
1001        if keep == 'left':
1002            if ic == 0: 
1003                dpts = icut[ic] + Nadds[ic] + (N - ecut[ic])
1004            else:
1005                dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
1006
1007            # Adding end of the polygon in 'left' keeps
1008            if ic == Ncuts - 1: dpts = dpts + N-ecut[ic]
1009        else:
1010            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
1011
1012        Ntotpts = Ntotpts + dpts
1013        Ncpts.append(ecut[ic] - icut[ic])
1014
1015    cutpolygon = []
1016    iipc = 0
1017    for ic in range(Ncuts):
1018        dcpt = Ncpts[ic]
1019        cutpolygon.append(ipt[ic])
1020        if keep == 'left':
1021            if ic == 0:
1022                add_secpolygon_list(cutpolygon,icut[ic]+1,N,polygon)
1023                add_secpolygon_list(cutpolygon,0,ecut[ic],polygon)
1024                iipc = icut[ic]
1025            else:
1026                add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
1027        else:
1028            add_secpolygon_list(cutpolygon,icut[ic]+1,ecut[ic],polygon)
1029            iipc = iipc+dcpt-1
1030        # cutting line
1031        cutline = np.zeros((Nadds[ic],2), dtype=np.float)
1032        dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
1033        dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
1034        cutline[0,:] = ipt[ic]
1035        for ip in range(1,Nadds[ic]-1):
1036            cutline[ip,:] = ipt[ic] + np.array([dy*ip,dx*ip])
1037        cutline[Nadds[ic]-1,:] = ept[ic]
1038        if keep == 'left':
1039            for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:])
1040            iipc = iipc+Nadds[ic]
1041            if ic == 0:
1042                add_secpolygon_list(cutpolygon,ecut[ic],N,polygon)
1043                cutpolygon.append(polygon[0,:])
1044                iipc = iipc + 1
1045        else:
1046            for ip in range(Nadds[ic]-1,-1,-1): cutpolygon.append(cutline[ip,:])
1047            iipc = iipc+Nadds[ic]
1048        cutpolygon.append([gen.fillValueF, gen.fillValueF])
1049        iipc = iipc + 1
1050
1051    cutpolygon = np.array(cutpolygon)
1052    rmpolygon = []
1053    Npts = cutpolygon.shape[0]
1054    if keep == 'left':
1055        for ip in range(Npts):
1056            if cutpolygon[ip,1] > xval:
1057                rmpolygon.append([gen.fillValueF, gen.fillValueF])
1058            else:
1059                rmpolygon.append(cutpolygon[ip,:])
1060    else:
1061        for ip in range(Npts):
1062            if cutpolygon[ip,1] < xval:
1063                rmpolygon.append([gen.fillValueF, gen.fillValueF])
1064            else:
1065                rmpolygon.append(cutpolygon[ip,:])
1066
1067    rmpolygon = np.array(rmpolygon)
1068    cutpolygon = rm_consecpt_polygon(rmpolygon)
1069    Npts = cutpolygon.shape[0]
1070
1071    cutpolygon = ma.masked_equal(cutpolygon, gen.fillValueF)
1072
1073    return Npts, cutpolygon
1074
1075def cut_between_ypolygon(polygon, yval1, yval2, Nadd=20):
1076    """ Function to cut a polygon between 2 given value of the y-axis
1077      polygon: polygon to cut
1078      yval1: first value to use to cut the polygon
1079      yval2: first value to use to cut the polygon
1080      Nadd: additional points to add to draw the line (20, default)
1081    """
1082    fname = 'cut_betwen_ypolygon'
1083
1084    N = polygon.shape[0]
1085
1086    if yval1 > yval2:
1087        print errormsg
1088        print '  ' + fname + ': wrong between cut values !!'
1089        print '     it is expected yval1 < yval2'
1090        print '     values provided yval1: (', yval1, ')> yval2 (', yval2, ')'
1091        quit(-1)
1092
1093    yvals = [yval1, yval2]
1094
1095    ipt = None
1096    ept = None
1097
1098    cuts = {}
1099    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
1100      type(gen.mamat.mask[1]):
1101        for ic in range(2):
1102            yval = yvals[ic]
1103            # There might be more than 1 cut ...
1104            icut = []
1105            ecut = []
1106            ipt = []
1107            ept = []
1108            Ncuts = 0
1109            # Assuming clockwise polygons
1110            for ip in range(N-1):
1111                if not polygon.mask[ip,0]:
1112                    eep = ip + 1
1113                    if eep == N: eep = 0
1114     
1115                    if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
1116                        icut.append(ip)
1117                        dx = polygon[eep,1] - polygon[ip,1]
1118                        dy = polygon[eep,0] - polygon[ip,0]
1119                        dd = yval - polygon[ip,0]
1120                        ipt.append([yval, polygon[ip,1]+dx*dd/dy])
1121
1122                    if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
1123                        ecut.append(ip)
1124                        dx = polygon[eep,1] - polygon[ip,1]
1125                        dy = polygon[eep,0] - polygon[ip,0]
1126                        dd = yval - polygon[ip,0]
1127                        ept.append([yval, polygon[ip,1]+dx*dd/dy])
1128                        Ncuts = Ncuts + 1
1129
1130            # Looking for repeated
1131            newicut = icut + []
1132            newecut = ecut + []
1133            newipt = ipt + []
1134            newept = ept + []
1135            newNcuts = Ncuts
1136            for icp in range(newNcuts-1):
1137                for ic2 in range(icp+1,newNcuts):
1138                    if newipt[icp] == newipt[ic2]:
1139                        Ncuts = Ncuts-1
1140                        icut.pop(ic2)
1141                        ecut.pop(ic2)
1142                        ipt.pop(ic2)
1143                        ept.pop(ic2)
1144                        newNcuts = Ncuts + 0
1145
1146            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
1147    else:
1148        for ic in range(2):
1149            yval = yvals[ic]
1150            # There might be more than 1 cut ...
1151            icut = []
1152            ecut = []
1153            ipt = []
1154            ept = []
1155            Ncuts = 0
1156            # Assuming clockwise polygons
1157            for ip in range(N-1):
1158                eep = ip + 1
1159                if eep == N: eep = 0
1160     
1161                if val_consec_between(polygon[ip,0], polygon[eep,0], yval):
1162                    icut.append(ip)
1163                    dx = polygon[eep,1] - polygon[ip,1]
1164                    dy = polygon[eep,0] - polygon[ip,0]
1165                    dd = yval - polygon[ip,0]
1166                    ipt.append([yval, polygon[ip,1]+dx*dd/dy])
1167
1168                if val_consec_between(polygon[eep,0], polygon[ip,0], yval):
1169                    ecut.append(ip)
1170                    dx = polygon[eep,1] - polygon[ip,1]
1171                    dy = polygon[eep,0] - polygon[ip,0]
1172                    dd = yval - polygon[ip,0]
1173                    ept.append([yval, polygon[ip,1]+dx*dd/dy])
1174                    Ncuts = Ncuts + 1
1175            # Looking for repeated
1176            newicut = icut + []
1177            newecut = ecut + []
1178            newipt = ipt + []
1179            newept = ept + []
1180            newNcuts = Ncuts
1181            for icp in range(newNcuts-1):
1182                for ic2 in range(icp+1,newNcuts):
1183                    if newipt[icp] == newipt[ic2]:
1184                        Ncuts = Ncuts-1
1185                        icut.pop(ic2)
1186                        ecut.pop(ic2)
1187                        ipt.pop(ic2)
1188                        ept.pop(ic2)
1189                        newNcuts = Ncuts + 0
1190
1191            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
1192
1193    Naddlines = {}
1194    for icc in range(2):
1195        cutv = cuts[icc]
1196        Ncuts = cutv[4]
1197        # Length of joining lines
1198        Nadds = []
1199        if Ncuts > 1:
1200            Naddc = (Nadd-Ncuts)/(Ncuts)
1201            if Naddc < 3:
1202                print errormsg
1203                print '  ' + fname + ': too few points for jioning lines !!'
1204                print '    increase Nadd at least to:', Ncuts*3+Ncuts
1205                quit(-1)
1206            for ic in range(Ncuts-1):
1207                Nadds.append(Naddc)
1208
1209            Nadds.append(Nadd-Naddc*(Ncuts-1))
1210        else:
1211            Nadds.append(Nadd)
1212
1213        # Total points cut polygon
1214        Ntotpts = 0
1215        Ncpts = []
1216        for ic in range(Ncuts):
1217            dpts = ecut[ic] - icut[ic] + Nadds[ic] - 1
1218
1219            Ntotpts = Ntotpts + dpts
1220            Ncpts.append(ecut[ic] - icut[ic])
1221
1222        Naddlines[icc] = [Nadds, Ntotpts, Ncpts]
1223
1224    cutv1 = cuts[0]
1225    addv1 = Naddlines[0]
1226    Nadds1 = addv1[0]
1227    Ncuts1 = cutv1[4]
1228
1229    cutv2 = cuts[1]
1230    addv2 = Naddlines[1]
1231    Nadds2 = addv2[0]
1232    Ncuts2 = cutv2[4]
1233
1234    if Ncuts1 != Ncuts2: 
1235        print errormsg
1236        print '  ' + fname + ": different number of cuts !!"
1237        print '    yval1:', yval1, 'Ncuts=', Ncuts1
1238        print '    yval2:', yval2, 'Ncuts=', Ncuts2
1239        print '      I am not prepare to deal with it'
1240        quit(-1)
1241    #else:
1242    #    print '  ' + fname + ' _______'
1243    #    print '    yval1:', yval1, 'Ncuts=', Ncuts1
1244    #    print '    yval2:', yval2, 'Ncuts=', Ncuts2
1245
1246    icut1 = cutv1[0]
1247    ecut1 = cutv1[1]
1248    ipt1 = cutv1[2]
1249    ept1 = cutv1[3]
1250    icut2 = cutv2[0]
1251    ecut2 = cutv2[1]
1252    ipt2 = cutv2[2]
1253    ept2 = cutv2[3]
1254
1255    # Looking for pairs of cuts. Grouping for smallest x distance between initial   
1256    #   points of each cut
1257    cutpolygons = []
1258    for iic1 in range(Ncuts1):
1259        iip = 0
1260        cutpolygon = []
1261        ic1 = icut1[iic1]
1262        ec1 = ecut1[iic1]
1263        ip1 = ipt1[iic1]
1264        ep1 = ept1[iic1]
1265
1266        ipx2s = []
1267        for ip in range(Ncuts2): 
1268            ip2 = ipt2[ip]
1269            ipx2s.append(ip2[1])
1270        dxps = ipx2s - ip1[1]
1271        dxps = np.where(dxps < 0., gen.fillValueF, dxps)
1272        ndxps = np.min(dxps)
1273        iic12 = gen.index_vec(dxps,ndxps)
1274
1275        ic2 = icut2[iic12]
1276        ec2 = ecut2[iic12]
1277        ip2 = ipt2[iic12]
1278        ep2 = ept2[iic12]
1279
1280        #print 'Lluis iic1', iic1, 'ic1', ic1, 'ec1', ec1, 'ipt1', ip1, 'ept1', ep1, 'Nadds1', Nadds1
1281        #print '    iic12', iic12, 'ic2', ic2, 'ec2', ec2, 'ipt2', ip2, 'ept2', ep2, 'Nadds2', Nadds2
1282
1283        cutpolygon.append(ip1)
1284        for ip in range(ic1+1,ic2-1):
1285            cutpolygon.append(polygon[ip,:])
1286        iip = ic2-ic1
1287        # cutting line 1
1288        Nadd2 = Nadds1[iic1]
1289        cutlines = np.zeros((Nadd2,2), dtype=np.float)
1290        dx = (ep2[1] - ip2[1])/(Nadd2-2)
1291        dy = (ep2[0] - ip2[0])/(Nadd2-2)
1292        cutlines[0,:] = ip2
1293        for ip in range(1,Nadd2-1):
1294            cutlines[ip,:] = ip2 + np.array([dy*ip,dx*ip])
1295        cutlines[Nadd2-1,:] = ep2
1296        for ip in range(Nadd2): cutpolygon.append(cutlines[ip,:])
1297        iip = iip + Nadd2
1298
1299        for ip in range(ec2,ec1):
1300            cutpolygon.append(polygon[ip,:])
1301        iip = iip + ec1-ec2
1302        # cutting line 2
1303        Nadd2 = Nadds2[iic12]
1304        cutlines = np.zeros((Nadd2,2), dtype=np.float)
1305        dx = (ep1[1] - ip1[1])/(Nadd2-2)
1306        dy = (ep1[0] - ip1[0])/(Nadd2-2)
1307        cutlines[0,:] = ip1
1308        for ip in range(1,Nadd2-1):
1309            cutlines[ip,:] = ip1 + np.array([dy*ip,dx*ip])
1310        cutlines[Nadd2-1,:] = ep1
1311        for ip in range(Nadd2-1,0,-1):
1312            cutpolygon.append(cutlines[ip,:])
1313
1314        cutpolygon.append(ip1)
1315
1316        cutpolygon.append([gen.fillValueF,gen.fillValueF])
1317        if len(cutpolygons) == 0: cutpolygons = cutpolygon
1318        else: cutpolygons = cutpolygons + cutpolygon
1319
1320    cutpolygons = np.array(cutpolygons)
1321    cutpolygons = rm_consecpt_polygon(cutpolygons)
1322    cutpolygons = ma.masked_equal(cutpolygons, gen.fillValueF)
1323
1324    Npts = cutpolygons.shape[0]
1325
1326    return Npts, cutpolygons
1327
1328def cut_between_xpolygon(polygon, xval1, xval2, Nadd=20):
1329    """ Function to cut a polygon between 2 given value of the x-axis
1330      polygon: polygon to cut
1331      xval1: first value to use to cut the polygon
1332      xval2: first value to use to cut the polygon
1333      Nadd: additional points to add to draw the line (20, default)
1334    """
1335    fname = 'cut_betwen_xpolygon'
1336
1337    N = polygon.shape[0]
1338
1339    if xval1 > xval2:
1340        print errormsg
1341        print '  ' + fname + ': wrong between cut values !!'
1342        print '     it is expected xval1 < xval2'
1343        print '     values provided xval1: (', xval1, ')> xval2 (', xval2, ')'
1344        quit(-1)
1345
1346    xvals = [xval1, xval2]
1347
1348    ipt = None
1349    ept = None
1350
1351    cuts = {}
1352    if type(polygon) == type(gen.mamat) and type(polygon.mask) !=                    \
1353      type(gen.mamat.mask[1]):
1354        for ic in range(2):
1355            xval = xvals[ic]
1356            # There might be more than 1 cut ...
1357            icut = []
1358            ecut = []
1359            ipt = []
1360            ept = []
1361            Ncuts = 0
1362            # Assuming clockwise polygons
1363            for ip in range(N-1):
1364                if not polygon.mask[ip,0]:
1365                    eep = ip + 1
1366                    if eep == N: eep = 0
1367     
1368                    if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
1369                        icut.append(ip)
1370                        dx = polygon[eep,1] - polygon[ip,1]
1371                        dy = polygon[eep,0] - polygon[ip,0]
1372                        dd = xval - polygon[ip,1]
1373                        ipt.append([polygon[ip,0]+dy*dd/dx, xval])
1374
1375                    if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
1376                        ecut.append(ip)
1377                        dx = polygon[eep,1] - polygon[ip,1]
1378                        dy = polygon[eep,0] - polygon[ip,0]
1379                        dd = xval - polygon[ip,1]
1380                        ept.append([polygon[ip,0]+dy*dd/dx, xval])
1381                        Ncuts = Ncuts + 1
1382
1383            # Looking for repeated
1384            newicut = icut + []
1385            newecut = ecut + []
1386            newipt = ipt + []
1387            newept = ept + []
1388            newNcuts = Ncuts
1389            for icp in range(newNcuts-1):
1390                for ic2 in range(icp+1,newNcuts):
1391                    if newipt[icp] == newipt[ic2]:
1392                        Ncuts = Ncuts-1
1393                        icut.pop(ic2)
1394                        ecut.pop(ic2)
1395                        ipt.pop(ic2)
1396                        ept.pop(ic2)
1397                        newNcuts = Ncuts + 0
1398
1399            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
1400    else:
1401        for ic in range(2):
1402            xval = xvals[ic]
1403            # There might be more than 1 cut ...
1404            icut = []
1405            ecut = []
1406            ipt = []
1407            ept = []
1408            Ncuts = 0
1409            # Assuming clockwise polygons
1410            for ip in range(N-1):
1411                eep = ip + 1
1412                if eep == N: eep = 0
1413     
1414                if val_consec_between(polygon[ip,1], polygon[eep,1], xval):
1415                    icut.append(ip)
1416                    dx = polygon[eep,1] - polygon[ip,1]
1417                    dy = polygon[eep,0] - polygon[ip,0]
1418                    dd = xval - polygon[ip,1]
1419                    ipt.append([polygon[ip,0]+dy*dd/dx, xval])
1420
1421                if val_consec_between(polygon[eep,1], polygon[ip,1], xval):
1422                    ecut.append(ip)
1423                    dx = polygon[eep,1] - polygon[ip,1]
1424                    dy = polygon[eep,0] - polygon[ip,0]
1425                    dd = xval - polygon[ip,1]
1426                    ept.append([polygon[ip,0]+dy*dd/dx, xval])
1427                    Ncuts = Ncuts + 1
1428            # Looking for repeated
1429            newicut = icut + []
1430            newecut = ecut + []
1431            newipt = ipt + []
1432            newept = ept + []
1433            newNcuts = Ncuts
1434            for icp in range(newNcuts-1):
1435                for ic2 in range(icp+1,newNcuts):
1436                    if newipt[icp] == newipt[ic2]:
1437                        Ncuts = Ncuts-1
1438                        icut.pop(ic2)
1439                        ecut.pop(ic2)
1440                        ipt.pop(ic2)
1441                        ept.pop(ic2)
1442                        newNcuts = Ncuts + 0
1443
1444            cuts[ic] = [icut, ecut, ipt, ept, Ncuts]
1445
1446    for iic in range(1):
1447        cutvs = cuts[iic]
1448        icut = cutvs[0]
1449        ecut = cutvs[1]
1450        ipt = cutvs[2]
1451        ept = cutvs[3]
1452        Ncuts = cutvs[4]
1453        if Ncuts > 0:
1454            # Re-shifting cuts by closest heigth.
1455            yis = []
1456            yes = []
1457            for ic in range(Ncuts):
1458                yp = ipt[ic]
1459                yis.append(yp[0])
1460                yp = ept[ic]
1461                yes.append(yp[0])
1462            ys = yis + yes
1463            ys.sort()
1464            newicut = icut + []
1465            newecut = ecut + []
1466            newipt = ipt + []
1467            newept = ept + []
1468            icut = []
1469            ecut = []
1470            ipt = []
1471            ept = []
1472            for yv in ys:
1473                ic = yis.count(yv)
1474                if ic != 0: 
1475                    icc = yis.index(yv)
1476                    if len(icut) > len(ecut): 
1477                        ecut.append(newicut[icc])
1478                        ept.append(newipt[icc])
1479                    else: 
1480                        icut.append(newicut[icc])
1481                        ipt.append(newipt[icc])
1482                else:
1483                    icc = yes.index(yv)
1484                    if len(icut) > len(ecut): 
1485                        ecut.append(newecut[icc])
1486                        ept.append(newept[icc])
1487                    else: 
1488                        icut.append(newecut[icc])
1489                        ipt.append(newept[icc])
1490
1491        cuts[iic] = [icut, ecut, ipt, ept, Ncuts]
1492
1493    Naddlines = {}
1494    for icc in range(2):
1495        cutv = cuts[icc]
1496        Ncuts = cutv[4]
1497        if Ncuts == 0:
1498            print errormsg
1499            print '  ' + fname + ": no cuts for xval=", xvals[icc], '!!'
1500            quit(-1)
1501        #print '  icc:', icc, 'ic ec ipt ept _______'
1502        #for ic in range(Ncuts):
1503        #    print ic, ':', cutv[0][ic], cutv[1][ic], cutv[2][ic], cutv[3][ic]
1504
1505        # Length of joining lines
1506        Nadds = []
1507        if Ncuts > 1:
1508            Naddc = (Nadd-Ncuts)/(Ncuts)
1509            if Naddc < 3:
1510                print errormsg
1511                print '  ' + fname + ': too few points for jioning lines !!'
1512                print '    increase Nadd at least to:', Ncuts*3+Ncuts
1513                quit(-1)
1514            for ic in range(Ncuts-1):
1515                Nadds.append(Naddc)
1516
1517            Nadds.append(Nadd-Naddc*(Ncuts-1))
1518        else:
1519            Nadds.append(Nadd)
1520
1521        Naddlines[icc] = Nadds
1522
1523    # sides
1524    sides = {}
1525    for iic in range(2):
1526        cutvs = cuts[iic]
1527        icut = cutvs[0]
1528        ecut = cutvs[1]
1529        ipt = cutvs[2]
1530        ept = cutvs[3]
1531        Ncuts = cutvs[4]
1532        Nadds = Naddlines[iic]
1533        cutpolygon = []
1534        # left side
1535        if iic == 0:
1536            for ic in range(Ncuts-1):
1537                cutpolygon.append(ipt[ic])
1538                dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
1539                dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
1540                for ip in range(1,Nadds[ic]-1):
1541                    cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
1542                cutpolygon.append(ept[ic])
1543                for ip in range(ecut[ic]+1,icut[ic+1]): cutpolygon.append(polygon[ip,:])
1544
1545            ic = Ncuts-1
1546            cutpolygon.append(ipt[ic])
1547            dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
1548            dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
1549            for ip in range(1,Nadds[ic]-1):
1550                cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
1551        # right side
1552        else:
1553            for ic in range(Ncuts-1):
1554                cutpolygon.append(ipt[ic])
1555
1556                # line
1557                dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
1558                dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
1559                for ip in range(1,Nadds[ic]-1):
1560                    cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
1561                cutpolygon.append(ept[ic])
1562                for ip in range(ecut[ic],icut[ic+1]): cutpolygon.append(polygon[ip,:])
1563
1564            ic = Ncuts-1
1565            cutpolygon.append(ipt[ic])
1566            # line
1567            dx = (ept[ic][1] - ipt[ic][1])/(Nadds[ic]-1)
1568            dy = (ept[ic][0] - ipt[ic][0])/(Nadds[ic]-1)
1569            for ip in range(1,Nadds[ic]-1):
1570                cutpolygon.append([ipt[ic][0]+dy*ip, ipt[ic][1]+dx*ip])
1571            cutpolygon.append(ept[ic])
1572        sides[iic] = cutpolygon
1573       
1574    # joining sides by e1[Ncuts1-1] --> i2[0]; e2[Ncuts2-1] --> i1[0]
1575    cutv1 = cuts[0]
1576    Ncuts1 = cutv1[4]
1577    ec1 = cutv1[1][np.max([0,Ncuts1-1])]
1578    ic1 = cutv1[0][0]
1579    ept1 = cutv1[3][np.max([0,Ncuts1-1])]
1580    ipt1 = cutv1[2][0]
1581
1582    cutv2 = cuts[1]
1583    Ncuts2 = cutv2[4]
1584    ec2 = cutv2[1][np.max([0,Ncuts2-1])]
1585    ic2 = cutv2[0][0]
1586    ept2 = cutv2[3][np.max([0,Ncuts2-1])]
1587    ipt2 = cutv2[2][0]
1588
1589    finalcutpolygon = sides[0]
1590    for ip in range(ec1+1,ic2): finalcutpolygon.append(polygon[ip,:])
1591    finalcutpolygon = finalcutpolygon + sides[1]
1592    for ip in range(ec2+1,ic1): finalcutpolygon.append(polygon[ip,:])
1593    finalcutpolygon.append(ipt1)
1594
1595    finalcutpolygon = np.array(finalcutpolygon)
1596
1597    finalcutpolygon = rm_consecpt_polygon(finalcutpolygon)
1598    finalcutpolygon = ma.masked_equal(finalcutpolygon, gen.fillValueF)
1599
1600    Npts = finalcutpolygon.shape[0]
1601
1602    return Npts, finalcutpolygon
1603
1604def pile_polygons(polyns, polygons):
1605    """ Function to pile polygons one over the following one
1606      polyns: ordered list of polygons. First over all. last below all
1607      polygons: dictionary with the polygons
1608    >>> pns = ['sqra', 'sqrb']
1609    >>> polya = np.array([[-0.5, -0.75], [0.5, -0.75], [0.5, 0.75], [-0.5, 0.75]])
1610    >>> polyb = np.array([[-0.75, -0.5], [0.75, -0.5], [0.75, 0.5], [-0.75, 0.5]])
1611    >>> plgs = {'sqra': polya, 'sqrb': polyb}
1612    >>> pile_polygons(pns, plgs)
1613    #  sqrb :
1614    [[-0.75 -0.5]
1615     [-0.5 -0.5]
1616     [-- --]
1617     [0.5 -0.5]
1618     [0.75 -0.5]
1619     [0.75 0.5]
1620     [0.5 0.5]
1621     [-- --]
1622     [-0.5 0.5]
1623     [-0.75 0.5]
1624     [-0.75 -0.5]]
1625    #  sqra :
1626     [[-0.5  -0.75]
1627      [ 0.5  -0.75]
1628      [ 0.5   0.75]
1629      [-0.5   0.75]
1630      [-0.5  -0.75]]
1631    """
1632    fname = 'pile_polygons'
1633    pilepolygons = dict(polygons)
1634    Npolys = len(polyns)
1635
1636    for ipolyp in range(Npolys-2,-1,-1):
1637        polyn = polyns[ipolyp]
1638        poly = pilepolygons[polyn]
1639        Npts = poly.shape[0]
1640        for ipolyi in range(ipolyp+1,Npolys,1):
1641            ipolyn = polyns[ipolyi]
1642            #print '  Lluis ' + polyn + ' above ' + ipolyn
1643            ipoly = pilepolygons[ipolyn]
1644            iNpts = ipoly.shape[0]
1645            newipoly = []
1646
1647            Nint, inti, intp, pts = fsci.module_scientific.crossingpoints_polys(     \
1648              nvertexa=iNpts, nvertexb=Npts, nvertexab=iNpts*Npts, polya=ipoly,      \
1649              polyb=poly)
1650            # We're in C-mode !
1651            inti = inti-1
1652            intp = intp-1
1653
1654            # Re-constructing below polygon looking for respective crossings
1655            linti = list(inti)
1656            for ip in range(iNpts):
1657                iip1 = ip+1
1658                if ip == iNpts-1: iip1 = 0
1659                #print ip, ipoly[ip,:], ':', ipoly[iip1,:]
1660                Nc = linti.count(ip)
1661                ldists = []
1662                ddists = {}
1663                if Nc > 0:
1664                    iic = gen.multi_index_vec(inti,ip)
1665                    mindist = 1000000.
1666                    # Sorting from distance respect the vertex ip
1667                    for ic in range(Nc):
1668                        ddists[iic[ic]] = dist_points(ipoly[ip,:], pts[iic[ic],:])
1669                        ldists.append(ddists[iic[ic]])
1670                        #print '  ', ic, ';', iic[ic], '=', pts[iic[ic],:], ddists[iic[ic]]
1671                    ldists.sort()
1672                    #print '    ldists', ldists
1673                    newipoly.append(ipoly[ip,:])
1674                    for ic in range(Nc):
1675                        iic = gen.dictionary_key(ddists,ldists[ic])
1676                        newipoly.append(pts[iic,:])
1677                        #print '  ', ic, '|', iic, ';', pts[iic,:]
1678                        if ic < Nc-1: newipoly.append([gen.fillValueF, gen.fillValueF])
1679                    newipoly.append(ipoly[iip1,:])
1680
1681                else:
1682                    newipoly.append(ipoly[ip,:])
1683   
1684            newipoly = np.array(newipoly)
1685            pilepolygons[polyns[Npolys-1]] = rm_consecpt_polygon(newipoly)
1686   
1687        for polyn in polyns:
1688            poly = pilepolygons[polyn]
1689            poly = ma.masked_equal(poly, gen.fillValueF)
1690            pilepolygons[polyn] = poly
1691
1692    return pilepolygons
1693
1694####### ###### ##### #### ### ## #
1695# Shapes/objects
1696
1697def surface_sphere(radii,Npts):
1698    """ Function to provide an sphere as matrix of x,y,z coordinates
1699      radii: radii of the sphere
1700      Npts: number of points to discretisize longitues (half for latitudes)
1701    """
1702    fname = 'surface_sphere'
1703
1704    sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float)
1705    spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float)
1706    for ia in range(Npts):
1707        alpha = ia*2*np.pi/(Npts-1)
1708        for ib in range(Npts/2):
1709            beta = ib*np.pi/(2.*(Npts/2-1))
1710            sphereup[:,ib,ia] = position_sphere(radii, alpha, beta)
1711        for ib in range(Npts/2):
1712            beta = -ib*np.pi/(2.*(Npts/2-1))
1713            spheredown[:,ib,ia] = position_sphere(radii, alpha, beta)
1714
1715    return sphereup, spheredown
1716
1717def ellipse_polar(c, a, b, Nang=100):
1718    """ Function to determine an ellipse from its center and polar coordinates
1719        FROM: https://en.wikipedia.org/wiki/Ellipse
1720      c= coordinates of the center
1721      a= distance major axis
1722      b= distance minor axis
1723      Nang= number of angles to use
1724    """
1725    fname = 'ellipse_polar'
1726
1727    if np.mod(Nang,2) == 0: Nang=Nang+1
1728 
1729    dtheta = 2*np.pi/(Nang-1)
1730
1731    ellipse = np.zeros((Nang,2), dtype=np.float)
1732    for ia in range(Nang):
1733        theta = dtheta*ia
1734        rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 )
1735        x = rad*np.cos(theta)
1736        y = rad*np.sin(theta)
1737        ellipse[ia,:] = [y+c[0],x+c[1]]
1738
1739    return ellipse
1740
1741def hyperbola_polar(a, b, Nang=100):
1742    """ Fcuntion to determine an hyperbola in polar coordinates
1743        FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates
1744          x^2/a^2 - y^2/b^2 = 1
1745      a= x-parameter
1746      y= y-parameter
1747      Nang= number of angles to use
1748      DOES NOT WORK!!!!
1749    """
1750    fname = 'hyperbola_polar'
1751
1752    dtheta = 2.*np.pi/(Nang-1)
1753
1754    # Positive branch
1755    hyperbola_p = np.zeros((Nang,2), dtype=np.float)
1756    for ia in range(Nang):
1757        theta = dtheta*ia
1758        x = a*np.cosh(theta)
1759        y = b*np.sinh(theta)
1760        hyperbola_p[ia,:] = [y,x]
1761
1762    # Negative branch
1763    hyperbola_n = np.zeros((Nang,2), dtype=np.float)
1764    for ia in range(Nang):
1765        theta = dtheta*ia
1766        x = -a*np.cosh(theta)
1767        y = b*np.sinh(theta)
1768        hyperbola_n[ia,:] = [y,x]
1769
1770    return hyperbola_p, hyperbola_n
1771
1772def circ_sec(ptA, ptB, radii, arc='short', pos='left', Nang=100):
1773    """ Function union of point A and B by a section of a circle
1774      ptA= coordinates od the point A [yA, xA]
1775      ptB= coordinates od the point B [yB, xB]
1776      radii= radi of the circle to use to unite the points
1777      arc= which arc to be used ('short', default)
1778        'short': shortest angle between points
1779        'long': largest angle between points
1780      pos= orientation of the arc following clockwise union of points ('left', default)
1781        'left': to the left of union
1782        'right': to the right of union
1783      Nang= amount of angles to use
1784    """
1785    fname = 'circ_sec'
1786    availarc = ['short', 'long']
1787    availpos = ['left', 'right']
1788
1789    distAB = dist_points(ptA,ptB)
1790
1791    if distAB > radii:
1792        print errormsg
1793        print '  ' + fname + ': radii=', radii, " too small for the distance " +     \
1794          "between points !!"
1795        print '    distance between points:', distAB
1796        quit(-1)
1797
1798    # Coordinate increments
1799    dAB = np.abs(ptA-ptB)
1800
1801    # angle of the circular section joining points
1802    alpha = 2.*np.arcsin((distAB/2.)/radii)
1803
1804    # center along coincident bisection of the union
1805    xcc = -radii
1806    ycc = 0.
1807
1808    # Getting the arc of the circle at the x-axis
1809    if arc == 'short':
1810        dalpha = alpha/(Nang-1)
1811    elif arc == 'long':
1812        dalpha = (2.*np.pi - alpha)/(Nang-1)
1813    else:
1814        print errormsg
1815        print '  ' + fname + ": arc '" + arc + "' not ready !!" 
1816        print '    available ones:', availarc
1817        quit(-1)
1818    if pos == 'left': sign=-1.
1819    elif pos == 'right': sign=1.
1820    else:
1821        print errormsg
1822        print '  ' + fname + ": position '" + pos + "' not ready !!" 
1823        print '     available ones:', availpos
1824        quit(-1)
1825
1826    circ_sec = np.zeros((Nang,2), dtype=np.float)
1827    for ia in range(Nang):
1828        alpha = sign*dalpha*ia
1829        x = radii*np.cos(alpha)
1830        y = radii*np.sin(alpha)
1831
1832        circ_sec[ia,:] = [y+ycc,x+xcc]
1833
1834    # Angle of the points
1835    theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
1836
1837    # rotating angle of the circ
1838    if pos == 'left': 
1839        rotangle = theta + np.pi/2. - alpha/2.
1840    elif pos == 'right':
1841        rotangle = theta + 3.*np.pi/2. - alpha/2.
1842    else:
1843        print errormsg
1844        print '  ' + fname + ": position '" + pos + "' not ready !!" 
1845        print '     available ones:', availpos
1846        quit(-1)
1847
1848    #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi
1849 
1850    # rotating the arc along the x-axis
1851    rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle)
1852
1853    # Moving arc to the ptA
1854    circ_sec = rotcirc_sec + ptA
1855
1856    return circ_sec
1857
1858def p_square(face, N=5):
1859    """ Function to get a polygon square
1860      face: length of the face of the square
1861      N: number of points of the polygon
1862    """
1863    fname = 'p_square'
1864
1865    square = np.zeros((N,2), dtype=np.float)
1866
1867    f2 = face/2.
1868    N4 = N/4
1869    df = face/(N4)
1870    # SW-NW
1871    for ip in range(N4):
1872        square[ip,:] = [-f2+ip*df,-f2]
1873    # NW-NE
1874    for ip in range(N4):
1875        square[ip+N4,:] = [f2,-f2+ip*df]
1876    # NE-SE
1877    for ip in range(N4):
1878        square[ip+2*N4,:] = [f2-ip*df,f2]
1879    N42 = N-3*N4-1
1880    df = face/(N42)
1881    # SE-SW
1882    for ip in range(N42):
1883        square[ip+3*N4,:] = [-f2,f2-ip*df]
1884    square[N-1,:] = [-f2,-f2]
1885
1886    return square
1887
1888def p_prism(base, height, N=5):
1889    """ Function to get a polygon prism
1890      base: length of the base of the prism
1891      height: length of the height of the prism
1892      N: number of points of the polygon
1893    """
1894    fname = 'p_prism'
1895
1896    prism = np.zeros((N,2), dtype=np.float)
1897
1898    b2 = base/2.
1899    h2 = height/2.
1900    N4 = N/4
1901    dh = height/(N4)
1902    db = base/(N4)
1903
1904    # SW-NW
1905    for ip in range(N4):
1906        prism[ip,:] = [-h2+ip*dh,-b2]
1907    # NW-NE
1908    for ip in range(N4):
1909        prism[ip+N4,:] = [h2,-b2+ip*db]
1910    # NE-SE
1911    for ip in range(N4):
1912        prism[ip+2*N4,:] = [h2-ip*dh,b2]
1913    N42 = N-3*N4-1
1914    db = base/(N42)
1915    # SE-SW
1916    for ip in range(N42):
1917        prism[ip+3*N4,:] = [-h2,b2-ip*db]
1918    prism[N-1,:] = [-h2,-b2]
1919
1920    return prism
1921
1922def p_circle(radii, N=50):
1923    """ Function to get a polygon of a circle
1924      radii: length of the radii of the circle
1925      N: number of points of the polygon
1926    """
1927    fname = 'p_circle'
1928
1929    circle = np.zeros((N,2), dtype=np.float)
1930
1931    dangle = 2.*np.pi/(N-1)
1932
1933    for ia in range(N):
1934        circle[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
1935
1936    circle[N-1,:] = [0., radii]
1937
1938    return circle
1939
1940def p_triangle(p1, p2, p3, N=4):
1941    """ Function to provide the polygon of a triangle from its 3 vertices
1942      p1: vertex 1 [y,x]
1943      p2: vertex 2 [y,x]
1944      p3: vertex 3 [y,x]
1945      N: number of vertices of the triangle
1946    """
1947    fname = 'p_triangle'
1948
1949    triangle = np.zeros((N,2), dtype=np.float)
1950
1951    N3 = N / 3
1952    # 1-2
1953    dx = (p2[1]-p1[1])/N3
1954    dy = (p2[0]-p1[0])/N3
1955    for ip in range(N3):
1956        triangle[ip,:] = [p1[0]+ip*dy,p1[1]+ip*dx]
1957    # 2-3
1958    dx = (p3[1]-p2[1])/N3
1959    dy = (p3[0]-p2[0])/N3
1960    for ip in range(N3):
1961        triangle[ip+N3,:] = [p2[0]+ip*dy,p2[1]+ip*dx]
1962    # 3-1
1963    N32 = N - 2*N/3
1964    dx = (p1[1]-p3[1])/N32
1965    dy = (p1[0]-p3[0])/N32
1966    for ip in range(N32):
1967        triangle[ip+2*N3,:] = [p3[0]+ip*dy,p3[1]+ip*dx]
1968
1969    triangle[N-1,:] = p1
1970
1971    return triangle
1972
1973def p_spiral(loops, eradii, N=1000):
1974    """ Function to provide a polygon of an Archimedean spiral
1975        FROM: https://en.wikipedia.org/wiki/Spiral
1976      loops: number of loops of the spiral
1977      eradii: length of the radii of the final spiral
1978      N: number of points of the polygon
1979    """
1980    fname = 'p_spiral'
1981
1982    spiral = np.zeros((N,2), dtype=np.float)
1983
1984    dangle = 2.*np.pi*loops/(N-1)
1985    dr = eradii*1./(N-1)
1986
1987    for ia in range(N):
1988        radii = dr*ia
1989        spiral[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
1990
1991    return spiral
1992
1993def p_reg_polygon(Nv, lf, N=50):
1994    """ Function to provide a regular polygon of Nv vertices
1995      Nv: number of vertices
1996      lf: length of the face
1997      N: number of points
1998    """
1999    fname = 'p_reg_polygon'
2000
2001    reg_polygon = np.zeros((N,2), dtype=np.float)
2002
2003    # Number of points per vertex
2004    Np = N/Nv
2005    # Angle incremental between vertices
2006    da = 2.*np.pi/Nv
2007    # Radii of the circle according to lf
2008    radii = lf*Nv/(2*np.pi)
2009
2010    iip = 0
2011    for iv in range(Nv-1):
2012        # Characteristics between vertices iv and iv+1
2013        av1 = da*iv
2014        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
2015        av2 = da*(iv+1)
2016        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
2017        dx = (v2[1]-v1[1])/Np
2018        dy = (v2[0]-v1[0])/Np
2019        for ip in range(Np):
2020            reg_polygon[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
2021
2022    # Characteristics between vertices Nv and 1
2023
2024    # Number of points per vertex
2025    Np2 = N - Np*(Nv-1)
2026
2027    av1 = da*Nv
2028    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
2029    av2 = 0.
2030    v2 = [radii*np.sin(av2), radii*np.cos(av2)]
2031    dx = (v2[1]-v1[1])/Np2
2032    dy = (v2[0]-v1[0])/Np2
2033    for ip in range(Np2):
2034        reg_polygon[ip+(Nv-1)*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
2035
2036    return reg_polygon
2037
2038def p_reg_star(Nv, lf, freq, vs=0, N=50):
2039    """ Function to provide a regular star of Nv vertices
2040      Nv: number of vertices
2041      lf: length of the face of the regular polygon
2042      freq: frequency of union of vertices ('0', for just centered to zero arms)
2043      vs: vertex from which start (0 being first [0,lf])
2044      N: number of points
2045    """
2046    fname = 'p_reg_star'
2047
2048    reg_star = np.zeros((N,2), dtype=np.float)
2049
2050    # Number of arms of the star
2051    if freq != 0 and np.mod(Nv,freq) == 0: 
2052        Na = Nv/freq + 1
2053    else:
2054        Na = Nv
2055
2056    # Number of points per arm
2057    Np = N/Na
2058    # Angle incremental between vertices
2059    da = 2.*np.pi/Nv
2060    # Radii of the circle according to lf
2061    radii = lf*Nv/(2*np.pi)
2062
2063    iip = 0
2064    av1 = vs*da
2065    for iv in range(Na-1):
2066        # Characteristics between vertices iv and iv+1
2067        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
2068        if freq != 0:
2069            av2 = av1 + da*freq
2070            v2 = [radii*np.sin(av2), radii*np.cos(av2)]
2071        else:
2072            v2 = [0., 0.]
2073            av2 = av1 + da
2074        dx = (v2[1]-v1[1])/(Np-1)
2075        dy = (v2[0]-v1[0])/(Np-1)
2076        for ip in range(Np):
2077            reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
2078        if av2 > 2.*np.pi: av1 = av2 - 2.*np.pi
2079        else: av1 = av2 + 0.
2080
2081    iv = Na-1
2082    # Characteristics between vertices Na and 1
2083    Np2 = N-Np*iv
2084    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
2085    if freq != 0:
2086        av2 = vs*da
2087        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
2088    else:
2089        v2 = [0., 0.]
2090    dx = (v2[1]-v1[1])/(Np2-1)
2091    dy = (v2[0]-v1[0])/(Np2-1)
2092    for ip in range(Np2):
2093        reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
2094
2095    return reg_star
2096
2097def p_sinusoide(length=10., amp=5., lamb=3., ival=0., func='sin', N=100):
2098    """ Function to get coordinates of a sinusoidal curve
2099      length: length of the line (default 10.)
2100      amp: amplitude of the peaks (default 5.)
2101      lamb: wave longitude (defalult 3.)
2102      ival: initial angle (default 0. in degree)
2103      func: function to use: (default sinus)
2104        'sin': sinus
2105        'cos': cosinus
2106      N: number of points (default 100)
2107    """
2108    fname = 'p_sinusoide'
2109    availfunc = ['sin', 'cos']
2110
2111    dx = length/(N-1)
2112    ia = ival*np.pi/180.
2113    da = 2*np.pi*dx/lamb
2114
2115    sinusoide = np.zeros((N,2), dtype=np.float)
2116    if func == 'sin':
2117        for ix in range(N):
2118            sinusoide[ix,:] = [amp*np.sin(ia+da*ix),dx*ix]
2119    elif func == 'cos':
2120        for ix in range(N):
2121            sinusoide[ix,:] = [amp*np.cos(ia+da*ix),dx*ix]
2122    else:
2123        print errormsg
2124        print '  ' + fname + ": function '" + func + "' not ready !!"
2125        print '    available ones:', availfunc
2126        quit(-1)
2127
2128    sinusoidesecs = ['sinusoide']
2129    sinusoidedic = {'sinusoide': [sinusoide, '-', '#000000', 1.]}
2130
2131    return sinusoide, sinusoidesecs, sinusoidedic
2132
2133def p_doubleArrow(length=5., angle=45., width=1., alength=0.10, N=50):
2134    """ Function to provide an arrow with double lines
2135      length: length of the arrow (5. default)
2136      angle: angle of the head of the arrow (45., default)
2137      width: separation between the two lines (2., default)
2138      alength: length of the head (as percentage in excess of width, 0.1 default)
2139      N: number of points (50, default)
2140    """
2141    function = 'p_doubleArrow'
2142
2143    doubleArrow = np.zeros((50,2), dtype=np.float)
2144    N4 = int((N-3)/4)
2145
2146    doublearrowdic = {}
2147    ddy = width*np.tan(angle*np.pi/180.)/2.
2148    # Arms
2149    dx = (length-ddy)/(N4-1)
2150    for ix in range(N4):
2151        doubleArrow[ix,:] = [dx*ix,-width/2.]
2152    doublearrowdic['leftarm'] = [doubleArrow[0:N4,:], '-', '#000000', 2.]
2153    doubleArrow[N4,:] = [gen.fillValueF,gen.fillValueF]
2154    for ix in range(N4):
2155        doubleArrow[N4+1+ix,:] = [dx*ix,width/2.]
2156    doublearrowdic['rightarm'] = [doubleArrow[N4+1:2*N4+1,:], '-', '#000000', 2.]
2157    doubleArrow[2*N4+1,:] = [gen.fillValueF,gen.fillValueF]
2158
2159    # Head
2160    N42 = int((N-2 - 2*N4)/2)
2161    dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N42-1)
2162    dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N42-1)
2163    for ix in range(N42):
2164        doubleArrow[2*N4+2+ix,:] = [length-dy*ix,-dx*ix]
2165    doublearrowdic['lefthead'] = [doubleArrow[2*N4:2*N4+N42,:], '-', '#000000', 2.]
2166    doubleArrow[2*N4+2+N42,:] = [gen.fillValueF,gen.fillValueF]
2167
2168    N43 = N-3 - 2*N4 - N42 + 1
2169    dx = width*(1.+alength)*np.cos(angle*np.pi/180.)/(N43-1)
2170    dy = width*(1.+alength)*np.sin(angle*np.pi/180.)/(N43-1)
2171    for ix in range(N43):
2172        doubleArrow[2*N4+N42+2+ix,:] = [length-dy*ix,dx*ix]
2173    doublearrowdic['rightthead'] = [doubleArrow[2*N4+N42+2:51,:], '-', '#000000', 2.]
2174
2175    doubleArrow = ma.masked_equal(doubleArrow, gen.fillValueF)
2176    doublearrowsecs = ['leftarm', 'rightarm', 'lefthead', 'righthead']
2177
2178    return doubleArrow, doublearrowsecs, doublearrowdic
2179
2180def p_angle_triangle(pi=np.array([0.,0.]), angle1=60., length1=1., angle2=60., N=100):
2181    """ Function to draw a triangle by an initial point and two consecutive angles
2182        and the first length of face. The third angle and 2 and 3rd face will be
2183        computed accordingly the provided values:
2184           length1 / sin(angle1) = length2 / sin(angle2) = length3 / sin(angle3)
2185           angle1 + angle2 + angle3 = 180.
2186      pi: initial point ([0., 0.], default)
2187      angle1: first angle from pi clockwise (60., default)
2188      length1: length of face from pi by angle1 (1., default)
2189      angle2: second angle from second point (60., default)
2190      length2: length of face from p2 by angle2 (1., default)
2191      N: number of points (100, default)
2192    """
2193    fname = 'p_angle_triangle'
2194
2195    angle3 = 180. - angle1 - angle2
2196    length2 = np.sin(angle2*np.pi/180.)*length1/np.sin(angle1*np.pi/180.)
2197    length3 = np.sin(angle3*np.pi/180.)*length1/np.sin(angle1*np.pi/180.)
2198
2199    triangle = np.zeros((N,2), dtype=np.float)
2200
2201    N3 = int(N/3)
2202    # first face
2203    ix = pi[1]
2204    iy = pi[0]
2205    dx = length1*np.cos(angle1*np.pi/180.)/(N3-1)
2206    dy = length1*np.sin(angle1*np.pi/180.)/(N3-1)
2207    for ip in range(N3):
2208        triangle[ip,:] = [iy+dy*ip, ix+dx*ip]
2209
2210    # second face
2211    ia = -90. - (90.-angle1)
2212    ix = triangle[N3-1,1]
2213    iy = triangle[N3-1,0]
2214    dx = length2*np.cos((ia+angle2)*np.pi/180.)/(N3-1)
2215    dy = length2*np.sin((ia+angle2)*np.pi/180.)/(N3-1)
2216    for ip in range(N3):
2217        triangle[N3+ip,:] = [iy+dy*ip, ix+dx*ip]
2218
2219    # third face
2220    N32 = N - 2*N3
2221    ia = -180. - (90.-angle2)
2222    ix = triangle[2*N3-1,1]
2223    iy = triangle[2*N3-1,0]
2224    angle3 = np.arctan2(pi[0]-iy, pi[1]-ix)
2225    dx = (pi[1]-ix)/(N32-1)
2226    dy = (pi[0]-iy)/(N32-1)
2227    for ip in range(N32):
2228        triangle[2*N3+ip,:] = [iy+dy*ip, ix+dx*ip]
2229
2230    return triangle
2231
2232def p_cross_width(larm=5., width=1., Narms=4, N=200):
2233    """ Function to draw a cross with arms with a given width and an angle
2234      larm: legnth of the arms (5., default)
2235      width: width of the arms (1., default)
2236      Narms: Number of arms (4, default)
2237      N: number of points to us (200, default)
2238    """
2239    fname = 'p_cross_width'
2240
2241    Narm = int((N-Narms)/Narms)
2242
2243    larm2 = larm/2.
2244    width2 = width/2.
2245
2246    cross = np.ones((N,2), dtype=np.float)*gen.fillValueF
2247    da = np.pi/Narms
2248
2249    N1 = int(Narm*3./8.)
2250    N2 = int((Narm - 2*N1)/2.)
2251    N21 = Narm - 2*N1 - N2
2252
2253    if N2 < 3:
2254        print errormsg
2255        print '  ' + fname + ": too few points for ", Narms, " arms !!"
2256        print "    increase number 'N' at least up to '", 25*Narms
2257        quit(-1)
2258
2259    crosssecs = []
2260    crossdic = {}
2261    Npot = int(np.log10(Narms))+1
2262
2263    iip = 0
2264    for iarm in range(Narms-1):
2265
2266        a = da*iarm
2267        iip0 = iip
2268
2269        # bottom coordinate
2270        bx = larm*np.cos(a+np.pi)
2271        by = larm*np.sin(a+np.pi)
2272
2273        # upper coordinate
2274        ux = larm*np.cos(a)
2275        uy = larm*np.sin(a)
2276
2277        rela = a+np.pi*3./2.
2278        # SW-NW
2279        ix = bx + width2*np.cos(rela)
2280        iy = by + width2*np.sin(rela)
2281        ex = ux + width2*np.cos(rela)
2282        ey = uy + width2*np.sin(rela)
2283        dx = (ex-ix)/(N1-1)
2284        dy = (ey-iy)/(N1-1)
2285        for ip in range(N1):
2286            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2287        iip = iip + N1
2288
2289        # NW-NE
2290        ix = ex + 0.
2291        iy = ey + 0.
2292        ex = ux - width2*np.cos(rela)
2293        ey = uy - width2*np.sin(rela)
2294        dx = (ex-ix)/(N2-1)
2295        dy = (ey-iy)/(N2-1)
2296        for ip in range(N2):
2297            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2298        iip = iip + N2
2299
2300        # NW-SW
2301        ix = ex + 0.
2302        iy = ey + 0.
2303        ex = bx - width2*np.cos(rela)
2304        ey = by - width2*np.sin(rela)
2305        dx = (ex-ix)/(N1-1)
2306        dy = (ey-iy)/(N1-1)
2307        for ip in range(N1):
2308            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2309        iip = iip + N1
2310
2311        # SW-SE
2312        ix = ex + 0.
2313        iy = ey + 0.
2314        ex = bx + width2*np.cos(rela)
2315        ey = by + width2*np.sin(rela)
2316        dx = (ex-ix)/(N21-1)
2317        dy = (ey-iy)/(N21-1)
2318        for ip in range(N21):
2319            cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2320        iip = iip + N21 + 1
2321
2322        iarmS = str(iarm).zfill(Npot)
2323        crosssecs.append(iarmS)
2324        crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.']
2325
2326    iip0 = iip
2327
2328    Narm = N - Narm*(Narms-1) - Narms
2329
2330    N1 = int(Narm*3./8.)
2331    N2 = int((Narm - 2*N1)/2.)
2332    N21 = Narm - 2*N1 - N2
2333
2334    iarm = Narms-1
2335    a = da*iarm
2336
2337    # bottom coordinate
2338    bx = larm*np.cos(a+np.pi)
2339    by = larm*np.sin(a+np.pi)
2340
2341    # upper coordinate
2342    ux = larm*np.cos(a)
2343    uy = larm*np.sin(a)
2344
2345    rela = a+np.pi*3./2.
2346    # SW-NW
2347    ix = bx + width2*np.cos(rela)
2348    iy = by + width2*np.sin(rela)
2349    ex = ux + width2*np.cos(rela)
2350    ey = uy + width2*np.sin(rela)
2351    dx = (ex-ix)/(N1-1)
2352    dy = (ey-iy)/(N1-1)
2353    for ip in range(N1):
2354      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2355    iip = iip + N1
2356
2357    # NW-NE
2358    ix = ex + 0.
2359    iy = ey + 0.
2360    ex = ux - width2*np.cos(rela)
2361    ey = uy - width2*np.sin(rela)
2362    dx = (ex-ix)/(N2-1)
2363    dy = (ey-iy)/(N2-1)
2364    for ip in range(N2):
2365      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2366    iip = iip + N2
2367
2368    # NW-SW
2369    ix = ex + 0.
2370    iy = ey + 0.
2371    ex = bx - width2*np.cos(rela)
2372    ey = by - width2*np.sin(rela)
2373    dx = (ex-ix)/(N1-1)
2374    dy = (ey-iy)/(N1-1)
2375    for ip in range(N1):
2376      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2377    iip = iip + N1
2378
2379    # SW-SE
2380    ix = ex + 0.
2381    iy = ey + 0.
2382    ex = bx + width2*np.cos(rela)
2383    ey = by + width2*np.sin(rela)
2384    dx = (ex-ix)/(N21-1)
2385    dy = (ey-iy)/(N21-1)
2386    for ip in range(N21):
2387      cross[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
2388    iip = iip + N21
2389
2390    iarmS = str(iarm).zfill(Npot)
2391    crosssecs.append(iarmS)
2392    crossdic[iarmS] = [cross[iip0:iip0+iip-1], '-', 'k', '1.']
2393
2394    cross = ma.masked_equal(cross, gen.fillValueF)
2395
2396    return cross, crosssecs, crossdic
2397
2398####### ####### ##### #### ### ## #
2399# Plotting
2400
2401def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10,                   \
2402  drwsfc=[True,True], colsfc=['#AAAAAA','#646464'],                                  \
2403  drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.],          \
2404  drwzline = True, linez=['-.','g',2.], drwxcline=[True,True],                       \
2405  linexc=[['-','#646400',1.],['--','#646400',1.]],                                   \
2406  drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]],           \
2407  drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]):
2408    """ Function to plot an sphere and determine which standard lines will be also
2409        drawn
2410      iazm: azimut of the camera form the sphere
2411      iele: elevation of the camera form the sphere
2412      dist: distance of the camera form the sphere
2413      Npts: Resolution for the sphere
2414      radii: radius of the sphere
2415      drwsfc: whether 'up' and 'down' portions of the sphere should be drawn
2416      colsfc: colors of the surface of the sphere portions ['up', 'down']
2417      drwxline: whether x-axis line should be drawn
2418      linex: properties of the x-axis line ['type', 'color', 'wdith']
2419      drwyline: whether y-axis line should be drawn
2420      liney: properties of the y-axis line ['type', 'color', 'wdith']
2421      drwzline: whether z-axis line should be drawn
2422      linez: properties of the z-axis line ['type', 'color', 'wdith']
2423      drwequator: whether 'front' and 'back' portions of the Equator should be drawn
2424      lineeq: properties of the lines 'front' and 'back' of the Equator
2425      drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn
2426      linegw: properties of the lines 'front' and 'back' Greenwhich
2427      drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn
2428      linexc: properties of the lines 'front' and 'back' for the 90 line
2429    """
2430    fname = 'plot_sphere'
2431
2432    iazmrad = iazm*np.pi/180.
2433    ielerad = iele*np.pi/180.
2434
2435    # 3D surface Sphere
2436    sfcsphereu, sfcsphered = surface_sphere(radii,Npts)
2437   
2438    # greenwhich
2439    if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.:
2440        ia=np.pi-ielerad
2441    else:
2442        ia=0.-ielerad
2443    ea=ia+np.pi
2444    da = (ea-ia)/(Npts-1)
2445    beta = np.arange(ia,ea+da,da)[0:Npts]
2446    alpha = np.zeros((Npts), dtype=np.float)
2447    greenwhichc = spheric_line(radii,alpha,beta)
2448    ia=ea+0.
2449    ea=ia+np.pi
2450    da = (ea-ia)/(Npts-1)
2451    beta = np.arange(ia,ea+da,da)[0:Npts]
2452    greenwhichd = spheric_line(radii,alpha,beta)
2453
2454    # Equator
2455    ia=np.pi-iazmrad/2.
2456    ea=ia+np.pi
2457    da = (ea-ia)/(Npts-1)
2458    alpha = np.arange(ia,ea+da,da)[0:Npts]
2459    beta = np.zeros((Npts), dtype=np.float)
2460    equatorc = spheric_line(radii,alpha,beta)
2461    ia=ea+0.
2462    ea=ia+np.pi
2463    da = (ea-ia)/(Npts-1)
2464    alpha = np.arange(ia,ea+da,da)[0:Npts]
2465    equatord = spheric_line(radii,alpha,beta)
2466
2467    # 90 line
2468    if iazmrad > np.pi and iazmrad < 2.*np.pi:
2469        ia=3.*np.pi/2. + ielerad
2470    else:
2471        ia=np.pi/2. - ielerad
2472    if ielerad < 0.:
2473        ia = ia + np.pi
2474    ea=ia+np.pi
2475    da = (ea-ia)/(Npts-1)
2476    beta = np.arange(ia,ea+da,da)[0:Npts]
2477    alpha = np.ones((Npts), dtype=np.float)*np.pi/2.
2478    xclinec = spheric_line(radii,alpha,beta)
2479    ia=ea+0.
2480    ea=ia+np.pi
2481    da = (ea-ia)/(Npts-1)
2482    beta = np.arange(ia,ea+da,da)[0:Npts]
2483    xclined = spheric_line(radii,alpha,beta)
2484
2485    # x line
2486    xline = np.zeros((2,3), dtype=np.float)
2487    xline[0,:] = position_sphere(radii, 0., 0.)
2488    xline[1,:] = position_sphere(radii, np.pi, 0.)
2489
2490    # y line
2491    yline = np.zeros((2,3), dtype=np.float)
2492    yline[0,:] = position_sphere(radii, np.pi/2., 0.)
2493    yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.)
2494
2495    # z line
2496    zline = np.zeros((2,3), dtype=np.float)
2497    zline[0,:] = position_sphere(radii, 0., np.pi/2.)
2498    zline[1,:] = position_sphere(radii, 0., -np.pi/2.)
2499
2500    fig = plt.figure()
2501    ax = fig.gca(projection='3d')
2502
2503    # Sphere surface
2504    if drwsfc[0]:
2505        ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:],     \
2506          color=colsfc[0])
2507    if drwsfc[1]:
2508        ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:],     \
2509          color=colsfc[1])
2510
2511    # greenwhich
2512    linev = linegw[0]
2513    if drwgreeenwhich[0]:
2514        ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0],      \
2515          color=linev[1], linewidth=linev[2],  label='Greenwich')
2516    linev = linegw[1]
2517    if drwgreeenwhich[1]:
2518        ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0],      \
2519          color=linev[1], linewidth=linev[2])
2520
2521    # Equator
2522    linev = lineeq[0]
2523    if drwequator[0]:
2524        ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0],               \
2525          color=linev[1], linewidth=linev[2], label='Equator')
2526    linev = lineeq[1]
2527    if drwequator[1]:
2528        ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0],               \
2529          color=linev[1], linewidth=linev[2])
2530
2531    # 90line
2532    linev = linexc[0]
2533    if drwxcline[0]:
2534        ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1],  \
2535          linewidth=linev[2], label='90-line')
2536    linev = linexc[1]
2537    if drwxcline[1]:
2538        ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1],  \
2539          linewidth=linev[2])
2540
2541    # x line
2542    linev = linex
2543    if drwxline:
2544        ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]],                    \
2545          [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
2546          label='xline')
2547
2548    # y line
2549    linev = liney
2550    if drwyline:
2551        ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]],                    \
2552          [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
2553          label='yline')
2554
2555    # z line
2556    linev = linez
2557    if drwzline:
2558        ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]],                    \
2559          [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2],     \
2560          label='zline')
2561
2562    plt.legend()
2563
2564    return fig, ax
2565
2566def draw_secs(objdic):
2567    """ Function to draw an object according to its dictionary
2568      objdic: dictionary with the parts to draw [polygon, ltype, lcol, lw]
2569    """
2570    fname = 'draw_secs'
2571
2572    for secn in objdic.keys():
2573        secv = objdic[secn]
2574        poly = secv[0]
2575        lt = secv[1]
2576        lc = secv[2]
2577        lw = secv[3]
2578
2579        plt.plot(poly[:,1], poly[:,0], lt, color=lc, linewidth=lw)
2580
2581    return
2582
2583def paint_filled(objdic, fillsecs):
2584    """ Function to draw an object filling given sections
2585      objdic: dictionary of the object
2586      filesecs: list of sections to be filled
2587    """
2588    fname = 'paint_filled'
2589
2590    Nsecs = len(fillsecs)
2591
2592    for secn in fillsecs:
2593        secvals=objdic[secn]
2594        pvals = secvals[0]
2595        fillsecs = []
2596        Nvals = pvals.shape[0]
2597        # re-sectionning to plot without masked values
2598        for ip in range(Nvals-1):
2599            if type(pvals[ip][0]) == type(gen.mamat[1]): fillsecs.append(ip) 
2600        Nsecs = len(fillsecs)
2601        iisc = 0
2602        for isc in range(Nsecs):
2603            plt.fill(pvals[iisc:fillsecs[isc],1], pvals[iisc:fillsecs[isc],0],       \
2604              color=secvals[2])
2605            iisc = fillsecs[isc]+1
2606        plt.fill(pvals[iisc:Nvals-1,1], pvals[iisc:Nvals-1,0], color=secvals[2])
2607
2608    return
2609
Note: See TracBrowser for help on using the repository browser.