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

Last change on this file since 2640 was 2630, checked in by lfita, 6 years ago

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