Fractal map

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Post by Luxan »

123
Last edited by Luxan on Jan 11, 2017 7:40, edited 1 time in total.
Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Post by Luxan »

Oops , forgot to put the previous code in brackets .



file fract03.par

Code: Select all


test {
reset=2004 type=formula formulafile=fractint.frm
formulaname=d1iMandelbrot
corners=-1.599499/1.008761/-0.9269799/1.029215 float=y maxiter=1024
inside=bof60 colors=@blues.map
}

formula , append this to fractint.frm

Code: Select all

d1iMandelbrot(XAXIS) {; Edward Montague (c) 2016
;
; Mandelbrot series = z
; First Derivative of Mandelbrot Series = z1
; Second Derivative of Mandelbrot Series = z2
;
; Differential Equation is , diff(f(z),z) + 3*f(z)*z = sin(z)
; Actually w.r.t c , where c == Pixel.
;

z = Pixel
z1=1
z2 = 0
u = 0
ed = 0:

z2 = 2*z*z2+2*z1^2
z1 = 2*z*z1+1
u = z
z = z*z + Pixel
ed = 3*z*Pixel + z1 - sin(u)
|ed| < 10000 && |z| > 0

}


Now the FreeBASIC code to use with this

Code: Select all


'
'
' (c) Copyright 2017 sciwise@ihug.co.nz
'
' FracD1M.bas
'
'
' ---------------------------------------------------------------------------------
'
' differential equation via mandelbrot series .
'
' ed:diff(f(z),z) + 3*f(z)*z = sin(z) , maxima cas
'
' ed = 3*z*Px + z1 - sin(u) , fractint
'
' ----------------------------------------------------------------------------------
'
'
' To select using the mouse , comment out fg=attractors(a2() ,a3(),Image)
' and uncomment fg = mousey() .
'
' I used Maxima CAS fairly extensively to produce the equations for the
' main iteration loop[s] , this avoids , for now , the use of the double precision
' complex functions that I've written .
'
'
Type dcomplex
x As Double
y As Double
End Type
'
'
Const LEFTBUTTON = 1
Const MIDDLEBUTTON = 4 ' UNUSED IN THIS DEMO
Const RIGHTBUTTON = 2 ' UNUSED IN THIS DEMO
Const SHOWMOUSE = 1
Const HIDEMOUSE = 0
'
'
Declare Function sinh(x As Double) As Double
Declare Function cosh(x As Double) As Double
'
Declare Function dcmul(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcadd(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcsub(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcmulr(v As dcomplex,w As Double) As dcomplex
Declare Function dcaddr(v As dcomplex,w As Double) As dcomplex
Declare Function dcabs(v As dcomplex) As Double
Declare Function dcsin(v As dcomplex) As dcomplex

declare function store(Byref a1 as double,Byref b1 as double,Byref c1 as double,Byref d1 as double) as integer
'
declare function Mandelbrot(a1() as double,Px as dcomplex ,n As Integer) As integer
Declare function d1iMandelbrot(a1() as double,Px as dcomplex ,n As Integer) As Double

declare function mag(x as double,y as double) as double
declare function mousey() as integer
declare function plot2d(a1() as double, steps as integer ) as integer
declare function scanner2(xImage As Any Ptr ,a2() as double,a3() as integer ) as integer
declare function attractors(a2() as double,a3() as integer,Image As Any Ptr ) as integer
'
'
declare function pal() as integer
declare function getputxy(i as integer,j as integer , Image As Any Ptr , flag as integer) as integer
'
declare function transi2x(i as integer,a as double,b as double) as double
declare function transj2y(j as integer,c as double,d as double) as double
'
declare function transx2i(x as double , a as double , b as double) as integer
declare function transy2j(y as double,c as double,d as double) as integer
'
' -----------------------------------------------------------------------------
'
' ScreenRes 820, 690, 8 ' all fractint image files are 8 bit == 256 colours
'
Const W1 = 820, H1 = 690
ScreenRes W1, H1,8
'
dim as double a2()
dim as integer a3()
Dim As Integer twid, tw, th
'
dim Image As Any Ptr = ImageCreate( 21, 21 )
Dim myImage As Any Ptr = ImageCreate( 800, 600 )

Dim fg As double
'
'
' ================ main =======================
'
'
' Load an 800x600 bitmap into an image
'
'
BLoad "fract003.bmp", myImage
'
' -----------------------------------------------------------------------------
'
Width W1\8, H1\16 '' Use 8*16 font
'
twid = Width()
tw = LoWord(twid): th = HiWord(twid)
'
fg= scanner2(myImage ,a2(),a3())
'
window screen (0,0)-(W1,H1)
Put (10,10), myImage
'
' Press space bar to step through attractors .
'
fg=attractors(a2() ,a3(),Image)
'fg = mousey()
'
ImageDestroy( Image )
ImageDestroy( myImage )


Sleep
'
end
'
' ========================== fin ===============
'
' functions
'
'
'
Function sinh(x As Double) As Double
'
' Hyperbolic sine function
'
'
Static y As Double

y = (Exp(x) - Exp(-x))/2


Return y
End Function
'
'
'
Function cosh(x As Double) As Double
'
' Hyperbolic cosine function
'
'
Static y As Double

y = (Exp(x) + Exp(-x))/2


Return y
End Function
'
'
'
Function mousey() as integer
'
' Use mouse to select a point from the fractal .
'
'
'
Dim CurrentX As Integer
Dim CurrentY As Integer
Dim MouseButtons As Integer
Dim CanExit As Integer
Dim n As Integer
n=256
Dim as double ax,bx,cy,dy,a1(n,1),x,y
dim as double fg
Dim Px As Dcomplex
'

fg= store(ax , bx ,cy,dy)
'
'
SetMouse 1, 1, SHOWMOUSE

Do

window screen (0,0)-( W1, H1)
view (0,0)-(W1,H1)


GetMouse CurrentX, CurrentY, , MouseButtons

If MouseButtons And LEFTBUTTON Then
'
If (CurrentX >=10 and CurrentX <=810) and (CurrentY >=10 and CurrentY <= 610 ) then

' circle(CurrentX,CurrentY), 10,20


x = ax+(bx-ax)*(CurrentX-10)/800
y = cy+(dy-cy)*(CurrentY-10)/600
Px.x = x
Px.y = y

fg = d1iMandelbrot(a1() ,Px ,n )

fg = plot2d(a1() , n )
color 150,0
locate 42,2
print " ";
locate 42,2
print "px ";x;
locate 43,2
print " ";
locate 43,2
print "py ";y;


End If
'
End If

Loop While Inkey$ = ""

return 0

end function
'
' ------------------------------------------------------------------------------
'
function mag(x as double ,y as double) as double
'
' sqrt(x*x+y*y)
'
static as double w

w=x*x+y*y
if (w>0) then w=sqr(w)

return (w)
'
end function
'
' ------------------------------------------------------------------------------
'
function plot2d(a1() as double, steps as integer ) as integer
'
' Plot sequence generated
'
'
static as integer i,j
static as double maxx,maxy,x,y,u,v
static as double minx,miny
'
window (0,1)-(steps,-1)
view (10,612)-(810,688)
line (0,1)-(steps,-1),0,bf
'
i=0
minx = 1000000
miny=1000000
for i=0 to steps-8
x = a1(i,0)
y = a1(i,1)
x = abs(x)
y= abs(y)
if (x < minx) then minx=x
if(y < miny) then miny=y
next i
'
for i=0 to steps/4
a1(i,0) = 0
a1(i,1) = 0
next i
'
i=2
maxx = 0
maxy = 0
for i=0 to steps
x = a1(i,0)
y = a1(i,1)
x = abs(x)
y= abs(y)
if (x>maxx) then maxx=x
if(y>maxy) then maxy=y
next i
'
' color 10,0 ' bright green
color 96,0
locate 42,32
print " ";
locate 42,32
print "max z ";maxx;
locate 42,65
print " ";
locate 42,65
print "min z ";minx;

' color 11,0 ' cyan
color 150,0
locate 43,32
print " ";
locate 43,32
print "max d ";maxy;
locate 43,65
print " ";
locate 43,65
print "min d ";miny;


'
if (maxx=0) then maxx=1
if (maxy=0) then maxy=1
'
x = -a1(0,0)/maxx
y = -a1(0,1)/maxy
j=0
'
for i=1 to steps
u = -a1(i,0)/maxx
v = -a1(i,1)/maxy
line(j,x)-(i,u), 96
line(j,y)-(i,v), 150
j = i
x = u
y = v
next i
'
window screen (0,0)-( W1, H1)
view (0,0)-(W1,H1)
'
return (i)

end function
'
' --------------------------------------------------------------------------------
'
function pal() as integer
'
' Examine palette associated with fractint image .
'
' Choose lower and upper limits for scanner2 function .
'
static as integer i,c
'
line(10,10)-(266,50),0,b
line(10,10)-(266,50),56,b
for i=1 to 255
line(i+10,10)-(i+10,50),i,bf
next i
line(10,10)-(266,50),56,b
'
line(10,70)-(266,110),0,bf
line(10,70)-(266,110),56,b
for i=1 to 255
c = point(i+10,20)
line(i+10,70)-(i+10,110),c,bf
if (c=90) then line(i+10,70)-(i+10,110),12,bf
if (c=102) then line(i+10,70)-(i+10,110),12,bf
next i
line(10,70)-(266,110),56,b
'
return (0)
'
end function
'
' -----------------------------------------------------------------------------
'
function transx2i(x as double , a as double , b as double) as integer
'
' Translate from map coordinate to screen coordinate.
'
static as double i
'
i = (800*x-800*a)/(b-a)
'
return i
'
end function
'
' ------------------------------------------------------------------------------
'
function transy2j(y as double,c as double,d as double) as integer
'
' Translate from map coordinates to screen coordinates.
'
static as integer j

j = (600*y-600*c)/(d-c)
'
return j
'
end function
'
' ----------------------------------------------------------------------------
'
function transi2x(i as integer,a as double,b as double) as double
'
' translate from screen coordinate to map coordinate
'
static as double x
'
x = a+(b-a)*i/800
'
return (x)
'
end function
'
' -----------------------------------------------------------------------------
'
function transj2y(j as integer,c as double,d as double) as double
'
' translate from screen coordinate to map coordinate
'
static as double y
'
y = c+(d-c)*(j)/600
'
return (y)
'
end function
'
' -----------------------------------------------------------------------------
'
function store(Byref a1 as double,Byref b1 as double,Byref c1 as double,Byref d1 as double) as integer
'
' Store coordinates for this fractal
'
' From fract03.par
'
' corners=-1.599499/1.008761/-0.9269799/1.029215
' ordering of corner data .
'
'
' [x, xmin , xmax ]
' [y, ymin , ymax ]
'
'
' Top left corner .
'

a1=-1.599499
c1=-0.9269799

'
' Bottom right corner .
'

b1=1.008761
d1=1.029215

'
return (0)
'
end function
'
' -----------------------------------------------------------------------------
'
function scanner2(xImage As Any Ptr ,a2() as double,a3() as integer ) as integer
'
' scan image , in memory , for stable points ; these
' are coloured white when using the blue color map.
'
' Use lower and upper limits selected from function pal().
'
'
' The dimensions of the image are : 800x600 , n x m
'
'
static as integer i,j,n,m,c1,k
static as double ax,bx,cy,dy,x,y
'
i= store(ax ,bx ,cy ,dy )
'
n=800
m=600
'
k = 0
for j=0 to m
for i=0 to n
c1=point(i,j,xImage)
' if (c1 > 86) and (c1<114) then k=k+1
if (c1 > 86) and (c1<106) then k=k+1

next i
next j
'
redim as double a2(k,1)
redim as integer a3(k,1)
'
k = 0
for j=0 to m
for i=0 to n
c1=point(i,j,xImage)
' if (c1 > 86) and (c1<114) then
if (c1 > 86) and (c1<106) then
a2(k,0)=i
a2(k,1)=j
k=k+1
end if
next i
next j
'
for c1 =0 to k
i=a2(c1,0)
j=a2(c1,1)
a3(c1,0)=i
a3(c1,1)=j
x = ax+((bx-ax)*i)/800
y = cy+((dy-cy)*j)/600
a2(c1,0) = x
a2(c1,1) = y
next c1
'
return (0)
'
'
end function
'
' --------------------------------------------------------------------------------
'
function attractors(a2() as double,a3() as integer,Image As Any Ptr ) as integer
'
' Waveforms from results of scanner2
'
' Note : a2() holds [x,y] values , a1() holds sequence values.
'
static as integer k ,i,j,g,fg
static as integer n
n=256
static as double a1(n,1),x,y
static as dcomplex Px
'
'
k=ubound(a2)

'print"k====";k

for g = 0 to k-1
x = a2(g,0)
y = a2(g,1)
Px.x = x
Px.y = y


fg = d1iMandelbrot(a1() ,Px ,n )
x=a1(fg-1,0)
y=a1(fg-1,1)
for i=j-1 to n
a1(i,0)=x
a1(i,1)=y
next i
fg = plot2d(a1() , n )
'
color 150,0
locate 42,2
print " ";
locate 42,2
print "px ";x;
locate 43,2
print " ";
locate 43,2
print "py ";y;
'
i = a3(g,0)
j = a3(g,1)

fg = getputxy(i ,j , Image , 1 )
sleep ' 800
fg = getputxy(i ,j , Image , 0 )
next g
'
return (k)
'
'
end function
'
' ------------------------------------------------------------------------------
'
function getputxy(i as integer,j as integer , Image As Any Ptr , flag as integer) as integer
'
' Selectively ,
'
' Draw circle around a chosen point .
' Return image to original instance .
'
' i == x
' j == y
'
'
select case flag
case 0
if (i>=0) and (j>=0) and (i<=780) and (j<=580) then
Put (i,j),image,pset
end if
'
case 1
'
if (i>=0) and (j>=0) and (i<=780) and (j<=580) then
Get (i,j)-(i+20,j+20), image
circle(i+10,j+10), 8,96
end if
'

case else

end select
'
'
return 0
'
end function
'
' --------------------------- complex functions -----------------------------
'
Function dcmul(v As dcomplex,w As dcomplex) As dcomplex
'
' complex multiplication , double precision .
'
'
Static dc As dcomplex


dc.x = v.x*w.x - v.y*w.y
dc.y= w.x*v.y + w.y*v.x

Return dc
End Function
'
'
'
Function dcmulr(v As dcomplex,w As Double) As dcomplex
'
' complex number multiplied by a double precision value.
'
'
Static dc As dcomplex


dc.x = v.x*w
dc.y= w*v.y

Return dc
End Function
'
'
'
Function dcadd(v As dcomplex,w As dcomplex) As dcomplex
'
' double complex value added to a double complex .
'
'
Static dc As dcomplex


dc.x = v.x + w.x
dc.y= v.y + w.y

Return dc
End Function
'
'
'
Function dcsub(v As dcomplex,w As dcomplex) As dcomplex
'
' double complex subtracted from a double complex .
'
' v - w
'
Static dc As dcomplex


dc.x = v.x - w.x
dc.y= v.y - w.y

Return dc
End Function
'
'
'
Function dcaddr(v As dcomplex,w As Double) As dcomplex
'
' double precision value added to a complex .
'
'
Static dc As dcomplex


dc.x = v.x + w
dc.y= v.y

Return dc
End Function
'
'
'
Function dcabs(v As dcomplex) As Double
'
' Absolute value of a complex number .
'
'
Static dc As Double

dc = v.x*v.x + v.y*v.y

if ( dc > 0 ) then dc = sqr(dc)

Return dc

End Function
'
'
'
Function dcsin(v As dcomplex) As dcomplex
'
' complex sin function .
'
'
' %i*cos(a)*sinh(b)+sin(a)*cosh(b) . maxima
' sinh(x):=(exp(x)-exp(-x))/2
' cosh(x):=(exp(x)+exp(-x))/2
'
'
Static dc As dcomplex

dc.x = sin(v.x)*cosh(v.y)
dc.y = cos(v.x)*sinh(v.y)

Return dc

End Function
'
' ------------------------------------------------------------------------------
'
function Mandelbrot(a1() as double,Px as dcomplex ,n As Integer) As integer
'
' Mandelbrot series = z
' First Derivative of Mandelbrot Series = z1
' Second Derivative of Mandelbrot Series = z2
'
' Differential Equation is , diff(f(z),z) + 3*f(z)*z = sin(z)
' Actually w.r.t c , where c == Px.
'
Dim z As dcomplex
Dim tc As dcomplex
Dim i As integer
Dim j As integer
z = Px


j=0
for i=0 to n - 1
tc.x=z.x*z.x - z.y*z.y
tc.y=z.y*z.x*2
tc.x = tc.x + Px.x
tc.y = tc.y + Px.y
z.x = tc.x
z.y = tc.y
' z = z*z + Px
a1(i,0) = z.x
a1(i,1) = dcabs(z)
if (dcabs(z) > 4 ) then exit for
j=j+1
next i

Return j

end function
'
' ------------------------------------------------------------------------------
'
function d1iMandelbrot(a1() as double,Px as dcomplex ,n As Integer) As Double
'
' Mandelbrot series = z
' First Derivative of Mandelbrot Series = z1
' Second Derivative of Mandelbrot Series = z2
'
' Differential Equation is , diff(f(z),z) + 3*f(z)*z = sin(z)
' Actually w.r.t c , where c == Px.
'
Dim As dcomplex z,z1,z2,u,ed,t,sc
Dim i As integer
Dim j As integer
'
'
for i=0 to n
a1(i,0)=0
a1(i,1)=0
next i
'
z = Px
z1.x=1
z1.y=0
z2.x = 0
z2.y=0
u.x = 0
u.y=0
ed.x = 0
ed.y = 0

j=0
for i=0 to n - 1
t.x = -2*(z . y)*(z2 . y)+2*(z . x)*(z2 . x)-2*(z1 . y^2)+2*(z1 . x^2)
t.y = 2*(z . x)*(z2 . y)+2*(z . y)*(z2 . x)+4*(z1 . x)*(z1 . y)
z2.x = t.x
z2.y = t.y
' z2 = 2*z*z2+2*z1^2
u.x = z.x
u.y = z.y
' u = z
t.x = -2*(z . y)*(z1 . y)+2*(z . x)*(z1 . x)+1
t.y = 2*(z . x)*(z1 . y)+2*(z . y)*(z1 . x)
z1.x = t.x
z1.y = t.y
' z1 = 2*z*z1+1
t.x = -z . y^2+z . x^2+Px . x
t.y = 2*(z . x)*(z . y)+Px . y
z.x = t.x
z.y = t.y
' z = z*z + Px
t.x=-2*(z . y)*(z1 . y)+2*(z . x)*(z1 . x)-sin(z . x)*cosh(z . y)-3*(Px . x)*(z . y^2)-6*(Px . y)*(z . x)*(z . y)+3*(Px . x)*(z . x^2)-3*(Px . y^2)+3*(Px . x^2)+1
t.y = 2*(z . x)*(z1 . y)+2*(z . y)*(z1 . x)-cos(z . x)*sinh(z . y)-3*(Px . y)*(z . y^2)+6*(Px . x)*(z . x)*(z . y)+3*(Px . y)*(z . x^2)+6*(Px . x)*(Px . y)
ed.x = t.x
ed.y = t.y
' ed = 3*z*Px + z1 - sin(u)
a1(i,0) = dcabs(z)
a1(i,1) = dcabs(ed)
if (dcabs(z) > 4 ) then exit for
j=j+1
next i



Return j

end function



Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Post by Luxan »

There's an error in one of the routines , hence I shall re post .

The correct coordinates should now appear in the lower left side of the screen .

I've also set the sleep interval from indefinite to less than a second , now you don't
need to press any key to get a new location that has been determined by the scanner2
routine.


Code: Select all


'
'
'   (c) Copyright 2017 sciwise@ihug.co.nz
'
'  FracD1M.bas
'
'
' ---------------------------------------------------------------------------------
'
'   differential equation via mandelbrot series .
'
'         ed:diff(f(z),z) + 3*f(z)*z = sin(z)    ,   maxima cas
'
'         ed = 3*z*Px  + z1 - sin(u)               ,    fractint
'
' ----------------------------------------------------------------------------------
'
'
'    To select using the mouse , comment out fg=attractors(a2() ,a3(),Image)
' and uncomment fg = mousey() .
'
'      I used Maxima CAS fairly extensively to produce the equations for the
' main iteration loop[s] , this avoids , for now , the use of the double precision 
' complex functions that I've written .
'
'
Type dcomplex
	x As Double
	y As Double
End Type
'
'
Const LEFTBUTTON   = 1
Const MIDDLEBUTTON = 4   ' UNUSED IN THIS DEMO
Const RIGHTBUTTON  = 2   ' UNUSED IN THIS DEMO
Const SHOWMOUSE    = 1
Const HIDEMOUSE    = 0
'
'
Declare  Function sinh(x As Double) As Double
Declare  Function cosh(x As Double) As Double
'
Declare Function dcmul(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcadd(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcsub(v As dcomplex,w As dcomplex) As dcomplex
Declare Function dcmulr(v As dcomplex,w As Double) As dcomplex
Declare Function dcaddr(v As dcomplex,w As Double) As dcomplex
Declare Function dcabs(v As dcomplex) As Double
Declare Function dcsin(v As dcomplex) As dcomplex

declare function store(Byref a1 as double,Byref b1 as double,Byref c1 as double,Byref d1 as double) as integer
'
declare function Mandelbrot(a1() as double,Px as dcomplex ,n As Integer) As integer
Declare function d1iMandelbrot(a1() as double,Px as dcomplex ,n As Integer) As Double

declare function mag(x as double,y as double) as double
declare function mousey() as integer
declare function plot2d(a1() as double, steps as integer ) as integer
declare function scanner2(xImage As Any Ptr ,a2() as double,a3() as integer ) as integer
declare function attractors(a2() as double,a3() as integer,Image As Any Ptr ) as integer
'
'
declare function pal() as integer
declare function getputxy(i as integer,j as integer , Image As Any Ptr , flag as integer) as integer
'
declare function transi2x(i as integer,a as double,b as double)  as double
declare function transj2y(j as integer,c as double,d as double) as double
'
declare function transx2i(x as double , a as double , b as double) as integer
declare function transy2j(y as double,c as double,d as double) as integer
'
' -----------------------------------------------------------------------------
'
' ScreenRes 820, 690, 8 ' all fractint image files are 8 bit == 256 colours
'
Const W1 = 820, H1 = 690
ScreenRes W1, H1,8
'
dim as double a2()
dim as integer a3()
Dim As Integer twid, tw, th
'
dim Image As Any Ptr = ImageCreate( 21, 21 )
Dim myImage As Any Ptr = ImageCreate( 800, 600 )

Dim fg As double
'
'
' ================  main =======================
'
'
'             Load an 800x600 bitmap into an image
'
'
'BLoad "fract003.bmp", myImage

'BLoad "fract004.bmp", myImage

'BLoad "fract005.bmp", myImage

BLoad "fract006.bmp", myImage
'
' -----------------------------------------------------------------------------
'
Width W1\8, H1\16 '' Use 8*16 font
'
twid = Width()
tw = LoWord(twid): th = HiWord(twid)
'
fg= scanner2(myImage ,a2(),a3())
'
window screen  (0,0)-(W1,H1)
Put (10,10), myImage
'
'  Press space bar to step through attractors .
'
fg=attractors(a2() ,a3(),Image)
'fg = mousey()
'
ImageDestroy( Image )
ImageDestroy( myImage )


Sleep
'
end
'
' ========================== fin ===============
'
'                      functions  
'
'
'
Function sinh(x As Double) As Double
'
'   Hyperbolic sine function
'
'	
Static y As Double

         y = (Exp(x) - Exp(-x))/2

	
	Return   y
End Function
'
'
'
Function cosh(x As Double) As Double
'
'   Hyperbolic cosine function
'
'	
Static y As Double

         y = (Exp(x) + Exp(-x))/2

	
	Return   y
End Function
'
'
'
Function mousey() as integer
	'
	'   Use mouse to select a point from the fractal .
	'
	'
	'
	Dim CurrentX     As Integer
	Dim CurrentY     As Integer
	Dim MouseButtons As Integer
	Dim CanExit      As Integer
	Dim n As Integer
	n=256
	Dim as double  ax,bx,cy,dy,a1(n,1),x,y
	dim as double fg
	Dim Px As Dcomplex
	'

	fg= store(ax , bx ,cy,dy)
	'
	'
	SetMouse 1, 1, SHOWMOUSE

	Do

		window screen (0,0)-( W1, H1)
		view (0,0)-(W1,H1)


		GetMouse CurrentX, CurrentY, , MouseButtons

		If MouseButtons And LEFTBUTTON Then
			'
			If (CurrentX >=10 and CurrentX <=810) and (CurrentY >=10 and CurrentY <= 610 ) then

				'      circle(CurrentX,CurrentY), 10,20


				x = ax+(bx-ax)*(CurrentX-10)/800
				y = cy+(dy-cy)*(CurrentY-10)/600
				Px.x = x
				Px.y = y
				
				fg = d1iMandelbrot(a1() ,Px  ,n ) 
				
				fg = plot2d(a1() , n )
				color 150,0
				locate 42,2
				print "                        ";
				locate 42,2
				print "px  ";x;
				locate 43,2
				print "                        ";
				locate 43,2
				print "py  ";y;


			End If
			'
		End If

	Loop While Inkey$ = ""

	return 0

end function
'
' ------------------------------------------------------------------------------
'
function mag(x as double ,y as double) as double
	'
	'                                        sqrt(x*x+y*y)
	'
	static as double w

	w=x*x+y*y
	if (w>0) then w=sqr(w)

	return (w)
	'
end function
'
' ------------------------------------------------------------------------------
'
function plot2d(a1() as double, steps as integer ) as integer
	'
	'   Plot sequence generated 
	'
	'
	static as integer i,j
	static as double maxx,maxy,x,y,u,v
	static as double  minx,miny
	'
	window  (0,1)-(steps,-1)
	view (10,612)-(810,688)
	line (0,1)-(steps,-1),0,bf
'
	i=0
	minx = 1000000
	miny=1000000
	for i=0 to steps-8
		x = a1(i,0)
		y = a1(i,1)
		x = abs(x)
		y= abs(y)
		if (x < minx) then minx=x
		if(y < miny) then miny=y	
	next i
'
for i=0 to steps/4
		 a1(i,0) = 0
		 a1(i,1) = 0
next i
'
	i=2
	maxx = 0
	maxy = 0
	for i=0 to steps
		x = a1(i,0)
		y = a1(i,1)
		x = abs(x)
		y= abs(y)
		if (x>maxx) then maxx=x
		if(y>maxy) then maxy=y
	next i
	'
'	color 10,0   '   bright  green
	color  96,0   
	locate 42,32
	print "                        ";
	locate 42,32
	print "max z  ";maxx;
	locate 42,65
	print "                        ";
	locate 42,65
	print "min z  ";minx;

'	color 11,0        '  cyan
	color 150,0       
	locate 43,32
	print "                        ";
	locate 43,32
	print "max d  ";maxy;
	locate 43,65
	print "                        ";
	locate 43,65
	print "min d  ";miny;
	

	'
	if (maxx=0) then maxx=1
	if (maxy=0) then maxy=1
	'
	x = -a1(0,0)/maxx
	y = -a1(0,1)/maxy
	j=0
	'
	for i=1 to steps
		u = -a1(i,0)/maxx
		v = -a1(i,1)/maxy
		line(j,x)-(i,u), 96
		line(j,y)-(i,v), 150
		j = i
		x = u
		y = v
	next i
	'
	window screen (0,0)-( W1, H1)
	view (0,0)-(W1,H1)
	'
	return (i)

end function
'
' --------------------------------------------------------------------------------
'
function pal() as integer
	'
	'  Examine palette associated with fractint image .
	'
	'  Choose lower and upper limits for scanner2 function .
	'
	static as integer i,c
	'
	line(10,10)-(266,50),0,b
	line(10,10)-(266,50),56,b
	for i=1 to 255
		line(i+10,10)-(i+10,50),i,bf
	next i
	line(10,10)-(266,50),56,b
	'
	line(10,70)-(266,110),0,bf
	line(10,70)-(266,110),56,b
	for i=1 to 255
		c = point(i+10,20)
		line(i+10,70)-(i+10,110),c,bf
		if (c=90) then line(i+10,70)-(i+10,110),12,bf
		if (c=102) then line(i+10,70)-(i+10,110),12,bf
	next i
	line(10,70)-(266,110),56,b
	'
	return (0)
	'
end function
'
' -----------------------------------------------------------------------------
'
function transx2i(x as double , a as double , b as double) as integer
	'
	'    Translate from map coordinate to screen coordinate.
	'
	static as double i
	'
	i = cint((800*x-800*a)/(b-a))
	'
	return i
	'
end function
'
' ------------------------------------------------------------------------------
'
function transy2j(y as double,c as double,d as double) as integer
	'
	'   Translate from map coordinates to screen coordinates.
	'
	static as integer j

	j = cint((600*y-600*c)/(d-c))
	'
	return j
	'
end function
'
'  ----------------------------------------------------------------------------
'
function transi2x(i as integer,a as double,b as double)  as double
	'
	'  translate from screen coordinate to map coordinate
	'
	static as double x
	'
	x = a+(b-a)*cdbl(i)/800
	'
	return (x)
	'
end function
'
' -----------------------------------------------------------------------------
'
function transj2y(j as integer,c as double,d as double) as double
	'
	'  translate from screen coordinate to map coordinate
	'
	static as double y
	'
	y = c+(d-c)*cdbl(j)/600
	'
	return (y)
	'
end function
'
' -----------------------------------------------------------------------------
'
function store(Byref a1 as double,Byref b1 as double,Byref c1 as double,Byref d1 as double) as integer
	'
	'   Store coordinates for this fractal
	'
	'    From fract03.par
	'
    '  corners=-1.599499/1.008761/-0.9269799/1.029215 
     '
  '     From fract04.par
  ' 
   'corners=-0.37797247/0.19274093/-0.94139064/-0.51335559
   '
    '  ordering  of corner data .
   '
   '   [x, xmin  , xmax  ]
   '   [y, ymin  , ymax  ]
   '
	'
	'   Top left corner .
	'
	a1=-1.599499
	c1=-0.9269799
	'
	'  Bottom right corner .
	'
	b1=1.008761
	d1=1.029215 
'
'  .............................   fract04.par values .
'
	'   Top left corner .
	'
	a1=-0.37797247
	c1=-0.94139064
	'
	'  Bottom right corner .
	'
	b1=0.19274093
	d1=-0.51335559
	'
	 '     From fract05.par
  ' 
  ' corners=-0.15582996/-0.011544468/-0.93384235/-0.82562824 float=y
    ' 
  '   Top left corner .
	'
	a1=-0.15582996
	c1=-0.011544468
	'
	'  Bottom right corner .
	'
	b1=-0.93384235
	d1=-0.82562824
	'
	 '     From fract06.par
  ' 
 '  corners=-0.115379457/-0.0911813897/-0.900322728/-0.882174178 float=y
  
    ' 
  '   Top left corner .
	'
	a1=-0.115379457
	c1=-0.900322728
	'
	'  Bottom right corner .
	'
	b1=-0.0911813897
	d1=-0.882174178
	
	return (0)
	'
end function
'
' -----------------------------------------------------------------------------
'
function scanner2(xImage As Any Ptr ,a2() as double,a3() as integer ) as integer
	'
	'    scan  image , in memory , for stable points ; these
	'  are coloured white when using the blue color map.
	'
	'  Use lower and upper limits selected from function pal().
	'
	'
	'  The dimensions of the image are :  800x600 , n x m
	'
	'
	static as integer i,j,n,m,c1,k
	static as double ax,bx,cy,dy,x,y
	'
	i= store(ax ,bx ,cy ,dy )
	'
	n=800
	m=600
	'
	k = 0
	for j=0 to m
		for i=0 to n
			c1=point(i,j,xImage)
			'     if (c1 > 86) and (c1<114) then   k=k+1
			if (c1 > 86) and (c1<106) then   k=k+1

		next i
	next j
	'
	redim as double a2(k,1)
	redim as integer a3(k,1)
	'
	k = 0
	for j=0 to m
		for i=0 to n
			c1=point(i,j,xImage)
			'      if (c1 > 86) and (c1<114) then
			if (c1 > 86) and (c1<106) then
				a2(k,0)=i
				a2(k,1)=j
				k=k+1
			end if
		next i
	next j
	'
	for c1 =0 to k
		i=a2(c1,0)
		j=a2(c1,1)
		a3(c1,0)=i
		a3(c1,1)=j
		x = ax+((bx-ax)*i)/800
		y = cy+((dy-cy)*j)/600
		a2(c1,0) = x
		a2(c1,1) = y
	next c1
	'
	return (0)
	'
	'
end function
'
' --------------------------------------------------------------------------------
'
function attractors(a2() as double,a3() as integer,Image As Any Ptr ) as integer
	'
	'        Waveforms from results of scanner2
	'
	'   Note :  a2() holds [x,y] values , a1() holds sequence values.
	'
	static as integer k ,i,j,g,fg
	static as integer n
	n=256
	static as double a1(n,1),x,y
	static as dcomplex Px
	'
	'
	k=ubound(a2)

	'print"k====";k

	for g = 0 to k-1
		x = a2(g,0)
		y = a2(g,1)
       Px.x = x
       Px.y = y

		
		fg = d1iMandelbrot(a1() ,Px  ,n ) 
		      x=a1(fg-1,0)
		      y=a1(fg-1,1)
		for i=j-1 to n
		     a1(i,0)=x
		     a1(i,1)=y
		next i
		fg = plot2d(a1() , n )
		'
		color 150,0
		locate 42,2
		print "                        ";
		locate 42,2
		print "px  ";Px.x;
		locate 43,2
		print "                        ";
		locate 43,2
		print "py  ";Px.y;
		'
		i = a3(g,0)
		j = a3(g,1)

		fg = getputxy(i ,j  , Image  , 1 )
		sleep 200
		'sleep  ' 800
		fg = getputxy(i ,j  , Image  , 0 )
	next g
	'
	return (k)
	'
	'
end function
'
' ------------------------------------------------------------------------------
'
function getputxy(i as integer,j as integer , Image As Any Ptr , flag as integer) as integer
	'
	'                   Selectively ,
	'
	'                   Draw circle around a chosen point .
	'                   Return image to original instance .
	'
	'       i == x
	'       j == y
	'
	'
	select case flag
		case 0
			if (i>=0) and (j>=0) and (i<=780) and (j<=580) then
				Put (i,j),image,pset
			end if
			'
		case 1
			'
			if (i>=0) and (j>=0) and (i<=780) and (j<=580) then
				Get (i,j)-(i+20,j+20), image
				circle(i+10,j+10), 8,96
			end if
			'

		case else

	end select
	'
	'
	return 0
	'
end function
'
' --------------------------- complex functions -----------------------------
'
Function dcmul(v As dcomplex,w As dcomplex) As dcomplex
'
'      complex multiplication , double precision .
'
'	
Static  dc As  dcomplex

          
           dc.x = v.x*w.x - v.y*w.y
           dc.y= w.x*v.y + w.y*v.x	
	
	Return  dc
End Function
'
'
'
Function dcmulr(v As dcomplex,w As Double) As dcomplex
'
'      complex number multiplied by a  double precision value.
'
'	
Static  dc As  dcomplex

           
           dc.x = v.x*w
           dc.y= w*v.y 
	
	Return  dc
End Function
'
'
'
Function dcadd(v As dcomplex,w As dcomplex) As dcomplex
'
'      double  complex value added to a double complex   .
'
'	
Static  dc As  dcomplex

           
           dc.x = v.x + w.x
           dc.y= v.y  + w.y
	
	Return  dc
End Function
'
'
'
Function dcsub(v As dcomplex,w As dcomplex) As dcomplex
'
'      double  complex  subtracted from  a double complex   .
'
'                                           v - w
'	
Static  dc As  dcomplex

           
           dc.x = v.x - w.x
           dc.y= v.y  - w.y
	
	Return  dc
End Function
'
'
'
Function dcaddr(v As dcomplex,w As Double) As dcomplex
'
'      double precision value added to a complex   .
'
'	
Static  dc As  dcomplex

           
           dc.x = v.x + w
           dc.y= v.y 
	
	Return  dc
End Function
'
'
'
Function dcabs(v As dcomplex) As Double
'
'        Absolute value of a complex number .
'
'	
Static dc As Double

      dc = v.x*v.x + v.y*v.y
      
      if ( dc > 0 ) then dc = sqr(dc)

     Return dc
	
End Function
'
'
'
Function dcsin(v As dcomplex) As dcomplex
'
'         complex sin function .
'
'	
'            %i*cos(a)*sinh(b)+sin(a)*cosh(b)        .     maxima
'            sinh(x):=(exp(x)-exp(-x))/2
'           cosh(x):=(exp(x)+exp(-x))/2
'
'
Static dc As dcomplex

       dc.x = sin(v.x)*cosh(v.y) 
       dc.y = cos(v.x)*sinh(v.y)

     Return dc
	
End Function
'
' ------------------------------------------------------------------------------
'
function Mandelbrot(a1() as double,Px as dcomplex ,n As Integer) As integer
'
' Mandelbrot series = z
' First Derivative of Mandelbrot Series  = z1
' Second Derivative of Mandelbrot Series  = z2
'
'  Differential Equation is , diff(f(z),z) + 3*f(z)*z = sin(z)
'  Actually w.r.t  c , where c == Px.
'
Dim z As dcomplex 
Dim tc As dcomplex 
Dim i As integer
Dim j As integer
 z = Px 
 
  
  j=0
for i=0 to n - 1  
               tc.x=z.x*z.x - z.y*z.y
               tc.y=z.y*z.x*2
               tc.x = tc.x + Px.x
               tc.y = tc.y + Px.y 
               z.x = tc.x
               z.y = tc.y
'               z =  z*z + Px
        a1(i,0) = z.x
        a1(i,1) = dcabs(z)
   if (dcabs(z) > 4 ) then exit for   
       j=j+1  
next i        
        
  Return j
  
 end function 
'
' ------------------------------------------------------------------------------
'
function d1iMandelbrot(a1() as double,Px as dcomplex ,n As Integer) As Double
'
' Mandelbrot series = z
' First Derivative of Mandelbrot Series  = z1
' Second Derivative of Mandelbrot Series  = z2
'
'  Differential Equation is , diff(f(z),z) + 3*f(z)*z = sin(z)
'  Actually w.r.t  c , where c == Px.
'
Dim As  dcomplex z,z1,z2,u,ed,t,sc
Dim i As integer
Dim j As integer
'
'
for i=0 to n 
     a1(i,0)=0
     a1(i,1)=0
next i
'
 z = Px 
 z1.x=1
 z1.y=0
 z2.x = 0
 z2.y=0
 u.x = 0 
 u.y=0
 ed.x = 0
 ed.y = 0 
  
  j=0
for i=0 to n - 1  
              t.x = -2*(z . y)*(z2 . y)+2*(z . x)*(z2 . x)-2*(z1 . y^2)+2*(z1 . x^2)
              t.y = 2*(z . x)*(z2 . y)+2*(z . y)*(z2 . x)+4*(z1 . x)*(z1 . y)
            z2.x = t.x
            z2.y = t.y
  '           z2 = 2*z*z2+2*z1^2
             u.x = z.x
             u.y = z.y
 '             u = z
             t.x = -2*(z . y)*(z1 . y)+2*(z . x)*(z1 . x)+1
             t.y = 2*(z . x)*(z1 . y)+2*(z . y)*(z1 . x)
           z1.x = t.x
           z1.y = t.y
'            z1 = 2*z*z1+1
             t.x = -z . y^2+z . x^2+Px . x
             t.y = 2*(z . x)*(z . y)+Px . y
             z.x = t.x
             z.y = t.y
'              z =  z*z + Px
t.x=-2*(z . y)*(z1 . y)+2*(z . x)*(z1 . x)-sin(z . x)*cosh(z . y)-3*(Px . x)*(z . y^2)-6*(Px . y)*(z . x)*(z . y)+3*(Px . x)*(z . x^2)-3*(Px . y^2)+3*(Px . x^2)+1
t.y = 2*(z . x)*(z1 . y)+2*(z . y)*(z1 . x)-cos(z . x)*sinh(z . y)-3*(Px . y)*(z . y^2)+6*(Px . x)*(z . x)*(z . y)+3*(Px . y)*(z . x^2)+6*(Px . x)*(Px . y)
          ed.x = t.x
          ed.y = t.y
'           ed = 3*z*Px  + z1 - sin(u)
       a1(i,0) = dcabs(z)
       a1(i,1) = dcabs(ed)
   if (dcabs(z) > 4 ) then exit for   
       j=j+1  
next i        
 
  
  
  Return j
  
 end function 


Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Post by Luxan »

This is a general program for determining if a series is increasing or decreasing.

For intended use with the fractal map scanner , as mentioned in previous posts to this thread.

Code: Select all


'
'
'        signal2.bas
'
'         (c)  copyright 2017  sciwise@ihug.co.nz
'
'
'
'          Determine if series is increasing or decreasing .
'
'
'
declare function f1(x as double) as double
declare function incdec(a1() as double,n as integer) as integer
declare function maxmin(a1() as double,byref max as double,byref min as double,n as integer) as integer
'
' ------------------------------------------------------------
'
Dim As Integer twid, tw, th , fg , n
Const W1 = 820, H1 = 690

dim as integer i
dim as double x,y,u,v

dim as double max,min

'
' ------------------------------------------------------------
'
ScreenRes W1, H1,8
'
Width W1\8, H1\16 '' Use 8*16 font
'
twid = Width()
tw = LoWord(twid): th = HiWord(twid)
'
window screen  (0,0)-(W1,H1)
'
line (0,0)-(W1-1,H1-1),13,b
line (0,320)-(W1-1,320),13
'
' ============ main ==============
'
n=800
dim a1(0 to n) as double

for i=0 to n
      x=cdbl(i)/100
      y=f1(x)
      a1(i)=y
next i
'
fg = maxmin(a1() , max , min ,n ) 
'
if max =0 then max=1
'
i=0
u=cdbl(i)/100
v=(a1(i)/max)*300 +320
'
for i=1 to n
      x=cdbl(i)/100
      y=(a1(i)/max)*300 +320
      line (i-1,v)-(i,y),11  '  cyan
      u=x
      v=y
      pset(i-1,v),14  '   yellow
next i
'
'                Finding trend .
'

print
print "  max = ";max;"       min = ";min

fg = incdec(a1() ,n )
if (fg=0) then print "  decreasing series "
if (fg=1) then print " i ncreasing series "
'  print " fg = ",fg
line (0,0)-(W1-1,H1-1),13,b

sleep

end
'
' ==============================
'

function f1(x as double) as double
'
'      Test function .
'
static y as double
'
      y = sin(x)*exp(-x)

   return y
'
'
end function
'
' ------------------------------------------------------------
'
function incdec(a1() as double,n as integer) as integer
'
'   trend , if  |y(i+1)| < |y(i)| and (|y(i+1)| < lmin) then mini = i+1 : lmin = v
'               if  |y(i+1)| > |y(i)|  and (|y(i+1)| > lmax) then maxi = i+1 : lmax = v
'
'              if maxi < mini then series is decreasing .
'             if maxi  > mini then series is increasing .
'
'
'      index to successive max and min 
'
dim as double maxi,mini
dim as double lmax,lmin
dim as double v,y
static as integer i,flg
flg=0

maxi =0
mini = 0
lmax =0
lmin = 10^20

for i=0 to n-1
     y=a1(i)
     y=abs(y)
     v=a1(i+1)
     v=abs(v)    
     
if ( v > y ) and (v>lmax) then maxi = i+1:lmax=v
if ( v < y ) and (v<lmin) then mini = i+1 : lmin=v
next i
'
if (maxi < mini) then flg=0 '  decreasing
if (maxi > mini) then flg=1 '  increasing
'
     return flg
'
'
end function
'
' ------------------------------------------------------------
'
function maxmin(a1() as double,byref max as double,byref min as double,n as integer) as integer
'
'   Global max  & min
'
dim as double y
dim as integer i ,flg
'
    flg = 0

'
max =0
min = 10^10
for i=0 to 800
     y=a1(i)
     y=abs(y)
if ( y > max ) then max = y
if ( y < min )  then min = y
next i

'
     return flg
'
'
end function
'
' ------------------------------------------------------------
'





Luxan
Posts: 222
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Post by Luxan »

Just revising these files for the Volterra Lottka fractal map.

Haven't investigated zooming into the map as this likely requires some
recoding of frac_eh2.bas .

You'll require the program fractint, xfractint or equivalent.

To generate the image file used by frac_eh2.bas append this code to
the fractint , fractint.frm file :

Code: Select all


comment={received from Ramiro Perez <RPEREZ@EARN.UTPVM1> 18 Aug 93
}
V-Euler{
x=real(pixel),
y=imag(pixel),
h=real(p1)/2:
u=x-x*y,
w=-y+x*y,
c=x+h*(u+u),
d=y+h*(w+w),
x=c,
y=d,
z=x+flip(y),
|z|<=p2
}
V-Heun{
x=real(pixel),
y=imag(pixel),
p=real(p1),
h=imag(p1)/2:
u=x-x*y,
w=-y+x*y,
a=x+p*u,
b=y+p*w,
c=x+h*(u+(a-a*b)),
d=y+h*(w+(-b+a*b)),
x=c,
y=d,
z=x+flip(y),
|z|<=p2
}




Also append this code to the fractint.par file :

Code: Select all


VHeun { ; Volterra-Lottka
reset=2004 type=formula formulafile=fractint.frm formulaname=V-Heun
corners=0.000150015/6.00015/0/4.5
params=0.73899999999999999/0.73899999999999999/64/0 float=y
maxiter=2048 inside=bof60 logmap=yes colors=@blues.map
}


Start fractint, press Esc to get to main menu, then press @ , then
press F6 ; select fractint.par.
From there scroll to VHuen, select that; the iterations start and a
blue shaded image appears. Let this run to completion, if you have
your computers system sound on an audible beep upon completion is heard.

Once complete press s , to save, a popup window lets you know what the
file was saved as; close that window and fractint.

Open the saved file with an image editor, save from there as fig005.bmp
to the directory where frac_eh2.bas resides.

Now running frac_eh2.bas is going to open, display , and use that image
to either automatically find attractors or let you manually explore
the fractal map by pointing and clicking upon a region.
A graph of the relevant waveforms is displayed below the map.

Now for the code, recently tested by the author, and found to function
as expected.

Firstly the include file EHStr.bas :

Code: Select all


'
'   This is  EHStr.bas
'
'   Various functions used with the Euler Huen iterator .
'
'
'
Const LEFTBUTTON   = 1
Const MIDDLEBUTTON = 4   ' UNUSED IN THIS DEMO
Const RIGHTBUTTON  = 2   ' UNUSED IN THIS DEMO
Const SHOWMOUSE    = 1
Const HIDEMOUSE    = 0
'
declare function mag(x as single,y as single) as single
declare function mousey() as integer
declare function plot2d(a1() as single, steps as integer ) as integer  
'declare function scanner2(xImage As Any Ptr ,a2() as single) as integer
'declare function attractors(a2() as single) as integer
'
'declare function pal() as integer
'
declare function transi2x(i as integer,a as single,b as single)  as single
declare function transj2y(j as integer,c as single,d as single) as single
'
' =========================================
'
function mousey() as integer
'
'   Use mouse to select a point from the fractal .
'
'
'
Dim CurrentX     As Integer
Dim CurrentY     As Integer
Dim MouseButtons As Integer
Dim CanExit      As Integer
Dim as single  ax,bx,cy,dy,a1(256,1),x,y
dim as integer fg
'
fg= store(ax , bx ,cy,dy) 
'
'
SetMouse 1, 1, SHOWMOUSE

Do
   GetMouse CurrentX, CurrentY, , MouseButtons
   
   If MouseButtons And LEFTBUTTON Then
'    
    If (CurrentX >=10 and CurrentX <=810) and (CurrentY >=10 and CurrentY <= 610 ) then
                      x = ax+(bx-ax)*(CurrentX-10)/800
                      y = cy+(dy-cy)*(CurrentY-10)/600
                     fg = EH2a(a1() , 12 , 256,x,y ) 
                     fg = plot2d(a1() , 256 )
   End If
'
   End If
     
Loop While Inkey$ = ""

return 0

end function
'
' ------------------------------------------------------------------------------
'
function mag(x as single ,y as single) as single
'
'                                        sqrt(x*x+y*y)
'
static as single w

     w=x*x+y*y
     if (w>0) then w=sqr(w)
     
     return (w)
'
end function
'
 '
' ------------------------------------------------------------------------------
'
 function plot2d(a1() as single, steps as integer ) as integer  
 '
 '   Plot sequence generated from EH2a()
 '
 '
 'window (0,-1)-(steps,1)
' view (10,612)-(810,688)
 static as integer i,j
 static as single maxx,maxy,x,y,u,v
 
 line (0,-1)-(steps,1),0,bf
 
    i=2
   maxx = 0 
   maxy = 0
  for i=0 to steps
        x = a1(i,0)
        y = a1(i,1)
        x = abs(x)
        y= abs(y)
   if (x>maxx) then maxx=x
   if(y>maxy) then maxy=y
  next i
'  
if (maxx=0) then maxx=1
if (maxy=0) then maxy=1
'
        x = a1(0,0)/maxx
        y = a1(0,1)/maxy
        j=0
' 
 for i=1 to steps
        u = a1(i,0)/maxx
        v = a1(i,1)/maxy
        line(j,x)-(i,u), 96
        line(j,y)-(i,v), 53
         j = i
       x = u
       y = v
  next i
'     
    return (i) 

 end function
'
' ------------------------------------------------------------------------------
'
 function transi2x(i as integer,a as single,b as single)  as single
'
'  translate from screen coordinate to map coordinate
'
static as single x
'
                  x = a+(b-a)*i/800
'
     return (x)
'
end function
'
' -----------------------------------------------------------------------------
'
function transj2y(j as integer,c as single,d as single) as single
'
'  translate from screen coordinate to map coordinate
'
static as single y
'
                              y = c+(d-c)*(j)/600
'
                return (y)
'
end function
'


Now the main file frac_eh2.bas :

Code: Select all


'
' -----------------------------------------------------------------------------
'
'    frac_eh2.bas
'
'
'      (c) Edward . Q . Montague . 
'
'                  [ alias ]
'
'        quintin9g@gmail.com
'
'
'  Modified Euler Heun method for Volterra Lottka
'    differential equations .
'
'     x' =  a x - b x y      = f1(x,y)
'     y' = -g y + c x y    = g1(x,y)
'
' In this instance a=b=c=d = 1 , therefore :
'
'     u = x - x*y            = f1(x,y)
'     w = -y + x*y         = g1(x,y)
'
' -----------------------------------------------------------------------------
'
'
'             Load a 800x600 bitmap into an image
'
'   This is from a fractint iteration .
'
' -----------------------------------------------------------------------------
'
declare function f1(x as single,y as single) as single
declare function g1(x as single,y as single) as single
'
declare function store(Byref a1 as single,Byref b1 as single,Byref c1 as single,Byref d1 as single) as integer
'
declare function EH2a(a1() as single,xfinal as single,steps as integer,x as single,y as single) as integer
'
declare function pal() as integer
declare function scanner2(xImage As Any Ptr ,a2() as single) as integer
declare function attractors(a2() as single) as integer
'
' -----------------------------------------------------------------------------
'
#include "EHStr.bas"
'
' -----------------------------------------------------------------------------
'
dim as integer fg
dim as single a2()
'
' -----------------------------------------------------------------------------
'
ScreenRes 820, 690, 8 ' all fractint image files are 8 bit == 256 colours
Dim myImage As Any Ptr = ImageCreate( 800, 600 )  
'
BLoad "fract005.bmp", myImage
'
'fg=pal()
'sleep
'end
 fg=scanner2(myImage  ,a2() ) 
Put (10,10), myImage
'
window (0,-1)-(128,1)
view (10,612)-(810,688) 
'
'
'   Function attractors will automatically scan for regions of attraction ,
' from the 'blue' coloured map  for  the Euler - Huen iterator.
'
'
fg=attractors(a2() )

'   mousey , to manually select from image, using mouse point and click.
'fg = mousey() 

sleep
ImageDestroy( myImage )
'
end
'
' ==============================
'
function EH2a(a1() as single,xfinal as single,steps as integer,x as single,y as single) as integer
'
'
'    Modified Euler Huen method ,
' applied to differential equations defined by f1() and g1() .
'
'  This is similar to a fractint frm file .
'
'
static as integer p2,iter
static as single p,h,u,w,a,b,c,d,z
'
'
for iter=0 to steps
     a1(iter,0)=0
     a1(iter,1)=0
next iter
'
'
    p = 0.739
    h = 0.739/2
  p2 = 256
iter = 0
    u = 0
   w = 0
    a = 0
    b = 0
    c = 0
    d = 0
    z = 0
'
    while(z<=p2 and (iter<256) )       
                     iter = iter+1
                        u = f1(x,y)
                        w = g1(x,y)
                         a = x+p*u
                         b = y+p*w
                         c = x+h*(u+f1(a,b))
                         d = y+h*(w+g1(a,b))
                         x = c
                         y = d
                         a1(iter-1,0)=x 
                         a1(iter-1,1)=y
  '                       z = x + -flip(y)
                         z = mag(x,y) 
                       wend
'
                     return ( iter)
'
'                         
  end function
' 
' --------------------------------------------------------------------------------
'
function pal() as integer
'
'  Examine palette associated with fractint image .
'
'  Choose lower and upper limits for scanner2 function .
'
static as integer i,c
'
line(10,10)-(266,50),0,b
line(10,10)-(266,50),56,b
for i=1 to 255
     line(i+10,10)-(i+10,50),i,bf
next i
line(10,10)-(266,50),56,b
'
line(10,70)-(266,110),0,bf
line(10,70)-(266,110),56,b
for i=1 to 255
      c = point(i+10,20)
     line(i+10,70)-(i+10,110),c,bf
     if (c=86) then line(i+10,70)-(i+10,110),12,bf
     if (c=114) then line(i+10,70)-(i+10,110),12,bf
next i
line(10,70)-(266,110),56,b
'
              return (0)
'
end function
'
' -----------------------------------------------------------------------------
'
function store(Byref a1 as single,Byref b1 as single,Byref c1 as single,Byref d1 as single) as integer
'
'   Store coordinates for this fractal 
'
'
'
'   Top left corner .
'
a1=0
c1=4.5
'
'  Bottom right corner .
'
b1=6
d1=0
'
              return (0)
'
end function
'
' -----------------------------------------------------------------------------
'
function scanner2(xImage As Any Ptr ,a2() as single) as integer
'
'    scan  image , in memory , for stable points ; these
'  are coloured white when using the blue color map. 
'
'  Use lower and upper limits selected from function pal().
'
'
'  The dimensions of the image are :  800x600 , n x m
'
'
static as integer i,j,n,m,c1,k
static as single ax,bx,cy,dy,x,y
'
i= store(ax ,bx ,cy ,dy ) 
'
   n=800
   m=600
' 
      k = 0
  for j=0 to m
    for i=0 to n
        c1=point(i,j,xImage)
      if (c1 > 86) and (c1<114) then   k=k+1
    next i
  next j  
'
redim as single a2(k,1) 
'
      k = 0
  for j=0 to m
    for i=0 to n
        c1=point(i,j,xImage)
      if (c1 > 86) and (c1<114) then   
         a2(k,0)=i
         a2(k,1)=j
          k=k+1
     end if
    next i
  next j  
'
for c1 =0 to k
      i=a2(c1,0)
      j=a2(c1,1)
     x = ax+(bx-ax)*(i)/800
     y = cy+(dy-cy)*(j)/600     
     a2(c1,0) = x
     a2(c1,1) = y
next c1
'
          return (0) 
'
'
end function
'
' --------------------------------------------------------------------------------
'
function attractors(a2() as single) as integer
'
'        Waveforms from results of scanner2
'
'
'
static as integer k ,i,fg
static as single a1(256,1),x,y
'
'
k=ubound(a2)

'print"k====";k

for i = 0 to k-1
     x = a2(i,0)
     y = a2(i,1)
    fg = EH2a(a1() , 12 , 256 ,x,y) 
    fg = plot2d(a1() , 256 )
   sleep 500
next i
'
return (k)
'
'
end function
'
' --------------------------------------------------------------------------------
'
function f1(x as single,y as single) as single
'
'                       x' = x - x*y            = f1(x,y)
'
'
static as single z


                         z = x - x*y
          return (z)

 end function
'
' -----------------------------------------------------------------------------
'
function g1(x as single,y as single) as single
'
'
'                    y' = -y + x*y         = g1(x,y)
'
'
static as single z


                         z =  -y + x*y  
          return (z)

 end function
'
' ------------------------------------------------------------------------------
'

Post Reply