Changeset 2434 in lmdz_wrf for trunk/tools

Timestamp:

Apr 14, 2019, 10:08:57 PM (6 years ago)

Author:

lfita

Message:

Working version of `circ_sec'

File:

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

trunk/tools/geometry_tools.py

@@ -20 +20 @@
 # deg_deci: Function to pass from degrees [deg, minute, sec] to decimal angles [rad]
 # dist_points: Function to provide the distance between two points
-# multi_rotate_2D: Function to rotate multiple vectors by a certain angle in the plain
+# rotate_polygon_2D: Function to rotate 2D plain the vertices of a polygon
 # position_sphere: Function to tranform fom a point in lon, lat deg coordinates to 
 #   cartesian coordinates over an sphere
@@ -122 +122 @@
     return rotvector
 
-def multi_rotate_2D(vectors, angle):
-    """ Function to rotate multiple vectors by a certain angle in the plain
+def rotate_polygon_2D(vectors, angle):
+    """ Function to rotate 2D plain the vertices of a polygon
       line= matrix of vectors to rotate
       angle= angle to rotate [rad]
@@ -131 +131 @@
     >>> square[2,:] = [0.5,0.5]
     >>> square[3,:] = [-0.5,0.5]
-    >>> multi_rotate_2D(square, np.pi/4.)
+    >>> rotate_polygon_2D(square, np.pi/4.)
     [[-0.70710678  0.        ]
      [ 0.         -0.70710678]
@@ -137 +137 @@
      [ 0.          0.70710678]]
     """
-    fname = 'multi_rotate_2D'
+    fname = 'rotate_polygon_2D'
 
     rotvecs = np.zeros(vectors.shape, dtype=np.float)
@@ -318 +318 @@
 
     # angle of the circular section joining points
-    tottheta = 2.*np.arcsin(distAB/2./radii)
+    alpha = 2.*np.arcsin((distAB/2.)/radii)
 
     # center along coincident bisection of the union
@@ -324 +324 @@
     ycc = 0.
 
-    # Angle of the points
-    pttheta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
-
-    # rotating the position of the center
-    rotang = np.pi/2.-pttheta
-
-    newcc = rotate_line2D(np.array([[ycc,xcc], [0.,0.]]), rotang)
-    print newcc
-
-    quit()
-
-    dtheta = np.abs(tottheta)/(Nang-1)
-    if sign == 1:
-        dtheta = dtheta*(-1.)
-
-    print 'Lluis tottheta:', tottheta*180./np.pi, 'dtheta:', dtheta*180./np.pi, 'c:', xc,yc
-
+    # Getting the arc of the circle at the x-axis
+    dalpha = alpha/(Nang-1)
     circ_sec = np.zeros((Nang,2), dtype=np.float)
     for ia in range(Nang):
-        theta = dtheta*ia
-        x = radii*np.cos(theta)
-        y = radii*np.sin(theta)
-
-        circ_sec[ia,:] = [y+yc,x+xc]
-        print ia, 'Lluis xy:', x,y, 'circ', circ_sec[ia,:]
+        alpha = dalpha*ia
+        x = radii*np.cos(alpha)
+        y = radii*np.sin(alpha)
+
+        circ_sec[ia,:] = [y+ycc,x+xcc]
     
+    # Angle of the points
+    theta = np.arctan2(ptB[0]-ptA[0],ptB[1]-ptA[1])
+
+    # rotating angle of the circ
+    rotangle = theta + 3.*np.pi/2. - alpha/2.
+
+    #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)
+
+    # Moving arc to the ptA
+    circ_sec = rotcirc_sec + ptA
+
     return circ_sec