Changeset 2508 in lmdz_wrf for trunk

Timestamp:

May 4, 2019, 10:14:47 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `buoy1': Function to draw a buoy as superposition of prism and section of ball
• 'arc' into `circ_sec' in order to determine which arc should be taken, the 'short' or the 'long' ('short', default)

File:

• trunk/tools/geometry_tools.py (modified) (10 diffs)

trunk/tools/geometry_tools.py

@@ -42 +42 @@
 
 ## Shapes/objects
+# buoy1: Function to draw a buoy as superposition of prism and section of ball
 # circ_sec: Function union of point A and B by a section of a circle
 # ellipse_polar: Function to determine an ellipse from its center and polar coordinates
@@ -310 +311 @@
         p2 = lpoints[ip+1]
         dps = dist_points(p1, p2)
-        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, Np)
+        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, 'short', Np)
 
     Np2 = N - (Npts-1)*Np
@@ -316 +317 @@
     p2 = lpoints[0]
     dps = dist_points(p1, p2)
-    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., Np2)
+    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., 'short', Np2)
 
     return jcirc_sec
@@ -346 +347 @@
         dps = dist_points(p1, p2)
         angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
-        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, Np)
+        jcirc_sec[Np*ip:Np*(ip+1),:] = circ_sec(p1, p2, dps*radfrac, 'short', Np)
         drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
         for iip in range(Np*ip,Np*(ip+1)):
@@ -356 +357 @@
     dps = dist_points(p1, p2)
     angle = np.arctan2(p2[0]-p1[0], p2[1]-p1[1]) + np.pi/2.
-    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., Np2)
+    jcirc_sec[(Npts-1)*Np:N,:] = circ_sec(p1, p2, dps*3., 'short', Np2)
     drand = Lrand*np.array([np.sin(angle), np.cos(angle)])
     for iip in range(Np*(Npts-1),N):
@@ -363 +364 @@
     return jcirc_sec
 
+def write_join_poly(polys, flname='join_polygons.dat'):
+    """ Function to write an ASCII file with the combination of polygons
+      polys: dictionary with the names of the different polygons
+      flname: name of the ASCII file
+    """
+    fname = 'write_join_poly'
+
+    of = open(flname, 'w')
+
+    for polyn in polys.keys():
+        vertices = polys[polyn]
+        Npts = vertices.shape[0]
+        for ip in range(Npts):
+            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')
+
+    of.close()
+
+    return
+
+def read_join_poly(flname='join_polygons.dat'):
+    """ Function to read an ASCII file with the combination of polygons
+      flname: name of the ASCII file
+    """
+    fname = 'read_join_poly'
+
+    of = open(flname, 'r')
+
+    polys = {}
+    polyn = ''
+    poly = []
+    for line in of:
+        if len(line) > 1: 
+            linevals = line.replace('\n','').split(' ')
+            if polyn != linevals[0]:
+                if len(poly) > 1:
+                    polys[polyn] = np.array(poly)
+                polyn = linevals[0]
+                poly = []
+                poly.append([np.float(linevals[2]), np.float(linevals[1])])
+            else:
+                poly.append([np.float(linevals[2]), np.float(linevals[1])])
+
+    of.close()
+    polys[polyn] = np.array(poly)
+
+    return polys
 
 ####### ###### ##### #### ### ## #
@@ -442 +489 @@
     return hyperbola_p, hyperbola_n
 
-def circ_sec(ptA, ptB, radii, Nang=100):
+def circ_sec(ptA, ptB, radii, arc='short', Nang=100):
     """ Function union of point A and B by a section of a circle
       ptA= coordinates od the point A [yA, xA]
       ptB= coordinates od the point B [yB, xB]
       radii= radi of the circle to use to unite the points
+      arc= which arc to be used ('short', default)
+        'short': shortest angle between points
+        'long': largest angle between points
       Nang= amount of angles to use
     """
@@ -471 +521 @@
 
     # Getting the arc of the circle at the x-axis
-    dalpha = alpha/(Nang-1)
+    if arc == 'short':
+        dalpha = alpha/(Nang-1)
+    else:
+        dalpha = (2.*np.pi - alpha)/(Nang-1)
     circ_sec = np.zeros((Nang,2), dtype=np.float)
     for ia in range(Nang):
