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

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

Adding:

  • `join_circ_sec': Function to join aa series of points by circular segments
File size: 39.1 KB
Line 
1# Python tools to manage netCDF files.
2# L. Fita, CIMA. Mrch 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
18
19errormsg = 'ERROR -- error -- ERROR -- error'
20infmsg = 'INFORMATION -- information -- INFORMATION -- information'
21
22####### Contents:
23# deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
24# dist_points: Function to provide the distance between two points
25# join_circ_sec: Function to join aa series of points by circular segments
26# max_coords_poly: Function to provide the extremes of the coordinates of a polygon
27# mirror_polygon: Function to reflex a polygon for a given axis
28# position_sphere: Function to tranform fom a point in lon, lat deg coordinates to
29#   cartesian coordinates over an sphere
30# read_join_poly: Function to read an ASCII file with the combination of polygons
31# rotate_2D: Function to rotate a vector by a certain angle in the plain
32# rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon
33# rotate_line2D: Function to rotate a line given by 2 pairs of x,y coordinates by a
34#   certain angle in the plain
35# rotate_lines2D: Function to rotate multiple lines given by mulitple pars of x,y
36#   coordinates by a certain angle in the plain
37# spheric_line: Function to transform a series of locations in lon, lat coordinates
38#   to x,y,z over an 3D spaceFunction to provide coordinates of a line  on a 3D space
39# write_join_poly: Function to write an ASCII file with the combination of polygons
40
41## Shapes/objects
42# circ_sec: Function union of point A and B by a section of a circle
43# ellipse_polar: Function to determine an ellipse from its center and polar coordinates
44# p_doubleArrow: Function to provide an arrow with double lines
45# p_circle: Function to get a polygon of a circle
46# p_reg_polygon: Function to provide a regular polygon of Nv vertices
47# p_reg_star: Function to provide a regular star of Nv vertices
48# p_sinusiode: Function to get coordinates of a sinusoidal curve
49# p_square: Function to get a polygon square
50# p_spiral: Function to provide a polygon of an Archimedean spiral
51# p_triangle: Function to provide the polygon of a triangle from its 3 vertices
52# surface_sphere: Function to provide an sphere as matrix of x,y,z coordinates
53# z_boat: Function to define an schematic boat from the z-plane
54# zsailing_boat: Function to define an schematic sailing boat from the z-plane with sails
55# zisland1: Function to draw an island from z-axis as the union of a series of points by
56#   circular segments
57## Plotting
58# plot_sphere: Function to plot an sphere and determine which standard lines will be
59#   also drawn
60
61def deg_deci(angle):
62    """ Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
63      angle: list of [deg, minute, sec] to pass
64    >>> deg_deci([41., 58., 34.])
65    0.732621346072
66    """
67    fname = 'deg_deci'
68
69    deg = np.abs(angle[0]) + np.abs(angle[1])/60. + np.abs(angle[2])/3600.
70
71    if angle[0] < 0.: deg = -deg*np.pi/180.
72    else: deg = deg*np.pi/180.
73
74    return deg
75
76def position_sphere(radii, alpha, beta):
77    """ Function to tranform fom a point in lon, lat deg coordinates to cartesian 
78          coordinates over an sphere
79      radii: radii of the sphere
80      alpha: longitude of the point
81      beta: latitude of the point
82    >>> position_sphere(10., 30., 45.)
83    (0.81031678432964027, -5.1903473778327376, 8.5090352453411846
84    """
85    fname = 'position_sphere'
86
87    xpt = radii*np.cos(beta)*np.cos(alpha)
88    ypt = radii*np.cos(beta)*np.sin(alpha)
89    zpt = radii*np.sin(beta)
90
91    return xpt, ypt, zpt
92
93def spheric_line(radii,lon,lat):
94    """ Function to transform a series of locations in lon, lat coordinates to x,y,z
95          over an 3D space
96      radii: radius of the sphere
97      lon: array of angles along longitudes
98      lat: array of angles along latitudes
99    """
100    fname = 'spheric_line'
101
102    Lint = lon.shape[0]
103    coords = np.zeros((Lint,3), dtype=np.float)
104
105    for iv in range(Lint):
106        coords[iv,:] = position_sphere(radii, lon[iv], lat[iv])
107
108    return coords
109
110def rotate_2D(vector, angle):
111    """ Function to rotate a vector by a certain angle in the plain
112      vector= vector to rotate [y, x]
113      angle= angle to rotate [rad]
114    >>> rotate_2D(np.array([1.,0.]), np.pi/4.)
115    [ 0.70710678 -0.70710678]
116    """
117    fname = 'rotate_2D'
118
119    rotmat = np.zeros((2,2), dtype=np.float)
120
121    rotmat[0,0] = np.cos(angle)
122    rotmat[0,1] = -np.sin(angle)
123    rotmat[1,0] = np.sin(angle)
124    rotmat[1,1] = np.cos(angle)
125
126    rotvector = np.zeros((2), dtype=np.float)
127
128    vecv = np.zeros((2), dtype=np.float)
129
130    # Unifying vector
131    modvec = vector[0]**2+vector[1]**2
132    if modvec != 0: 
133        vecv[0] = vector[1]/modvec
134        vecv[1] = vector[0]/modvec
135
136        rotvec = np.matmul(rotmat, vecv)
137        rotvec = np.where(np.abs(rotvec) < 1.e-7, 0., rotvec)
138
139        rotvector[0] = rotvec[1]*modvec
140        rotvector[1] = rotvec[0]*modvec
141
142    return rotvector
143
144def rotate_polygon_2D(vectors, angle):
145    """ Function to rotate 2D plain the vertices of a polygon
146      line= matrix of vectors to rotate
147      angle= angle to rotate [rad]
148    >>> square = np.zeros((4,2), dtype=np.float)
149    >>> square[0,:] = [-0.5,-0.5]
150    >>> square[1,:] = [0.5,-0.5]
151    >>> square[2,:] = [0.5,0.5]
152    >>> square[3,:] = [-0.5,0.5]
153    >>> rotate_polygon_2D(square, np.pi/4.)
154    [[-0.70710678  0.        ]
155     [ 0.         -0.70710678]
156     [ 0.70710678  0.        ]
157     [ 0.          0.70710678]]
158    """
159    fname = 'rotate_polygon_2D'
160
161    rotvecs = np.zeros(vectors.shape, dtype=np.float)
162
163    Nvecs = vectors.shape[0]
164    for iv in range(Nvecs):
165        rotvecs[iv,:] = rotate_2D(vectors[iv,:], angle)
166
167    return rotvecs
168
169def rotate_line2D(line, angle):
170    """ Function to rotate a line given by 2 pairs of x,y coordinates by a certain
171          angle in the plain
172      line= line to rotate as couple of points [[y0,x0], [y1,x1]]
173      angle= angle to rotate [rad]
174    >>> rotate_line2D(np.array([[0.,0.], [1.,0.]]), np.pi/4.)
175    [[ 0.          0.        ]
176     [0.70710678  -0.70710678]]
177    """
178    fname = 'rotate_2D'
179
180    rotline = np.zeros((2,2), dtype=np.float)
181    rotline[0,:] = rotate_2D(line[0,:], angle)
182    rotline[1,:] = rotate_2D(line[1,:], angle)
183
184    return rotline
185
186def rotate_lines2D(lines, angle):
187    """ Function to rotate multiple lines given by mulitple pars of x,y coordinates 
188          by a certain angle in the plain
189      line= matrix of N couples of points [N, [y0,x0], [y1,x1]]
190      angle= angle to rotate [rad]
191    >>> square = np.zeros((4,2,2), dtype=np.float)
192    >>> square[0,0,:] = [-0.5,-0.5]
193    >>> square[0,1,:] = [0.5,-0.5]
194    >>> square[1,0,:] = [0.5,-0.5]
195    >>> square[1,1,:] = [0.5,0.5]
196    >>> square[2,0,:] = [0.5,0.5]
197    >>> square[2,1,:] = [-0.5,0.5]
198    >>> square[3,0,:] = [-0.5,0.5]
199    >>> square[3,1,:] = [-0.5,-0.5]
200    >>> rotate_lines2D(square, np.pi/4.)
201    [[[-0.70710678  0.        ]
202      [ 0.         -0.70710678]]
203
204     [[ 0.         -0.70710678]
205      [ 0.70710678  0.        ]]
206
207     [[ 0.70710678  0.        ]
208      [ 0.          0.70710678]]
209
210     [[ 0.          0.70710678]
211      [-0.70710678  0.        ]]]
212    """
213    fname = 'rotate_lines2D'
214
215    rotlines = np.zeros(lines.shape, dtype=np.float)
216
217    Nlines = lines.shape[0]
218    for il in range(Nlines):
219        line = np.zeros((2,2), dtype=np.float)
220        line[0,:] = lines[il,0,:]
221        line[1,:] = lines[il,1,:]
222
223        rotlines[il,:,:] = rotate_line2D(line, angle)
224
225    return rotlines
226
227def dist_points(ptA, ptB):
228    """ Function to provide the distance between two points
229      ptA: coordinates of the point A [yA, xA]
230      ptB: coordinates of the point B [yB, xB]
231    >>> dist_points([1.,1.], [-1.,-1.])
232    2.82842712475
233    """
234    fname = 'dist_points'
235
236    dist = np.sqrt( (ptA[0]-ptB[0])**2 + (ptA[1]-ptB[1])**2)
237
238    return dist
239
240def max_coords_poly(polygon):
241    """ Function to provide the extremes of the coordinates of a polygon
242      polygon: coordinates [Nvertexs, 2] of a polygon
243    >>> square = np.zeros((4,2), dtype=np.float)
244    >>> square[0,:] = [-0.5,-0.5]
245    >>> square[1,:] = [0.5,-0.5]
246    >>> square[2,:] = [0.5,0.5]
247    >>> square[3,:] = [-0.5,0.5]
248    >>> max_coords_poly(square)
249    [-0.5, 0.5], [-0.5, 0.5], [0.5, 0.5], 0.5
250    """
251    fname = 'max_coords_poly'
252
253    # x-coordinate min/max
254    nx = np.min(polygon[:,1])
255    xx = np.max(polygon[:,1])
256
257    # y-coordinate min/max
258    ny = np.min(polygon[:,0])
259    xy = np.max(polygon[:,0])
260
261    # x/y-coordinate maximum of absolute values
262    axx = np.max(np.abs(polygon[:,1]))
263    ayx = np.max(np.abs(polygon[:,0]))
264
265    # absolute maximum
266    xyx = np.max([axx, ayx])
267
268    return [nx, xx], [ny, xy], [ayx, axx], xyx
269
270def mirror_polygon(polygon,axis):
271    """ Function to reflex a polygon for a given axis
272      polygon: polygon to mirror
273      axis: axis at which mirror is located ('x' or 'y')
274    """
275    fname = 'mirror_polygon'
276
277    reflex = np.zeros(polygon.shape, dtype=np.float)
278
279    N = polygon.shape[0]
280    if axis == 'x':
281        for iv in range(N):
282            reflex[iv,:] = [-polygon[iv,0], polygon[iv,1]]
283    elif axis == 'y':
284        for iv in range(N):
285            reflex[iv,:] = [polygon[iv,0], -polygon[iv,1]]
286
287    return reflex
288
289def join_circ_sec(points, radfrac=3., N=200):
290    """ Function to join aa series of points by circular segments
291      points: main points of the island (clockwise ordered, to be joined by circular
292        segments of radii as the radfrac factor of the distance between
293        consecutive points)
294      radfrac: multiplicative factor of the distance between consecutive points to
295        draw the circular segment (3., default)
296      N: number of points (200, default)
297    """
298    fname = 'join_circ_sec'
299
300    jcirc_sec = np.ones((N,2), dtype=np.float)
301
302    # main points
303    lpoints = list(points)
304    Npts = len(lpoints)
305    Np = int(N/(Npts+1))
306    for ip in range(Npts-1):
307        dps = dist_points(lpoints[ip], lpoints[ip+1])
308        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(lpoints[ip], lpoints[ip+1],          \
309          dps*radfrac, Np)
310
311    Np2 = N - (Npts-1)*Np
312    dps = dist_points(lpoints[Npts-1], lpoints[0])
313    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(lpoints[Npts-1], lpoints[0], dps*3., Np2)
314
315    return jcirc_sec
316
317####### ###### ##### #### ### ## #
318# Shapes/objects
319
320def surface_sphere(radii,Npts):
321    """ Function to provide an sphere as matrix of x,y,z coordinates
322      radii: radii of the sphere
323      Npts: number of points to discretisize longitues (half for latitudes)
324    """
325    fname = 'surface_sphere'
326
327    sphereup = np.zeros((3,Npts/2,Npts), dtype=np.float)
328    spheredown = np.zeros((3,Npts/2,Npts), dtype=np.float)
329    for ia in range(Npts):
330        alpha = ia*2*np.pi/(Npts-1)
331        for ib in range(Npts/2):
332            beta = ib*np.pi/(2.*(Npts/2-1))
333            sphereup[:,ib,ia] = position_sphere(radii, alpha, beta)
334        for ib in range(Npts/2):
335            beta = -ib*np.pi/(2.*(Npts/2-1))
336            spheredown[:,ib,ia] = position_sphere(radii, alpha, beta)
337
338    return sphereup, spheredown
339
340def ellipse_polar(c, a, b, Nang=100):
341    """ Function to determine an ellipse from its center and polar coordinates
342        FROM: https://en.wikipedia.org/wiki/Ellipse
343      c= coordinates of the center
344      a= distance major axis
345      b= distance minor axis
346      Nang= number of angles to use
347    """
348    fname = 'ellipse_polar'
349
350    if np.mod(Nang,2) == 0: Nang=Nang+1
351 
352    dtheta = 2*np.pi/(Nang-1)
353
354    ellipse = np.zeros((Nang,2), dtype=np.float)
355    for ia in range(Nang):
356        theta = dtheta*ia
357        rad = a*b/np.sqrt( (b*np.cos(theta))**2 + (a*np.sin(theta))**2 )
358        x = rad*np.cos(theta)
359        y = rad*np.sin(theta)
360        ellipse[ia,:] = [y+c[0],x+c[1]]
361
362    return ellipse
363
364def hyperbola_polar(a, b, Nang=100):
365    """ Fcuntion to determine an hyperbola in polar coordinates
366        FROM: https://en.wikipedia.org/wiki/Hyperbola#Polar_coordinates
367          x^2/a^2 - y^2/b^2 = 1
368      a= x-parameter
369      y= y-parameter
370      Nang= number of angles to use
371      DOES NOT WORK!!!!
372    """
373    fname = 'hyperbola_polar'
374
375    dtheta = 2.*np.pi/(Nang-1)
376
377    # Positive branch
378    hyperbola_p = np.zeros((Nang,2), dtype=np.float)
379    for ia in range(Nang):
380        theta = dtheta*ia
381        x = a*np.cosh(theta)
382        y = b*np.sinh(theta)
383        hyperbola_p[ia,:] = [y,x]
384
385    # Negative branch
386    hyperbola_n = np.zeros((Nang,2), dtype=np.float)
387    for ia in range(Nang):
388        theta = dtheta*ia
389        x = -a*np.cosh(theta)
390        y = b*np.sinh(theta)
391        hyperbola_n[ia,:] = [y,x]
392
393    return hyperbola_p, hyperbola_n
394
395def circ_sec(ptA, ptB, radii, Nang=100):
396    """ Function union of point A and B by a section of a circle
397      ptA= coordinates od the point A [yA, xA]
398      ptB= coordinates od the point B [yB, xB]
399      radii= radi of the circle to use to unite the points
400      Nang= amount of angles to use
401    """
402    fname = 'circ_sec'
403
404    distAB = dist_points(ptA,ptB)
405
406    if distAB > radii:
407        print errormsg
408        print '  ' + fname + ': radii=', radii, " too small for the distance " +     \
409          "between points !!"
410        print '    distance between points:', distAB
411        quit(-1)
412
413    # Coordinate increments
414    dAB = np.abs(ptA-ptB)
415
416    # angle of the circular section joining points
417    alpha = 2.*np.arcsin((distAB/2.)/radii)
418
419    # center along coincident bisection of the union
420    xcc = -radii
421    ycc = 0.
422
423    # Getting the arc of the circle at the x-axis
424    dalpha = alpha/(Nang-1)
425    circ_sec = np.zeros((Nang,2), dtype=np.float)
426    for ia in range(Nang):
427        alpha = dalpha*ia
428        x = radii*np.cos(alpha)
429        y = radii*np.sin(alpha)
430
431        circ_sec[ia,:] = [y+ycc,x+xcc]
432   
433    # Angle of the points
434    theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
435
436    # rotating angle of the circ
437    rotangle = theta + 3.*np.pi/2. - alpha/2.
438
439    #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi
440 
441
442    # rotating the arc along the x-axis
443    rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle)
444
445    # Moving arc to the ptA
446    circ_sec = rotcirc_sec + ptA
447
448    return circ_sec
449
450def p_square(face, N=5):
451    """ Function to get a polygon square
452      face: length of the face of the square
453      N: number of points of the polygon
454    """
455    fname = 'p_square'
456
457    square = np.zeros((N,2), dtype=np.float)
458
459    f2 = face/2.
460    N4 = N/4
461    df = face/(N4)
462    # SW-NW
463    for ip in range(N4):
464        square[ip,:] = [-f2+ip*df,-f2]
465    # NW-NE
466    for ip in range(N4):
467        square[ip+N4,:] = [f2,-f2+ip*df]
468    # NE-SE
469    for ip in range(N4):
470        square[ip+2*N4,:] = [f2-ip*df,f2]
471    N42 = N-3*N4-1
472    df = face/(N42)
473    # SE-SW
474    for ip in range(N42):
475        square[ip+3*N4,:] = [-f2,f2-ip*df]
476    square[N-1,:] = [-f2,-f2]
477
478    return square
479
480def p_circle(radii, N=50):
481    """ Function to get a polygon of a circle
482      radii: length of the radii of the circle
483      N: number of points of the polygon
484    """
485    fname = 'p_circle'
486
487    circle = np.zeros((N,2), dtype=np.float)
488
489    dangle = 2.*np.pi/(N-1)
490
491    for ia in range(N):
492        circle[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
493
494    circle[N-1,:] = [0., radii]
495
496    return circle
497
498def p_triangle(p1, p2, p3, N=4):
499    """ Function to provide the polygon of a triangle from its 3 vertices
500      p1: vertex 1 [y,x]
501      p2: vertex 2 [y,x]
502      p3: vertex 3 [y,x]
503      N: number of vertices of the triangle
504    """
505    fname = 'p_triangle'
506
507    triangle = np.zeros((N,2), dtype=np.float)
508
509    N3 = N / 3
510    # 1-2
511    dx = (p2[1]-p1[1])/N3
512    dy = (p2[0]-p1[0])/N3
513    for ip in range(N3):
514        triangle[ip,:] = [p1[0]+ip*dy,p1[1]+ip*dx]
515    # 2-3
516    dx = (p3[1]-p2[1])/N3
517    dy = (p3[0]-p2[0])/N3
518    for ip in range(N3):
519        triangle[ip+N3,:] = [p2[0]+ip*dy,p2[1]+ip*dx]
520    # 3-1
521    N32 = N - 2*N/3
522    dx = (p1[1]-p3[1])/N32
523    dy = (p1[0]-p3[0])/N32
524    for ip in range(N32):
525        triangle[ip+2*N3,:] = [p3[0]+ip*dy,p3[1]+ip*dx]
526
527    triangle[N-1,:] = p1
528
529    return triangle
530
531def p_spiral(loops, eradii, N=1000):
532    """ Function to provide a polygon of an Archimedean spiral
533        FROM: https://en.wikipedia.org/wiki/Spiral
534      loops: number of loops of the spiral
535      eradii: length of the radii of the final spiral
536      N: number of points of the polygon
537    """
538    fname = 'p_spiral'
539
540    spiral = np.zeros((N,2), dtype=np.float)
541
542    dangle = 2.*np.pi*loops/(N-1)
543    dr = eradii*1./(N-1)
544
545    for ia in range(N):
546        radii = dr*ia
547        spiral[ia,:] = [radii*np.sin(ia*dangle), radii*np.cos(ia*dangle)]
548
549    return spiral
550
551def p_reg_polygon(Nv, lf, N=50):
552    """ Function to provide a regular polygon of Nv vertices
553      Nv: number of vertices
554      lf: length of the face
555      N: number of points
556    """
557    fname = 'p_reg_polygon'
558
559    reg_polygon = np.zeros((N,2), dtype=np.float)
560
561    # Number of points per vertex
562    Np = N/Nv
563    # Angle incremental between vertices
564    da = 2.*np.pi/Nv
565    # Radii of the circle according to lf
566    radii = lf*Nv/(2*np.pi)
567
568    iip = 0
569    for iv in range(Nv-1):
570        # Characteristics between vertices iv and iv+1
571        av1 = da*iv
572        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
573        av2 = da*(iv+1)
574        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
575        dx = (v2[1]-v1[1])/Np
576        dy = (v2[0]-v1[0])/Np
577        for ip in range(Np):
578            reg_polygon[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
579
580    # Characteristics between vertices Nv and 1
581
582    # Number of points per vertex
583    Np2 = N - Np*(Nv-1)
584
585    av1 = da*Nv
586    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
587    av2 = 0.
588    v2 = [radii*np.sin(av2), radii*np.cos(av2)]
589    dx = (v2[1]-v1[1])/Np2
590    dy = (v2[0]-v1[0])/Np2
591    for ip in range(Np2):
592        reg_polygon[ip+(Nv-1)*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
593
594    return reg_polygon
595
596def p_reg_star(Nv, lf, freq, vs=0, N=50):
597    """ Function to provide a regular star of Nv vertices
598      Nv: number of vertices
599      lf: length of the face of the regular polygon
600      freq: frequency of union of vertices ('0', for just centered to zero arms)
601      vs: vertex from which start (0 being first [0,lf])
602      N: number of points
603    """
604    fname = 'p_reg_star'
605
606    reg_star = np.zeros((N,2), dtype=np.float)
607
608    # Number of arms of the star
609    if freq != 0 and np.mod(Nv,freq) == 0: 
610        Na = Nv/freq + 1
611    else:
612        Na = Nv
613
614    # Number of points per arm
615    Np = N/Na
616    # Angle incremental between vertices
617    da = 2.*np.pi/Nv
618    # Radii of the circle according to lf
619    radii = lf*Nv/(2*np.pi)
620
621    iip = 0
622    av1 = vs*da
623    for iv in range(Na-1):
624        # Characteristics between vertices iv and iv+1
625        v1 = [radii*np.sin(av1), radii*np.cos(av1)]
626        if freq != 0:
627            av2 = av1 + da*freq
628            v2 = [radii*np.sin(av2), radii*np.cos(av2)]
629        else:
630            v2 = [0., 0.]
631            av2 = av1 + da
632        dx = (v2[1]-v1[1])/(Np-1)
633        dy = (v2[0]-v1[0])/(Np-1)
634        for ip in range(Np):
635            reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
636        if av2 > 2.*np.pi: av1 = av2 - 2.*np.pi
637        else: av1 = av2 + 0.
638
639    iv = Na-1
640    # Characteristics between vertices Na and 1
641    Np2 = N-Np*iv
642    v1 = [radii*np.sin(av1), radii*np.cos(av1)]
643    if freq != 0:
644        av2 = vs*da
645        v2 = [radii*np.sin(av2), radii*np.cos(av2)]
646    else:
647        v2 = [0., 0.]
648    dx = (v2[1]-v1[1])/(Np2-1)
649    dy = (v2[0]-v1[0])/(Np2-1)
650    for ip in range(Np2):
651        reg_star[ip+iv*Np,:] = [v1[0]+dy*ip,v1[1]+dx*ip]
652
653    return reg_star
654
655def p_sinusiode(length=10., amp=5., lamb=3., ival=0., func='sin', N=100):
656    """ Function to get coordinates of a sinusoidal curve
657      length: length of the line (default 10.)
658      amp: amplitude of the peaks (default 5.)
659      lamb: wave longitude (defalult 3.)
660      ival: initial angle (default 0. in degree)
661      func: function to use: (default sinus)
662        'sin': sinus
663        'cos': cosinus
664      N: number of points (default 100)
665    """
666    fname = 'p_sinusiode'
667    availfunc = ['sin', 'cos']
668
669    dx = length/(N-1)
670    ia = ival*np.pi/180.
671    da = 2*np.pi*dx/lamb
672
673    sinusoide = np.zeros((N,2), dtype=np.float)
674    if func == 'sin':
675        for ix in range(N):
676            sinusoide[ix,:] = [amp*np.sin(ia+da*ix),dx*ix]
677    elif func == 'cos':
678        for ix in range(N):
679            sinusoide[ix,:] = [amp*np.cos(ia+da*ix),dx*ix]
680    else:
681        print errormsg
682        print '  ' + fname + ": function '" + func + "' not ready !!"
683        print '    available ones:', availfunc
684        quit(-1)
685
686    sinusoidesecs = ['sinusoide']
687    sinusoidedic = {'sinusoide': [sinusoide, '-', '#000000', 1.]}
688
689    return sinusoide, sinusoidesecs, sinusoidedic
690
691def p_doubleArrow(length=5., angle=45., width=1., alength=0.10, N=50):
692    """ Function to provide an arrow with double lines
693      length: length of the arrow (5. default)
694      angle: angle of the head of the arrow (45., default)
695      width: separation between the two lines (2., default)
696      alength: length of the head (as percentage in excess of width, 0.1 default)
697      N: number of points (50, default)
698    """
699    import numpy.ma as ma
700    function = 'p_doubleArrow'
701
702    doubleArrow = np.zeros((50,2), dtype=np.float)
703    N4 = int((N-3)/4)
704
705    doublearrowdic = {}
706
707    # Arms
708    dx = length/(N4-1)
709    for ix in range(N4-1):
710        doubleArrow[ix,:] = [dx*ix,-width/2.]
711    doublearrowdic['leftarm'] = [doubleArrow[0:N4-1,:], '-', '#000000', 2.]
712    doubleArrow[N4-1,:] = [gen.fillValueF,gen.fillValueF]
713    for ix in range(N4-1):
714        doubleArrow[N4+ix,:] = [dx*ix,width/2.]
715    doublearrowdic['rightarm'] = [doubleArrow[N4:2*N4-1,:], '-', '#000000', 2.]
716    doubleArrow[2*N4-1,:] = [gen.fillValueF,gen.fillValueF]
717
718    # Head
719    N42 = int((N-2 - 2*N4)/2)
720    dx = width*(1.+alength)/(N42-1)
721    for ix in range(N42):
722        doubleArrow[2*N4+ix,:] = [length-dx*ix,-dx*ix]
723    doublearrowdic['lefthead'] = [doubleArrow[2*N4:2*N4+N42,:], '-', '#000000', 2.]
724    doubleArrow[2*N4+N42,:] = [gen.fillValueF,gen.fillValueF]
725
726    N43 = N-2 - 2*N4 - N42 + 1
727    dx = width*(1.+alength)/(N43-1)
728    for ix in range(N43):
729        doubleArrow[2*N4+N42+1+ix,:] = [length-dx*ix,dx*ix]
730    doublearrowdic['rightthead'] = [doubleArrow[2*N4+N42:51,:], '-', '#000000', 2.]
731
732    doubleArrow = ma.masked_equal(doubleArrow, gen.fillValueF)
733    doublearrowsecs = ['leftarm', 'rightarm', 'lefthead', 'righthead']
734
735    return doubleArrow, doublearrowsecs, doublearrowdic
736
737# Combined objects
738##
739
740# FROM: http://www.photographers1.com/Sailing/NauticalTerms&Nomenclature.html
741def zboat(length=10., beam=1., lbeam=0.4, sternbp=0.5):
742    """ Function to define an schematic boat from the z-plane
743      length: length of the boat (without stern, default 10)
744      beam: beam of the boat (default 1)
745      lbeam: length at beam (as percentage of length, default 0.4)
746      sternbp: beam at stern (as percentage of beam, default 0.5)
747    """
748    fname = 'zboat'
749
750    bow = np.array([length, 0.])
751    maxportside = np.array([length*lbeam, -beam])
752    maxstarboardside = np.array([length*lbeam, beam])
753    portside = np.array([0., -beam*sternbp])
754    starboardside = np.array([0., beam*sternbp])
755
756    # forward section
757    fportside = circ_sec(bow,maxportside, length*2)
758    fstarboardside = circ_sec(maxstarboardside, bow, length*2)
759    # aft section
760    aportside = circ_sec(maxportside, portside, length*2)
761    astarboardside = circ_sec(starboardside, maxstarboardside, length*2)
762    # stern
763    stern = circ_sec(portside, starboardside, length*2)
764
765    dpts = stern.shape[0]
766    boat = np.zeros((dpts*5,2), dtype=np.float)
767
768    boat[0:dpts,:] = fportside
769    boat[dpts:2*dpts,:] = aportside
770    boat[2*dpts:3*dpts,:] = stern
771    boat[3*dpts:4*dpts,:] = astarboardside
772    boat[4*dpts:5*dpts,:] = fstarboardside
773
774    fname = 'boat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) + '_lb' +  \
775      str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) + '.dat'
776    if not os.path.isfile(fname):
777        print infmsg
778        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
779        of = open(fname, 'w')
780        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
781          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
782          ' %\n')
783        for ip in range(dpts*5):
784            of.write(str(boat[ip,0]) + ' ' + str(boat[ip,1]) + '\n')
785       
786        of.close()
787        print fname + ": Successfull written '" + fname + "' !!"
788 
789    return boat
790
791def zsailing_boat(length=10., beam=1., lbeam=0.4, sternbp=0.5, lmast=0.6, wmast=0.1, \
792  hsd=5., msd=5., lheads=0.38, lmains=0.55):
793    """ Function to define an schematic sailing boat from the z-plane with sails
794      length: length of the boat (without stern, default 10)
795      beam: beam of the boat (default 1)
796      lbeam: length at beam (as percentage of length, default 0.4)
797      sternbp: beam at stern (as percentage of beam, default 0.5)
798      lmast: position of the mast (as percentage of length, default 0.6)
799      wmast: width of the mast (default 0.1)
800      hsd: head sail direction respect to center line (default 5., -999.99 for upwind)
801      msd: main sail direction respect to center line (default 5., -999.99 for upwind)
802      lheads: length of head sail (as percentage of legnth, defaul 0.38)
803      lmains: length of main sail (as percentage of legnth, defaul 0.55)
804    """
805    import numpy.ma as ma
806    fname = 'zsailing_boat'
807
808    bow = np.array([length, 0.])
809    maxportside = np.array([length*lbeam, -beam])
810    maxstarboardside = np.array([length*lbeam, beam])
811    portside = np.array([0., -beam*sternbp])
812    starboardside = np.array([0., beam*sternbp])
813
814    # forward section
815    fportside = circ_sec(bow,maxportside, length*2)
816    fstarboardside = circ_sec(maxstarboardside, bow, length*2)
817    dpts = fportside.shape[0]
818
819    # aft section
820    aportside = circ_sec(maxportside, portside, length*2)
821    astarboardside = circ_sec(starboardside, maxstarboardside, length*2)
822    # stern
823    stern = circ_sec(portside, starboardside, length*2)
824    # mast
825    mast = p_circle(wmast,N=dpts)
826    mast = mast + [length*lmast, 0.]
827    # head sails
828    lsail = lheads*length
829    if hsd != -999.99:
830        sailsa = np.pi/2. - np.pi*hsd/180.
831        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
832        endsail[0] = length - endsail[0]
833        if bow[1] < endsail[1]:
834            hsail = circ_sec(endsail, bow, lsail*2.15)
835        else:
836            hsail = circ_sec(bow, endsail, lsail*2.15)
837    else:
838        hsail0 = p_sinusiode(length=lsail, amp=0.2, lamb=0.75, N=dpts)
839        hsail = np.zeros((dpts,2), dtype=np.float)
840        hsail[:,0] = hsail0[:,1]
841        hsail[:,1] = hsail0[:,0]
842        hsail = bow - hsail
843
844    # main sails
845    lsail = lmains*length
846    if msd != -999.99:
847        sailsa = np.pi/2. - np.pi*msd/180.
848        begsail = np.array([length*lmast, 0.])
849        endsail = np.array([lsail*np.sin(sailsa), lsail*np.cos(sailsa)])
850        endsail[0] = length*lmast - endsail[0]
851        if endsail[1] < begsail[1]:
852            msail = circ_sec(begsail, endsail, lsail*2.15)
853        else:
854            msail = circ_sec(endsail, begsail, lsail*2.15)
855    else:
856        msail0 = p_sinusiode(length=lsail, amp=0.25, lamb=1., N=dpts)
857        msail = np.zeros((dpts,2), dtype=np.float)
858        msail[:,0] = msail0[:,1]
859        msail[:,1] = msail0[:,0]
860        msail = [length*lmast,0] - msail
861
862    sailingboat = np.zeros((dpts*8+4,2), dtype=np.float)
863
864    sailingboat[0:dpts,:] = fportside
865    sailingboat[dpts:2*dpts,:] = aportside
866    sailingboat[2*dpts:3*dpts,:] = stern
867    sailingboat[3*dpts:4*dpts,:] = astarboardside
868    sailingboat[4*dpts:5*dpts,:] = fstarboardside
869    sailingboat[5*dpts,:] = [gen.fillValueF, gen.fillValueF]
870    sailingboat[5*dpts+1:6*dpts+1,:] = mast
871    sailingboat[6*dpts+1,:] = [gen.fillValueF, gen.fillValueF]
872    sailingboat[6*dpts+2:7*dpts+2,:] = hsail
873    sailingboat[7*dpts+2,:] = [gen.fillValueF, gen.fillValueF]
874    sailingboat[7*dpts+3:8*dpts+3,:] = msail
875    sailingboat[8*dpts+3,:] = [gen.fillValueF, gen.fillValueF]
876
877    sailingboat = ma.masked_equal(sailingboat, gen.fillValueF)
878
879    # Center line extending [fcl] percentage from length on aft and stern
880    fcl = 0.15
881    centerline = np.zeros((dpts,2), dtype=np.float)
882    dl = length*(1.+fcl*2.)/(dpts-1)
883    centerline[:,0] = np.arange(-length*fcl, length*(1. + fcl)+dl, dl)
884
885    # correct order of sections
886    sailingboatsecs = ['fportside', 'aportside', 'stern', 'astarboardside',          \
887      'fstarboardside', 'mast', 'hsail', 'msail', 'centerline']
888    # dictionary with sections [polygon_vertices, line_type, line_color, line_width]
889    dicsailingboat = {'fportside': [fportside, '-', '#8A5900', 2.],                  \
890      'aportside': [aportside, '-', '#8A5900', 2.],                                  \
891      'stern': [stern, '-', '#8A5900', 2.],                                          \
892      'astarboardside': [astarboardside, '-', '#8A5900', 2.],                        \
893      'fstarboardside': [fstarboardside, '-', '#8A5900', 2.],                        \
894      'mast': [mast, '-', '#8A5900', 2.], 'hsail': [hsail, '-', '#AAAAAA', 1.],      \
895      'msail': [msail, '-', '#AAAAAA', 1.],                                          \
896      'centerline': [centerline, '-.', '#AA6464', 1.5]}
897   
898    fname = 'sailboat_L' + str(int(length*100.)) + '_B' + str(int(beam*100.)) +      \
899      '_lb' + str(int(lbeam*100.)) + '_sb' + str(int(sternbp*100.)) +                \
900      '_lm' + str(int(lmast*100.)) + '_wm' + str(int(wmast)) +                       \
901      '_hsd' + str(int(hsd)) + '_hs' + str(int(lheads*100.)) +                       \
902      '_ms' + str(int(lheads*100.)) + '_msd' + str(int(msd)) +'.dat'
903    if not os.path.isfile(fname):
904        print infmsg
905        print '  ' + fname + ": writting boat coordinates file '" + fname + "' !!"
906        of = open(fname, 'w')
907        of.write('# boat file with Length: ' + str(length) +' max_beam: '+str(beam)+ \
908          'length_at_max_beam:' + str(lbeam) + '% beam at stern: ' + str(sternbp)+   \
909          ' % mast position: '+ str(lmast) + ' % mast width: ' + str(wmast) + ' ' +  \
910          ' head sail direction:' + str(hsd) + ' head sail length: ' + str(lheads) + \
911          ' %' + ' main sail length' + str(lmains) + ' main sail direction:' +       \
912          str(msd) +'\n')
913        for ip in range(dpts*5):
914            of.write(str(sailingboat[ip,0]) + ' ' + str(sailingboat[ip,1]) + '\n')
915       
916        of.close()
917        print fname + ": Successfull written '" + fname + "' !!"
918 
919    return sailingboat, sailingboatsecs, dicsailingboat
920
921def zisland1(mainpts= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],\
922  [2.8, -0.1], [0.1,-0.6]], dtype=np.float), radfrac=3., N=200):
923    """ Function to draw an island from z-axis as the union of a series of points by
924        circular segments
925      mainpts: main points of the island (clockwise ordered, to be joined by
926        circular segments of radii as the radfrac factor of the distance between
927        consecutive points)
928          * default= np.array([[-0.1,0.], [-1.,1.], [-0.8,1.2], [0.1,0.6], [1., 0.9],
929            [2.8, -0.1], [0.1,-0.6]], dtype=np.float)
930      radfrac: multiplicative factor of the distance between consecutive points to
931        draw the circular segment (3., default)
932      N: number of points (200, default)
933    """
934    import numpy.ma as ma
935    fname = 'zisland1'
936
937    island1 = np.ones((N,2), dtype=np.float)*gen.fillValueF
938
939    # Coastline
940   
941    island1 = join_circ_sec(mainpts, radfrac=3., N=200)
942
943    islandsecs = ['coastline']
944    islanddic = {'coastline': [island1, '-', '#161616', 2.]}
945
946    island1 = ma.masked_equal(island1, gen.fillValueF)
947
948    return island1, islandsecs, islanddic
949
950def write_join_poly(polys, flname='join_polygons.dat'):
951    """ Function to write an ASCII file with the combination of polygons
952      polys: dictionary with the names of the different polygons
953      flname: name of the ASCII file
954    """
955    fname = 'write_join_poly'
956
957    of = open(flname, 'w')
958
959    for polyn in polys.keys():
960        vertices = polys[polyn]
961        Npts = vertices.shape[0]
962        for ip in range(Npts):
963            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')
964
965    of.close()
966
967    return
968
969def read_join_poly(flname='join_polygons.dat'):
970    """ Function to read an ASCII file with the combination of polygons
971      flname: name of the ASCII file
972    """
973    fname = 'read_join_poly'
974
975    of = open(flname, 'r')
976
977    polys = {}
978    polyn = ''
979    poly = []
980    for line in of:
981        if len(line) > 1: 
982            linevals = line.replace('\n','').split(' ')
983            if polyn != linevals[0]:
984                if len(poly) > 1:
985                    polys[polyn] = np.array(poly)
986                polyn = linevals[0]
987                poly = []
988                poly.append([np.float(linevals[2]), np.float(linevals[1])])
989            else:
990                poly.append([np.float(linevals[2]), np.float(linevals[1])])
991
992    of.close()
993    polys[polyn] = np.array(poly)
994
995    return polys
996
997####### ####### ##### #### ### ## #
998# Plotting
999
1000def plot_sphere(iazm=-60., iele=30., dist=10., Npts=100, radii=10,                   \
1001  drwsfc=[True,True], colsfc=['#AAAAAA','#646464'],                                  \
1002  drwxline = True, linex=[':','b',2.], drwyline = True, liney=[':','r',2.],          \
1003  drwzline = True, linez=['-.','g',2.], drwxcline=[True,True],                       \
1004  linexc=[['-','#646400',1.],['--','#646400',1.]],                                   \
1005  drwequator=[True,True], lineeq=[['-','#AA00AA',1.],['--','#AA00AA',1.]],           \
1006  drwgreeenwhich=[True,True], linegw=[['-','k',1.],['--','k',1.]]):
1007    """ Function to plot an sphere and determine which standard lines will be also
1008        drawn
1009      iazm: azimut of the camera form the sphere
1010      iele: elevation of the camera form the sphere
1011      dist: distance of the camera form the sphere
1012      Npts: Resolution for the sphere
1013      radii: radius of the sphere
1014      drwsfc: whether 'up' and 'down' portions of the sphere should be drawn
1015      colsfc: colors of the surface of the sphere portions ['up', 'down']
1016      drwxline: whether x-axis line should be drawn
1017      linex: properties of the x-axis line ['type', 'color', 'wdith']
1018      drwyline: whether y-axis line should be drawn
1019      liney: properties of the y-axis line ['type', 'color', 'wdith']
1020      drwzline: whether z-axis line should be drawn
1021      linez: properties of the z-axis line ['type', 'color', 'wdith']
1022      drwequator: whether 'front' and 'back' portions of the Equator should be drawn
1023      lineeq: properties of the lines 'front' and 'back' of the Equator
1024      drwgreeenwhich: whether 'front', 'back' portions of Greenqhich should be drawn
1025      linegw: properties of the lines 'front' and 'back' Greenwhich
1026      drwxcline: whether 'front', 'back' 90 line (lon=90., lon=270.) should be drawn
1027      linexc: properties of the lines 'front' and 'back' for the 90 line
1028    """
1029    fname = 'plot_sphere'
1030
1031    iazmrad = iazm*np.pi/180.
1032    ielerad = iele*np.pi/180.
1033
1034    # 3D surface Sphere
1035    sfcsphereu, sfcsphered = surface_sphere(radii,Npts)
1036   
1037    # greenwhich
1038    if iazmrad > np.pi/2. and iazmrad < 3.*np.pi/2.:
1039        ia=np.pi-ielerad
1040    else:
1041        ia=0.-ielerad
1042    ea=ia+np.pi
1043    da = (ea-ia)/(Npts-1)
1044    beta = np.arange(ia,ea+da,da)[0:Npts]
1045    alpha = np.zeros((Npts), dtype=np.float)
1046    greenwhichc = spheric_line(radii,alpha,beta)
1047    ia=ea+0.
1048    ea=ia+np.pi
1049    da = (ea-ia)/(Npts-1)
1050    beta = np.arange(ia,ea+da,da)[0:Npts]
1051    greenwhichd = spheric_line(radii,alpha,beta)
1052
1053    # Equator
1054    ia=np.pi-iazmrad/2.
1055    ea=ia+np.pi
1056    da = (ea-ia)/(Npts-1)
1057    alpha = np.arange(ia,ea+da,da)[0:Npts]
1058    beta = np.zeros((Npts), dtype=np.float)
1059    equatorc = spheric_line(radii,alpha,beta)
1060    ia=ea+0.
1061    ea=ia+np.pi
1062    da = (ea-ia)/(Npts-1)
1063    alpha = np.arange(ia,ea+da,da)[0:Npts]
1064    equatord = spheric_line(radii,alpha,beta)
1065
1066    # 90 line
1067    if iazmrad > np.pi and iazmrad < 2.*np.pi:
1068        ia=3.*np.pi/2. + ielerad
1069    else:
1070        ia=np.pi/2. - ielerad
1071    if ielerad < 0.:
1072        ia = ia + np.pi
1073    ea=ia+np.pi
1074    da = (ea-ia)/(Npts-1)
1075    beta = np.arange(ia,ea+da,da)[0:Npts]
1076    alpha = np.ones((Npts), dtype=np.float)*np.pi/2.
1077    xclinec = spheric_line(radii,alpha,beta)
1078    ia=ea+0.
1079    ea=ia+np.pi
1080    da = (ea-ia)/(Npts-1)
1081    beta = np.arange(ia,ea+da,da)[0:Npts]
1082    xclined = spheric_line(radii,alpha,beta)
1083
1084    # x line
1085    xline = np.zeros((2,3), dtype=np.float)
1086    xline[0,:] = position_sphere(radii, 0., 0.)
1087    xline[1,:] = position_sphere(radii, np.pi, 0.)
1088
1089    # y line
1090    yline = np.zeros((2,3), dtype=np.float)
1091    yline[0,:] = position_sphere(radii, np.pi/2., 0.)
1092    yline[1,:] = position_sphere(radii, 3*np.pi/2., 0.)
1093
1094    # z line
1095    zline = np.zeros((2,3), dtype=np.float)
1096    zline[0,:] = position_sphere(radii, 0., np.pi/2.)
1097    zline[1,:] = position_sphere(radii, 0., -np.pi/2.)
1098
1099    fig = plt.figure()
1100    ax = fig.gca(projection='3d')
1101
1102    # Sphere surface
1103    if drwsfc[0]:
1104        ax.plot_surface(sfcsphereu[0,:,:], sfcsphereu[1,:,:], sfcsphereu[2,:,:],     \
1105          color=colsfc[0])
1106    if drwsfc[1]:
1107        ax.plot_surface(sfcsphered[0,:,:], sfcsphered[1,:,:], sfcsphered[2,:,:],     \
1108          color=colsfc[1])
1109
1110    # greenwhich
1111    linev = linegw[0]
1112    if drwgreeenwhich[0]:
1113        ax.plot(greenwhichc[:,0], greenwhichc[:,1], greenwhichc[:,2], linev[0],      \
1114          color=linev[1], linewidth=linev[2],  label='Greenwich')
1115    linev = linegw[1]
1116    if drwgreeenwhich[1]:
1117        ax.plot(greenwhichd[:,0], greenwhichd[:,1], greenwhichd[:,2], linev[0],      \
1118          color=linev[1], linewidth=linev[2])
1119
1120    # Equator
1121    linev = lineeq[0]
1122    if drwequator[0]:
1123        ax.plot(equatorc[:,0], equatorc[:,1], equatorc[:,2], linev[0],               \
1124          color=linev[1], linewidth=linev[2], label='Equator')
1125    linev = lineeq[1]
1126    if drwequator[1]:
1127        ax.plot(equatord[:,0], equatord[:,1], equatord[:,2], linev[0],               \
1128          color=linev[1], linewidth=linev[2])
1129
1130    # 90line
1131    linev = linexc[0]
1132    if drwxcline[0]:
1133        ax.plot(xclinec[:,0], xclinec[:,1], xclinec[:,2], linev[0], color=linev[1],  \
1134          linewidth=linev[2], label='90-line')
1135    linev = linexc[1]
1136    if drwxcline[1]:
1137        ax.plot(xclined[:,0], xclined[:,1], xclined[:,2], linev[0], color=linev[1],  \
1138          linewidth=linev[2])
1139
1140    # x line
1141    linev = linex
1142    if drwxline:
1143        ax.plot([xline[0,0],xline[1,0]], [xline[0,1],xline[1,1]],                    \
1144          [xline[0,2],xline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='xline')
1145
1146    # y line
1147    linev = liney
1148    if drwyline:
1149        ax.plot([yline[0,0],yline[1,0]], [yline[0,1],yline[1,1]],                    \
1150          [yline[0,2],yline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='yline')
1151
1152    # z line
1153    linev = linez
1154    if drwzline:
1155        ax.plot([zline[0,0],zline[1,0]], [zline[0,1],zline[1,1]],                    \
1156          [zline[0,2],zline[1,2]], linev[0], color=linev[1], linewidth=linev[2],  label='zline')
1157
1158    plt.legend()
1159
1160    return fig, ax
Note: See TracBrowser for help on using the repository browser.