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

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

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