@@ -488 +541 @@
     #print 'alpha:', alpha*180./np.pi, 'theta:', theta*180./np.pi, 'rotangle:', rotangle*180./np.pi
   
-
     # rotating the arc along the x-axis
     rotcirc_sec = rotate_polygon_2D(circ_sec, rotangle)
@@ -996 +1048 @@
     return island1, islandsecs, islanddic
 
-def write_join_poly(polys, flname='join_polygons.dat'):
-    """ Function to write an ASCII file with the combination of polygons
-      polys: dictionary with the names of the different polygons
-      flname: name of the ASCII file
-    """
-    fname = 'write_join_poly'
-
-    of = open(flname, 'w')
-
-    for polyn in polys.keys():
-        vertices = polys[polyn]
-        Npts = vertices.shape[0]
-        for ip in range(Npts):
-            of.write(polyn+' '+str(vertices[ip,1]) + ' ' + str(vertices[ip,0]) + '\n')
-
-    of.close()
-
-    return
-
-def read_join_poly(flname='join_polygons.dat'):
-    """ Function to read an ASCII file with the combination of polygons
-      flname: name of the ASCII file
-    """
-    fname = 'read_join_poly'
-
-    of = open(flname, 'r')
-
-    polys = {}
-    polyn = ''
-    poly = []
-    for line in of:
-        if len(line) > 1: 
-            linevals = line.replace('\n','').split(' ')
-            if polyn != linevals[0]:
-                if len(poly) > 1:
-                    polys[polyn] = np.array(poly)
-                polyn = linevals[0]
-                poly = []
-                poly.append([np.float(linevals[2]), np.float(linevals[1])])
-            else:
-                poly.append([np.float(linevals[2]), np.float(linevals[1])])
-
-    of.close()
-    polys[polyn] = np.array(poly)
-
-    return polys
+def buoy1(height=5., width=10., bradii=1.75, bfrac=0.8, N=200):
+    """ Function to draw a buoy as superposition of prism and section of ball
+      height: height of the prism (5., default) 
+      width: width of the prism (10., default)
+      bradii: radii of the ball (1.75, default)
+      bfrac: fraction of the ball above the prism (0.8, default)
+      N: total number of points of the buoy (200, default)
+
+    """
+    fname = 'buoy1'
+
+    buoy = np.zeros((N,2), dtype=np.float)
+
+    NNp = 0
+    iip = 0
+    # Base
+    ix = -width/2.
+    iy = 0.
+    Np = 4
+    dx = width/(Np-1)
+    dy = 0.
+    for ip in range(Np):
+        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
+    NNp = NNp + Np
+    iip = NNp
+
+    # right lateral
+    ix = width/2.
+    iy = 0.
+    Np = 2
+    dx = 0.
+    dy = height/(Np-1)
+    for ip in range(Np):
+        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
+    NNp = NNp + Np
+    iip = NNp
+
+    # right upper
+    ix = width/2.
+    iy = height
+    Np = 4
+    dx = -(width/2.-bradii*bfrac)/(Np-1)
+    dy = 0.
+    for ip in range(Np):
+        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
+    NNp = NNp + Np
+    iip = NNp
+
+    # ball
+    p1 = np.array([height, -bradii*bfrac])
+    p2 = np.array([height, bradii*bfrac])
+    Np = N - 2*(NNp)
+    buoy[iip:iip+Np,:] = circ_sec(p2, p1, 2.*bradii, 'long', Np)
+    NNp = NNp + Np
+    iip = NNp
+
+    # left upper
+    ix = -bradii*bfrac
+    iy = height
+    Np = 4
+    dx = -(width/2.-bradii*bfrac)/(Np-1)
+    dy = 0.
+    for ip in range(Np):
+        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
+    NNp = NNp + Np
+    iip = NNp
+
+    # left lateral
+    ix = -width/2.
+    iy = height
+    Np = 2
+    dx = 0.
+    dy = -height/(Np-1)
+    for ip in range(Np):
+        buoy[iip+ip,:] = [iy+dy*ip,ix+dx*ip]
+    NNp = NNp + Np
+    iip = NNp
+
+    buoysecs = ['base']
+    buoydic = {'base': [buoy, '-', 'k', 1.5]}
+
+    return buoy, buoysecs, buoydic
 
 ####### ####### ##### #### ### ## #