Файл:Комплект базилики Юлия с позвоночником.svg

 /* пакетный файл для Maxima CAS     максима пакет ( "b.mac" ) */ начало:прошедшее_время_выполнения () ; убить ( все ) ; remvalue ( все ) ; /* -------------------------- определения функций ------------------------------*/ f ( z,c ) :=z*z+c ; /* Комплексное квадратичное отображение */ finverseplus ( z,c ) :=sqrt ( zc ) $ finverseminus ( z,c ) :=-sqrt ( zc ) $        /* */ fn ​​( p, z, c ) := если p=0 тогда z иначеесли p=1 тогда f ( z,c ) иначе f ( fn ( p-1, z, c ) , c ) $                / * Стандартный многочлен F_p \, корнями которого являются периодические z-точки периода p и его делители */ F ( p, z, c ) := fn ( p, z, c ) - z $                        /* Функция для вычисления редуцированного многочлена G_p\, корни которого являются периодическими z-точками периода p без его делителей * / G[p,z , c]:= block ( [f:divisors ( p ) , t:1], /* t является временной переменной = произведением Gn для ( делителей p ) отличных от p */ f:delete ( p,f ) , /* удалить p из списка делителей */ if p=1 then return ( F ( p,z,c )) , for i in f do t :t*G[i,z,c], g : F ( p,z,c ) /t, return ( ratsimp ( g )) ) $                                                   GiveRoots ( g ) := block ( [cc], cc:bfallroots ( expand ( %i*g ) =0 ) , cc:map ( rhs,cc ) , /* удалить строку "c=" */ cc:map ( 'float,cc ) , return ( cc ) ) $            /* endcons комплексная точка для вывода в формате для пакета draw */ endconsD ( point,list ) :=endcons ( [realpart ( point ) , imagpart ( point ) ] ,list ) $ consD ( point,list ) : =cons ( [realpart ( point ) , imagpart ( point ) ],list ) $              GiveForwardOrbit ( z0,c,iMax ) := /* вычисляет ( без проверки на выход ) прямую орбиту точки z0 и сохраняет ее в списке для пакета draw */ block ( [z,orbit,temp], z:z0, / * первая точка = критическая точка z :0+0*%i */ orbit:[[realpart ( z ) , imagpart ( z ) ]], for i:1 thru iMax step 1 do ( z:expand ( f ( z,c )) , orbit:endcons ( [realpart ( z ) , imagpart ( z ) ],orbit )) , return ( orbit ) ) $                                               /* Выдает точки обратной орбиты z =repellor */ GiveBackwardOrbit ( c ,repellor,zxMin,zxMax,zyMin,zyMax,iXmax,iYmax ) := block ( hit_limit:4, /* пропорционально количеству деталей и времени отрисовки * / PixelWidth : ( zxMax - zxMin ) / iXmax , PixelHeight : ( zyMax - zyMin ) / iYmax, /* 2D массив пикселей попаданий . Попадание > 0 означает , что точка была на орбите */ array ( Hits,fixnum,iXmax,iYmax ) , /* нет попаданий для начала * / /* выбираем repellor z=repellor в качестве начальной точки */ stack:[repellor], /* сохраняем repellor в стеке */ /* сохраняем первую точку в списке пикселей */ x_y :[ repellor ], /* обратная итерация repellor */ loop, /* pop = взять одну точку из стека */ z:last ( stack ) , stack : delete ( z, stack ) , /* обратная итерация - первый прообраз ( root ) */ z:finverseplus ( z, c ) , /* перевести из мировых координат в экранные */ iX:fix (( realpart ( z ) -zxMin ) /PixelWidth ) , iY:fix (( imagpart ( z ) -zyMin ) /PixelHeight ) , hit:Hits[iX,iY], if hit<hit_limit then ( Hits[iX,iY]:hit+1,                                                                                                                         stack :endcons ( z,stack ) , /* push = добавить z в конец стека списка * / if hit=0 then x_y:endcons ( z,x_y ) ) , /* обратная итерация - второй прообраз ( root ) */ z:-z, /* перевод из мировых координат в экранные , преобразование в целые числа */ iX:fix (( realpart ( z ) -zxMin ) /PixelWidth ) , iY:fix (( imagpart ( z ) -zyMin ) /PixelHeight ) , hit:Hits[iX,iY], if hit < hit_limit then ( Hits[iX,iY]:hit+1, stack:endcons ( z,stack ) , /* push = добавить z в конец стека списка для продолжения итерации */ if hit=0 then x_y :endcons ( z,x_y ) ) , if is ( not emptyp ( stack )) then go ( loop ) , return ( x_y ) / * список пикселей в форме [z1,z2] * / ) $                                                                                         /*----------------------------------*/ Psi_n ( r,t,z_last, Max_R ) := /* */ block ( [iMax:200, iMax2:0], /* ----- прямая итерация 2 точек : z_last и w -------------- * / array ( forward , iMax - 1 ) , / * прямая орбита z_last для сравнения * / forward [ 0 ] : z_last , i : 0, while cabs ( forward [ i] ) <Max_R and i< ( iMax-2 ) do ( /* прямая итерация z в плоскости fc и сохранение ее в прямом массиве * / forward[i+1]:forward[i]*forward[i] + c, /* z*z+c */ /* прямая итерация w в плоскости f0 : w ( n+1 ) : = wn^2 */ r:r*2, /* радиус квадрата = R^2=2^ ( 2*r ) потому что R=2^r */ t:mod ( 2*t,1 ) , /* */ iMax2:iMax2+1, i:i+1 ) , /* вычислить последнюю точку w ; она равна точке z */ R:2^r, /* w:R*exp ( 2*%pi*%i*t ) , z:w, */ array ( backward,iMax-1 ) , backward[iMax2]:rectform ( ev ( R*exp ( 2*%pi*%i*t ))) , /* использовать последнюю точку w как начальную точку для обратной итерации к новой z */ /* ----- точка обратной итерации z=w в плоскости fc --------------*/ для шага i:iMax2 -1                                                                                                                                  thru 1 do ( temp:float ( rectform ( sqrt ( backward[i]-c ))) , /* sqrt ( zc ) */ scalar_product:realpart ( temp ) *realpart ( forward[i-1] ) +imagpart ( temp ) *imagpart ( forward[i-1] ) , if ( 0>scalar_product ) then temp:-temp, /* выбор прообраза */ backward[i-1]:temp ) , return ( backward[0] ) ) $ GiveRay ( t , c ) := block ( [r], /* диапазон для рисования R=2^r ; поскольку r стремится к 0, R стремится к 1 */ rMin:1E-5, /* 1E-4 ; rMin > 0 ; если rMin=0, то программа имеет бесконечный цикл !!!!! */ rMax:2, caution:0.9330329915368074, /* r:r*caution ; это дает меньшее r */ /* верхний предел для итерации */ R_max:300, / * */ zz:[], /* массив для z - точек луча в плоскости fc */ /* некоторые w-точки внешнего луча в плоскости f0 */ r:rMax, while 2^r<R_max do r:2*r, /* найти точку w на луче около бесконечности ( R>=R_max ) в плоскости f0 */ R:2^r, w:rectform ( ev ( R*exp ( 2*%pi*%i*t ))) , z:w, /* около бесконечности z=w */ zz:cons ( z,zz ) , unless r<rMin do ( /* новое меньшее R */ r:r*caution, R:2^r, /* */ w:rectform ( ev                                                                                                                 ( R*exp ( 2*%pi*%i*t ))) , /* */ last_z:z, z:Psi_n ( r,t,last_z,R_max ) , /* z=Psi_n ( w ) */ zz:cons ( z,zz ) ) , return ( zz ) ) $            /* преобразует комплексное число z = x*y*%i в список в формате рисования : [x,y] * / d ( z ) :=[float ( realpart ( z )) , float ( imagpart ( z )) ]$                ToPoints ( myList ) := точки ( map ( d,myList )) $ /* дать список рисования из одной точки*/ ToPoint ( z ) :=points ( [d ( z ) ] ) $      GiveSpine ( zc, Alfa,c, iMax ) :=block ( [s, center, cut], /* первый центр компонента , только один !!! */ center: zc, s:[center], /* первая пара точек разреза */ s: cons ( Alfa,s ) , cut : -Alfa, s:endcons ( cut ,s ) , for i: 1 thru iMax step 1 do ( /* пара центров компонента * / center : finverseminus ( center,c ) , s: cons ( center, s ) , center : - center, s: endcons ( center, s ) , /* пара точек разреза * / cut : finverseminus ( cut,c ) , s: cons ( cut,s ) , cut : - cut , s:endcons ( cut,s ) ) , /* преобразование в формат чертежа и возврат списка */ s: ToPoints ( s ) ) $                                                                                                                  компилировать ( все ) $ /* ----------------------- основной ----------------------------------------------------*/   период:3$  /* внешний угол в поворотах */ /* разрешение пропорционально количеству деталей и времени рисования */ iX_max:1000$ iY_max:1000$ / * определяем z-плоскость ( динамическую ) */ ZxMin:-2.0$ ZxMax:2.0 $ ZyMin :-2.0$ ZyMax :2.0$                                /* присвоить c значение * / c:-1$ /* центр компонента периода 2 * /               /* вычислить фиксированные точки */ Beta:float ( rectform (( 1+sqrt ( 1-4*c )) /2 )) $ /* вычислить отталкивающую фиксированную точку beta */ alfa:float ( rectform (( 1-sqrt ( 1-4*c )) /2 )) $ /* другая фиксированная точка */                   / * вычислить обратную орбиту отталкивающей неподвижной точки */ xy: GiveBackwardOrbit ( c,Beta,ZxMin,ZxMax,ZyMin,ZyMax,iX_max,iY_max ) $            /* вычислить точки луча и сохранить в списке zz */ eRayZero:GiveRay ( 0/1 , c ) $ eRay1o2:GiveRay ( 1/2 , c ) $               позвоночник: GiveSpine ( 0 , alfa,c, 4 ) $    /* время вычислений * / time:fix ( elapsed_run_time () -start ) $       /* рисуем его с помощью пакета draw */ load ( draw ) $             путь: "~/maxima/batch/julia/spine/basilica/" $ / * если пусто , то файл находится в домашнем каталоге */            /* если графический файл пуст ( = 0 байт ), то снова запустите команду draw2d */ draw2d ( terminal = 'svg, file_name = sconcat ( path, "BasilicaSpine" ) , user_preamble= "set size square;set key top right" , title = concat ( " Динамическая плоскость для fc ( z ) = z * z " , string ( c ) ) , dimensions = [ iX_max , iY_max ] , yrange = [ ZyMin , ZyMax ], xrange = [ZxMin,ZyMax], xlabel = "Z.re " , ylabel = "Z.im" , point_type = filled_circle, points_joined =true, point_size = 0.2 , color = red, points_joined =false, color = black, key = "обратная орбита z=beta" , ToPoints ( xy ) , points_joined = true, point_size = 0.2 , цвет = красный, ключ = "внешний луч 0" , ToPoints ( eRayZero ) , ключ = "внешний луч 1/2" , цвет = пурпурный, ToPoints ( eRay1o2 ) , points_joined = true, point_size = 0.8 , цвет = серый, ключ = "позвоночник" , позвоночник,                                                                                                              points_joined = false, цвет = черный, point_size = 1.4 , ключ = "критическая точка z = 0.0" , ToPoint ( 0.0 ) , цвет = красный, point_size = 1.4 , ключ = "отталкивающая фиксированная точка z = бета" , ToPoint ( Beta ) , цвет = пурпурный, ключ = "минус бета" , ToPoint ( -Beta ) , цвет = желтый, ключ = "параболическая фиксированная точка z = альфа" , ToPoint ( alfa ) ) $