Fractal map

User projects written in or related to FreeBASIC.
Luxan
Posts: 46
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Postby Luxan » Jan 08, 2017 21:30

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

Re: Fractal map

Postby Luxan » Jan 09, 2017 13:46

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: 46
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Postby Luxan » Jan 11, 2017 7:47

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: 46
Joined: Feb 18, 2009 12:47
Location: New Zealand

Re: Fractal map

Postby Luxan » Jan 14, 2017 6:27

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
'
' ------------------------------------------------------------
'






Return to “Projects”

Who is online

Users browsing this forum: AWPStar and 2 guests