#! /usr/bin/env python3 # def equilateral_triangle ( ): #*****************************************************************************80 # ## equilateral_triangle() returns the vertices of an equilateral triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Output: # # real v[3,2]: the vertices. # import numpy as np v = np.array ( [ \ [ 0.5, 0.0 ], \ [ 0.0, 0.5 * np.sqrt ( 3.0 ) ], \ [ -0.5, 0.0 ] ] ) return v def left_right ( p, q, r ): #*****************************************************************************80 # ## left_right() reports whether point r is left or right of the line p->q. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real p[2], q[2]: points defining a line. # # real r[2]: a point whose orientation is desired. # # Output: # # real value: # * negative, to the left. # * zero, on the line (very rare) # * positive, to the right # import numpy as np # # Compute the vector cross product (r-p)x(q-p) # value = \ ( r[0] - p[0] ) * ( q[1] - p[1] ) \ - ( r[1] - p[1] ) * ( q[0] - p[0] ) return value def left_right_test ( ): #*****************************************************************************80 # ## left_right_test() tests left_right(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # import matplotlib.pyplot as plt import numpy as np from numpy.random import default_rng rng = default_rng ( ) p = rng.standard_normal ( 2 ) q = rng.standard_normal ( 2 ) n = 20 r = rng.standard_normal ( [ n, 2 ] ) plt.clf ( ) plt.plot ( [ p[0], q[0] ], [ p[1], q[1] ], 'k-', linewidth = 3 ) plt.plot ( p[0], p[1], 'g.', markersize = 30 ) plt.plot ( q[0], q[1], 'r.', markersize = 30 ) # # Use . and o, so the plot is still legible in black and white. # for i in range ( 0, n ): value = left_right ( p, q, r[i] ) if ( value < 0.0 ): plt.plot ( r[i,0], r[i,1], 'm*', markersize = 20 ) elif ( 0.0 < value ): plt.plot ( r[i,0], r[i,1], 'co', markersize = 10 ) else: plt.plot ( r[i,0], r[i,1], 'b.', markersize = 15 ) plt.grid ( True ) # # Stupid "axis('equal') doesn't... # # plt.axis ( 'equal' ) plt.gca().set_aspect ( 'equal', adjustable = 'box' ) filename = 'left_right.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + "'" ) plt.show ( block = False ) plt.close ( ) return def obtuse_triangle ( ): #*****************************************************************************80 # ## obtuse_triangle() returns the vertices of an obtuse triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Output: # # real v[3,2]: the vertices. # import numpy as np v = np.array ( [ \ [ 0.0, 0.0 ], \ [ 1.0, 0.0 ], \ [ -1.0, 1.0 ] ] ) return v def reference_triangle ( ): #*****************************************************************************80 # ## reference_triangle() returns the vertices of the reference triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Output: # # real v[3,2]: the vertices. # import numpy as np v = np.array ( [ \ [ 0.0, 0.0 ], \ [ 1.0, 0.0 ], \ [ 0.0, 1.0 ] ] ) return v def scalene_triangle ( ): #*****************************************************************************80 # ## scalene_triangle() returns the vertices of a scalene triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Output: # # real v[3,2]: the vertices. # import numpy as np v = np.array ( [ \ [ 6.0, 0.0 ], \ [ 1.0, 1.0 ], \ [ 5.0, 3.0 ] ] ) return v def triangle_angles ( v ): #*****************************************************************************80 # ## triangle_angles() returns angles of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real angle[3]: the angles. # import numpy as np s = triangle_sides ( v ) angle = np.zeros ( 3 ) for i in range ( 0, 3 ): angle[i] = - 0.5 * ( s[(i+1)%3]**2 - s[i]**2 - s[(i+2)%3]**2 ) \ / s[i] / s[(i+2)%3] angle = np.arccos ( angle ) return angle def triangle_area_cross_product ( v ): #*****************************************************************************80 # ## triangle_area_cross_product() returns the area of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real area: the area. # import numpy as np area = 0.5 * ( \ ( v[2,0] - v[1,0] ) * ( v[0,1] - v[1,1] ) \ - ( v[2,1] - v[1,1] ) * ( v[0,0] - v[1,0] ) ) return area def triangle_area_euclid ( v ): #*****************************************************************************80 # ## triangle_area_euclid() returns the area of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real area: the area. # import numpy as np base = np.linalg.norm ( v[1] - v[0] ) base_unit = ( v[1] - v[0] ) / base vector = ( v[2] - v[0] ) dot = np.dot ( vector, base_unit ) vector_parallel = dot * base_unit vector_perp = vector - vector_parallel height = np.linalg.norm ( vector_perp ) area = 0.5 * base * height return area def triangle_area_hero ( v ): #*****************************************************************************80 # ## triangle_area_hero() returns the area of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real area: the area. # import numpy as np s = triangle_sides ( v ) area = 0.25 * np.sqrt ( \ ( s[0] + s[1] + s[2] ) * \ ( - s[0] + s[1] + s[2] ) * \ ( s[0] - s[1] + s[2] ) * \ ( s[0] + s[1] - s[2] ) ) return area def triangle_centroid ( v ): #*****************************************************************************80 # ## triangle_centroid() returns the centroid of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real C: the centroid. # import numpy as np c = np.sum ( v[:], axis = 0 ) / 3.0 return c def triangle_contains ( v, pvec ): #*****************************************************************************80 # ## triangle_contains() determines if a triangle contains a point. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # real pvec[n,2]: the points. # # Output: # # bool contains: # True: the triangle contains the point # False: the triangle does not contain the point. # import numpy as np pvec = np.atleast_2d ( pvec ) n = pvec.shape[0] contains = np.zeros ( n, dtype = bool ) for i in range ( 0, n ): contains[i] = True for j in range ( 0, 3 ): if ( 0 < left_right ( v[j], v[(j+1)%3], pvec[i] ) ): contains[i] = False break return contains def triangle_contains_test ( ): #*****************************************************************************80 # ## triangle_contains_test() tests triangle_contains(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # import matplotlib.pyplot as plt import numpy as np from numpy.random import default_rng v = obtuse_triangle ( ) v = 4.0 * v rng = default_rng ( ) n = 100 r = 2.0 * rng.standard_normal ( [ n, 2 ] ) I = [ 0, 1, 2, 0 ] plt.clf ( ) plt.plot ( v[I,0], v[I,1], 'k-', linewidth = 3 ) # # Use . and o so the plot is still legible if printed in black and white. # for i in range ( 0, n ): value = triangle_contains ( v, r[i] ) if ( value ): plt.plot ( r[i,0], r[i,1], 'm*', markersize = 20 ) else: plt.plot ( r[i,0], r[i,1], 'co', markersize = 10 ) plt.grid ( True ) # # Stupid "axis('equal') doesn't... # # plt.axis ( 'equal' ) plt.gca().set_aspect ( 'equal', adjustable = 'box' ) filename = 'triangle_contains.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + "'" ) plt.show ( block = False ) plt.close ( ) return def triangle_integral_centroid ( v, f ): #*****************************************************************************80 # ## triangle_integral_centroid() estimates an integral over a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # function f(x,y): the function to be integrated. # # Output: # # real value: the estimated integral. # p = triangle_centroid ( v ) fp = f ( p[0], p[1] ) a = triangle_area_cross_product ( v ) value = a * fp return value def triangle_integral_centroid_test ( ): #*****************************************************************************80 # ## triangle_integral_centroid_test() tests triangle_integral_centroid(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'triangle_integral_centroid_test():' ) print ( ' triangle_integral_centroid() uses the centroid rule' ) print ( ' to estimate an integral f(x,y) over a triangle T.' ) exact = - 37.0 / 6.0 v = reference_triangle ( ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_centroid ( v, f ) print ( value, abs ( value - exact ) ) exact = - 935.0 / 3.0 v = np.array ( [ \ [ 0.0, 3.0 ], \ [ 1.0, 1.0 ], \ [ 5.0, 3.0 ] ] ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_centroid ( v, f ) print ( value, abs ( value - exact ) ) return def triangle_integral_estimate ( v, f, n ): #*****************************************************************************80 # ## triangle_integral_estimate() estimates an integral over a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # function f(x,y): the function to be integrated. # # integer n: the number of sample points. # # Output: # # real value: the estimated integral. # p = triangle_point_random ( v, n ) fp = f ( p[:,0], p[:,1] ) a = triangle_area_cross_product ( v ) value = a * sum ( fp ) / n return value def triangle_integral_estimate_test ( ): #*****************************************************************************80 # ## triangle_integral_estimate_test() tests triangle_integral_estimate(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'triangle_integral_estimate_test():' ) print ( ' triangle_integral_estimate() uses a Monte-Carlo approach' ) print ( ' to estimate an integral f(x,y) over a triangle T.' ) exact = - 935.0 / 3.0 v = np.array ( [ \ [ 0.0, 3.0 ], \ [ 1.0, 1.0 ], \ [ 5.0, 3.0 ] ] ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' n estimate error' ) print ( '' ) for n in [ 10, 100, 1000, 10000, 10000 ]: value = triangle_integral_estimate ( v, f, n ) print ( n, value, abs ( value - exact ) ) return def triangle_integral_strang ( v, f ): #*****************************************************************************80 # ## triangle_integral_strang() estimates an integral over a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # function f(x,y): the function to be integrated. # # Output: # # real value: the estimated integral. # import numpy as np p = np.zeros ( [ 3, 2 ] ) p[0] = 4.0 / 6.0 * v[0] + 1.0 / 6.0 * v[1] + 1.0 / 6.0 * v[2] p[1] = 1.0 / 6.0 * v[0] + 4.0 / 6.0 * v[1] + 1.0 / 6.0 * v[2] p[2] = 1.0 / 6.0 * v[0] + 1.0 / 6.0 * v[1] + 4.0 / 6.0 * v[2] fp = f ( p[:,0], p[:,1] ) a = triangle_area_cross_product ( v ) value = a * sum ( fp ) / 3.0 return value def triangle_integral_strang_test ( ): #*****************************************************************************80 # ## triangle_integral_strang_test() tests triangle_integral_strang(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'triangle_integral_strang_test():' ) print ( ' triangle_integral_strang() uses a Strang rule' ) print ( ' to estimate an integral f(x,y) over a triangle T.' ) exact = - 37.0 / 6.0 v = reference_triangle ( ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_strang ( v, f ) print ( value, abs ( value - exact ) ) exact = - 935.0 / 3.0 v = np.array ( [ \ [ 0.0, 3.0 ], \ [ 1.0, 1.0 ], \ [ 5.0, 3.0 ] ] ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_strang ( v, f ) print ( value, abs ( value - exact ) ) return def triangle_integral_vertex ( v, f ): #*****************************************************************************80 # ## triangle_integral_vertex() estimates an integral over a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # function f(x,y): the function to be integrated. # # Output: # # real value: the estimated integral. # fp = f ( v[:,0], v[:,1] ) a = triangle_area_cross_product ( v ) value = a * sum ( fp ) / 3.0 return value def triangle_integral_vertex_test ( ): #*****************************************************************************80 # ## triangle_integral_vertex_test() tests triangle_integral_vertex(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'triangle_integral_vertex_test():' ) print ( ' triangle_integral_vertex() uses the vertex rule' ) print ( ' to estimate an integral f(x,y) over a triangle T.' ) exact = - 37.0 / 6.0 v = reference_triangle ( ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_vertex ( v, f ) print ( value, abs ( value - exact ) ) exact = - 935.0 / 3.0 v = np.array ( [ \ [ 0.0, 3.0 ], \ [ 1.0, 1.0 ], \ [ 5.0, 3.0 ] ] ) f = lambda x, y: 6.0 * x**2 - 40.0 * y print ( '' ) print ( ' estimate error' ) print ( '' ) value = triangle_integral_vertex ( v, f ) print ( value, abs ( value - exact ) ) return def triangle_perimeter ( v ): #*****************************************************************************80 # ## triangle_perimeter() returns perimeter of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real p: the perimeter. # import numpy as np s = triangle_sides ( v ) p = np.sum ( s ) return p def triangle_plot ( v, filename ): #*****************************************************************************80 # ## triangle_plot ( v ) plots a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 29 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the triangle vertices # # string filename: a filename for the plot. # import matplotlib.pyplot as plt plt.clf ( ) I = [ 0, 1, 2, 0 ]; plt.plot ( v[I,0], v[I,1], 'b-', linewidth = 3 ) plt.plot ( v[:,0], v[:,1], 'r.', markersize = 30 ) plt.grid ( True ) plt.axis ( 'equal' ) plt.title ( filename ) plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( block = False ) plt.close ( ) return def triangle_point_random ( v, n = 1 ): #*****************************************************************************80 # ## triangle_point_random returns random points in a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # integer n: the number of points to generate. # # Output: # # real p[n,2]: the random points. # import numpy as np p = np.zeros ( [ n, 2 ] ) for i in range ( 0, n ): r = np.random.random ( 2 ) r1 = np.min ( r ) r2 = np.max ( r ) xi = [ r1, r2 - r1, 1.0 - r2 ] p[i] = np.dot ( xi, v ) return p def triangle_point_random_test ( ): #*****************************************************************************80 # ## triangle_point_random_test() tests triangle_point_random(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 05 July 2022 # # Author: # # John Burkardt # import matplotlib.pyplot as plt import numpy as np print ( '' ) print ( 'triangle_point_random_test():' ) print ( ' triangle_point_random() selects random points from a triangle.' ) v = equilateral_triangle ( ) n = 200 p = triangle_point_random ( v, n ) filename = 'triangle_point_random_test.png' plt.clf ( ) I = [ 0, 1, 2, 0 ] plt.plot ( v[I,0], v[I,1], 'b-', linewidth = 3 ) plt.plot ( p[:,0], p[:,1], 'r.', markersize = 15 ) plt.grid ( True ) plt.axis ( 'equal' ) plt.title ( filename ) plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( block = False ) plt.close ( ) return def triangle_sides ( v ): #*****************************************************************************80 # ## triangle_sides() returns side lengths of a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # # Input: # # real v[3,2]: the vertices. # # Output: # # real s[3]: the side lengths. # import numpy as np s = np.zeros ( 3 ) for i in range ( 0, 3 ): s[i] = np.linalg.norm ( v[(i+1)%3] - v[i] ) return s def triangle_test ( v, title, filename ): #*****************************************************************************80 # ## triangle_test() tests a triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # import numpy as np print ( '' ) print ( 'triangle_test:' ) print ( ' Measurements for ' + title ) print ( ' Vertices:' ) print ( v ) triangle_plot ( v, filename ) s = triangle_sides ( v ) print ( ' Sides:' ) print ( s ) p = triangle_perimeter ( v ) print ( ' Perimeter' ) print ( p ) c = triangle_centroid ( v ) print ( ' Centroid' ) print ( c ) area = triangle_area_euclid ( v ) print ( ' Area (Euclid)' ) print ( area ) area = triangle_area_hero ( v ) print ( ' Area (Hero)' ) print ( area ) area = triangle_area_cross_product ( v ) print ( ' Area (Cross product)' ) print ( area ) a = triangle_angles ( v ) print ( ' Angles in radians:' ) print ( a ) print ( ' Angles in degrees:' ) print ( a * 180.0 / np.pi ) return v def triangles_test ( ): #*****************************************************************************80 # ## triangles_test() tests triangles(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 30 June 2022 # # Author: # # John Burkardt # import numpy as np import platform print ( '' ) print ( 'triangles_test():' ) print ( ' Python version: ' + platform.python_version ( ) ) print ( ' triangles() does triangle computations.' ) v = equilateral_triangle ( ) triangle_test ( v, 'equilateral triangle', 'equilateral.png' ) v = obtuse_triangle ( ) triangle_test ( v, 'obtuse triangle', 'obtuse.png' ) v = reference_triangle ( ) triangle_test ( v, 'reference triangle', 'reference.png' ) left_right_test ( ) triangle_contains_test ( ) triangle_integral_centroid_test ( ) triangle_integral_estimate_test ( ) triangle_integral_strang_test ( ) triangle_integral_vertex_test ( ) triangle_point_random_test ( ) # # Terminate. # print ( '' ) print ( 'triangles_test():' ) print ( ' Normal end of execution.' ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 21 August 2019 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) triangles_test ( ) timestamp ( )