#! /usr/bin/env python3 # def julia_set ( h, w, xl, xr, yb, yt ): #*****************************************************************************80 # ## julia_set() returns points in a Julia set. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2022 # # Author: # # John Burkardt # # Input: # # integer H, W, the height and width of the region in pixels. # # real XL, XR, YB, YT, the left, right, bottom and top limits. # # Output: # # logical J(h,w): true if the corresponding point is in the Julia set. # # real X(h,w), Y(h,w): the coordinates of the points. # import numpy as np # # Create a hxw grid of X and Y coordinates. # x = np.linspace ( xl, xr, w ) y = np.linspace ( yb, yt, h ) X, Y = np.meshgrid ( x, y ) # # Construct a complex copy of X + Yi. # A = X + Y * 1j # # Repeatedly apply the following transformation: # A -> A * A + C # c = - 0.8 + 0.156 * 1j for k in range ( 0, 200 ): A = A**2 + c # # Record the points that didn't diverge. # J = np.abs ( A ) < 1000.0 return J, X, Y def julia_set_test ( ): #*****************************************************************************80 # ## julia_set_test() tests julia_set(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 23 December 2022 # # Author: # # John Burkardt # import matplotlib.pyplot as plt import platform import warnings print ( '' ) print ( 'julia_set_test():' ) print ( ' Python version: ' + platform.python_version ( ) ) print ( ' Compute and plot a Julia set.' ) print ( ' under the transformation A=A^2-0.8+0.156i' ) w = 1000 h = 1000 xl = -1.5 xr = 1.5 yb = -1.5 yt = 1.5 # # Before calling "julia_set()", we suppress warnings, which # report overflow. # with warnings.catch_warnings ( ): warnings.simplefilter ( "ignore" ) J, X, Y = julia_set ( w, h, xl, xr, yb, yt ) # # Use a plot() command directly on the X, Y values. # plt.figure ( ) plt.clf ( ) plt.plot ( X[J], Y[J], 'r.', markersize = 0.25 ) plt.grid ( True ) filename = 'julia_plot.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) # # Use a spy() command on the logical 0/1 J array. # plt.figure ( ) plt.clf ( ) plt.spy ( J, markersize = 0.25, color = 'b' ) filename = 'julia_spy.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) # # Terminate. # print ( '' ) print ( 'julia_set_test():' ) print ( ' Normal end of execution.' ) return def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 August 2019 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) julia_set_test ( ) timestamp ( )