3D surface plot

[dodicat]

z=f(x,y) or z=f(x,y,d) where d is distance of a point from the origin.

I think the dishelper.bi has a slightly different location in the latest fb version.
Just set up your scale with the arrows and key in a function.
trig functions and atn (not atan2) are available, and the usual + - / * ^
Abs is also available and sqr log exp.
Most likely some others, I havn't tried the vbscript with any others yet.
I don't think it'll work on Linux.

Code: Select all

'Surface plotter

#define UNICODE
#include Once "disphelper/disphelper.bi"
''#include once "/gd/disphelper/disphelper.bi" 

Dim Shared As IDISPATCH Ptr VBS

dhInitialize(TRUE)
dhToggleExceptions(FALSE)
dhCreateObject "MSScriptControl.ScriptControl",NULL,@VBS
dhPutValue VBS,".Language %s","VBScript"

Sub _END() 
    SAFE_RELEASE(VBS)
    dhUninitialize True
End Sub
Function eval(s As String) As Double
    Dim value As Zstring Ptr
    dhGetValue "%s",@value,VBS,".Eval %s",s
    Return Val(*value)
End Function
Sub Setvariable(s As String,REPLACE_THIS As String,WITHTHIS As Double)' As String
    var WITH_THIS=Str(WITHTHIS)
    var position=Instr(s,REPLACE_THIS)
    While position>0
        s=Mid(s,1,position-1) & WITH_THIS & Mid(s,position+Len(REPLACE_THIS))
        position=Instr(position+Len(WITH_THIS),s,REPLACE_THIS)
    Wend
End Sub
#macro _window(topleftX,topleftY,bottomrightX,bottomrightY) 
minx=topleftX
maxx=bottomrightX
miny=bottomrightY
maxy=topleftY
#endmacro
#macro _axis(colour)
If Sgn(minx)<>Sgn(maxx) Then
    Line(0,(yres-(miny/(miny-maxy))*yres))-(xres,(yres-(miny/(miny-maxy))*yres)),colour 'x axis
    Endif
    If Sgn(miny)<>Sgn(maxy) Then
        Line(((minx/(minx-maxx))*xres),0)-(((minx/(minx-maxx))*xres),yres),colour 'y axis
        Endif
        Draw String(0,yres/2),Str(minx),colour
        Draw String(xres-8-8*(Len(Str(maxx))),yres/2),Str(maxx),colour
        Draw String(xres/2,0),Str(maxy),colour
        Draw String(xres/2,yres-16),Str(miny),colour
        #endmacro
        Sub inspect
            Dim As Integer mx,my,mw
            mx=70:my=230
            mw=2
            Dim As Uinteger array(1 To 24641)
            Dim As Integer count
            For x As Integer=mx-00 To mx+600
                For y As Integer=my-20 To my+20
                    count=count+1
                    array(count)=Point(x,y)
                Next y
            Next x
            count=0
            'draw
            For x As Integer=mx-00 To mx+600
                For y As Integer=my-20 To my+20
                    count=count+1
                    var NewX=mw*(x-mx)+mx
                    var NewY=mw*(y-my)+my 
                    Line(NewX-mw/2,NewY-mw/2)-(NewX+mw/2,NewY+mw/2),array(count),BF
                Next y
            Next x
        End Sub
        'progress
        Type bar
            As Integer x,y,l,d,percent
            As Uinteger col
        End Type
        Dim As Integer percentage
        #define progress(value,lower,upper) 100*(value-lower)/(upper-lower)
        
        Sub progressbar(b As bar)
            Line(b.x+1,b.y+1)-( (b.l*b.percent/100+b.x),b.y+b.d-1),b.col,bf
            Line(b.x,b.y)-(b.x+b.l,b.y+b.d),,b
        End Sub
        Dim  As bar b
        b=Type<bar>(100,300,600,20,0,Rgb(0,0,200))
        
        Dim Shared As Integer xres,yres
        Dim  As Single minx,maxx,miny,maxy,PLOT_GRADE=15000
        
        Screen 20,32
        Screeninfo xres,yres
        Type vector3d
            As Single x,y,z
        End Type
        #define vct type<vector3d>
        #define distance(p1,p2) sqr((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y))
        #macro incircle(cx,cy,radius,x,y)
        (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
        #endmacro
        Dim Shared As vector3d eyepoint
        Dim Shared Rx(1 To 3,1 To 3) As Single
        Dim Shared Ry(1 To 3,1 To 3) As Single
        Dim Shared Rz(1 To 3,1 To 3) As Single
        Dim Shared pivot_vector(1 To 3) As Single
        Dim Shared new_pos(1 To 3) As Single
        Dim Shared temp1(1 To 3) As Single
        Dim Shared temp2(1 To 3) As Single
        Operator + (v1 As vector3d,v2 As vector3d) As vector3d
        Return Type<vector3d>(v1.x+v2.x,v1.y+v2.y,v1.z+v2.z)
        End Operator
        
        Operator -(v1 As vector3d,v2 As vector3d) As vector3d
        Return Type<vector3d>(v1.x-v2.x,v1.y-v2.y,v1.z-v2.z)
        End Operator
        
        Operator * (f As Single,v1 As vector3d) As vector3d
        Return Type<vector3d>(f*v1.x,f*v1.y,f*v1.z)
        End Operator
        Function r(first As Double, last As Double) As Double
            Function = Rnd * (last - first) + first
        End Function
        
        Function apply_perspective(p As vector3d) As vector3d
            Dim As Single   w=(p.z*(-1)/500+1)*.75'300
            Return vct((p.x-eyepoint.x)/w+eyepoint.x,(p.y-eyepoint.y)/w+eyepoint.y,(p.z-eyepoint.z)/w+eyepoint.z)
        End Function
        Sub framecounter
            Static As Double frame,fps
            frame=frame+1
            Static As Double t1,t2
            If frame>=fps Then
                t1 = Timer
                fps = frame/(t1-t2)
                Windowtitle "Frames per second = " & fps
                t2=Timer
                frame=0
            End If
        End Sub
        Function rotatepoint3d(Byval pivot As vector3d,_  
            Byval pt As vector3d,_  
            Byval angle As vector3d,_ 
            Byval dilator As Single=1) As vector3d  
            #macro mv(m1,v,ans)
            For i As Integer=1 To 3
                s=0
                For k As Integer = 1 To 3
                    s=s+m1(i,k)*v(k)
                Next k
                ans(i)=s
            Next i
            #endmacro
            #define cr 0.0174532925199433
            Dim angle_radians As vector3d=Type<vector3d>(cr*angle.x,cr*angle.y,cr*angle.z)
            
            Dim s As Single=Any
            pivot_vector(1)=(pt.x-pivot.x)*dilator
            pivot_vector(2)=(pt.y-pivot.y)*dilator
            pivot_vector(3)=(pt.z-pivot.z)*dilator
            
            'rotat1on matrices about the three axix
            Rx(1,1)=1:Rx(1,2)=0:Rx(1,3)=0
            Rx(2,1)=0:Rx(2,2)=Cos(angle_radians.x):Rx(2,3)=-Sin(angle_radians.x)
            Rx(3,1)=0:Rx(3,2)=Sin(angle_radians.x):Rx(3,3)=Cos(angle_radians.x)
            
            Ry(1,1)=Cos(angle_radians.y):Ry(1,2)=0:Ry(1,3)=Sin(angle_radians.y)
            Ry(2,1)=0:Ry(2,2)=1:Ry(2,3)=0
            Ry(3,1)=-Sin(angle_radians.y):Ry(3,2)=0:Ry(3,3)=Cos(angle_radians.y)
            
            Rz(1,1)=Cos(angle_radians.z):Rz(1,2)=-Sin(angle_radians.z):Rz(1,3)=0
            Rz(2,1)=Sin(angle_radians.z):Rz(2,2)=Cos(angle_radians.z):Rz(2,3)=0
            Rz(3,1)=0:Rz(3,2)=0:Rz(3,3)=1
            
            mv (Rx,pivot_vector,temp1)           
            mv (Ry,temp1,temp2)
            mv (Rz,temp2,new_pos)
            
            new_pos(1)=new_pos(1)+pivot.x
            new_pos(2)=new_pos(2)+pivot.y
            new_pos(3)=new_pos(3)+pivot.z
            
            Return Type<vector3d>(new_pos(1),new_pos(2),new_pos(3))
        End Function
        
        Sub blow(a() As vector3d,mag As Single)
            For z As Integer=1 To Ubound(a)-6
                a(z)=mag*a(z)
            Next z
        End Sub
        
        Sub translate(a() As vector3d,pt As vector3d)
            For z As Integer=1 To Ubound(a)
                a(z)=a(z)+vct(pt.x,pt.y,pt.z)
            Next z
        End Sub
        
        
        Function vertex(piv As vector3d,p1 As vector3d,ang As vector3d,dil As Single,col As Uinteger) As Single
            var _temp1=rotatepoint3d(piv,p1,ang,dil)
            _temp1=apply_perspective(_temp1)
            Pset (_temp1.x,_temp1.y),col
            Return _temp1.z
        End Function
        
        Sub set_perspective(x As Single,y As Single,z As Single,minz As Single,maxz As Single)
            eyepoint=vct(x,y,z)
        End Sub
        #macro combsort(array,comp)
        Scope
            var size=Ub,switch=0,j=0
            Dim As Single void=size
            Do
                void=void/1.3: If void<1 Then void=1
                switch=0
                For i As Integer =1 To size-void
                    j=i+void
                    If comp(i)>comp(j) Then 
                        Swap array(i),array(j): switch=1
                        Swap comp(i),comp(j)
                        Swap col(i),col(j)
                    End If
                Next
            Loop Until  switch =0 And void=1
        End Scope
        #endmacro
        
        Redim Shared As vector3d e(0)
        Dim count As Integer
        Redim Shared As Uinteger col(0)
        Dim As Single funct
        Dim As Single dist,scale=5
        _window(-scale,scale,scale,-scale)
        
        Dim As vector3d pt,cent=Type<vector3d>(0,0,0)
        
        Dim As String formula,worker
        
        Dim As String i,j
        Dim As Integer flag,ub
        Dim As Integer mx,my,mw,counter,mb
        start:
        
        Do
            i=Inkey
            If Left$(i,1)=Chr(08) Then j=Mid(j,1,Len(j)-1)
            Screenlock
            Cls 
            Draw String (.53*xres,.48*yres),"SCALERS"
            Circle(.53*xres,.53*yres),20,Rgb(100,100,100),,,,f
            Draw String(.525*xres,.525*yres),Chr(30)
            Circle(.58*xres,.53*yres),20,Rgb(100,100,100),,,,f
            Draw String(.575*xres,.525*yres),Chr(31)
            Getmouse mx,my,mw,mb
            If incircle(.53*xres,.53*yres,20,mx,my) Then   Circle(.53*xres,.53*yres),20,Rgb(255,255,255)
            If incircle(.58*xres,.53*yres,20,mx,my) Then   Circle(.58*xres,.53*yres),20,Rgb(255,255,255)
            
            If incircle(.53*xres,.53*yres,20,mx,my) And mb=1 Then 
                scale=scale+.1
                _window(-scale,scale,scale,-scale)
            End If
            If incircle(.58*xres,.53*yres,20,mx,my) And mb=1 Then 
                scale=scale-.1
                If scale<.1 Then scale=.1
                _window(-scale,scale,scale,-scale)
            End If
            
            Locate 3,3
            Print "  Example function sin(x)*cos(y)"
            Print "  You can also use d as a variable, which is distance from origin"
            Print "  E.G. (sin(2*d)/(2*d))*5"
            Print "  Use the scalers to adjust the X/Y plane"
            Print "  Enter a function in  x and y (or d) at the arrow "',formula
            If flag=1 Then
                Locate 10,12
                Print "Error in " & formula & " --please redo"
            End If
            Select Case Left$(i,1)
            Case Chr(0) To Chr(254)
                If Left$(i,1)<>Chr(08) Then
                    j=j+Left$(i,1)
                End If
            End Select
            
            Locate 15,5
            Print "----> " & j
            inspect
            _axis(Rgb(0,200,0))
            Screenunlock
            Sleep 1,1
        Loop Until i=Chr(13)
        j=Rtrim(j,Chr(13))
        formula=j
        
        For x As Single=-scale To scale Step 2*scale/150
            For y As Single=-scale To scale Step 2*scale/150
                
                pt=Type<vector3d>(x,y)
                dist=distance(pt,cent)
                count=count+1
                Redim Preserve e(count)
                Redim Preserve col(count)
                col(count)=Rgb(255*(x+scale)/(2*scale),155*(y+scale)/(2*scale),50)
                worker=formula
                setvariable(worker,"x",x)
                setvariable(worker,"y",y)
                setvariable(worker,"d",dist)
                funct=eval(worker)
                e(count)=Type<vector3d>(x,y,funct)
            Next y
            percentage=progress(x,-scale,scale)
            b.percent=percentage
            progressbar(b)
        Next x
        ub=Ubound(e)
        Redim Preserve e(ub+6)
        'axis ends
        e(ub+1)=Type<vector3d>(-.5*(xres/2),0,0)
        e(ub+2)=Type<vector3d>(.5*(xres/2),0,0)
        e(ub+3)=Type<vector3d>(0,-.5*(yres/2),0)
        e(ub+4)=Type<vector3d>(0,.5*(yres/2),0)
        e(ub+5)=Type<vector3d>(0,0,-.5*(yres/2))
        e(ub+6)=Type<vector3d>(0,0,.5*(yres/2))
        Dim As Single dummy
        For z As Integer=1 To ub
            If e(z).z=0  Then dummy=dummy+1
        Next z
        If dummy=ub Then
            Beep
            Erase(e)
            Cls
            flag=1
            j=""
            Goto start
        End If
        flag=0
        blow(e(),20*5/scale)
        translate(e(),vct(xres/2,yres/2,0))
        set_perspective(xres/2,yres/2,0,-100,100)
        Dim As Single dilation
        Dim As vector3d piv,ang
        piv=eyepoint
        Dim As Single zeds(1 To Ub),_mw
        Dim As vector3d axis(6)
        Dim As Uinteger colour
        Dim As Single startdilation=1
        dilation=startdilation
        Dim As Any Pointer im=imagecreate(xres,yres)
        Paint im,(0,0),Rgb(255,255,255)
        
        Do
            framecounter
            i=Inkey
            If i= Chr(255) + "K"  Then ang.x=ang.x+5
            If i= Chr(255) + "M"  Then ang.x=ang.x-5
            If i= Chr(255) + "P"  Then ang.y=ang.y-5
            If i= Chr(255) + "H"  Then ang.y=ang.y+5
            Screenlock
            Cls
            Put(0,0),im
            Draw String (20,20), "Use up/down keys and mouse wheel",Rgb(0,0,0)
            Draw String(20,50), "z= " & formula,Rgb(0,0,0)
            Draw String (.7*xres,.1*yres),"X Y plane = " & scale & " X " & scale,Rgb(0,0,0)
            Getmouse mx,my,mw
            _mw=mw/100
            combsort(e,zeds)
            dilation=startdilation+_mw
            counter=0
            
            'rotate axis
            For z As Integer=ub+1 To ub+6
                counter=counter+1
                axis(counter)=rotatepoint3d(piv,e(z),ang,dilation)
                axis(counter)=apply_perspective(axis(counter))
                
            Next z
            
            Line (axis(1).x,axis(1).y)-(axis(2).x,axis(2).y),Rgb(100,100,100)
            Draw String(axis(1).x,axis(1).y),"-X",Rgb(0,0,0)
            Draw String(axis(2).x,axis(2).y),"+X",Rgb(0,0,0)
            Line (axis(3).x,axis(3).y)-(axis(4).x,axis(4).y),Rgb(100,100,100)
            Draw String(axis(3).x,axis(3).y),"+Y",Rgb(0,0,0)
            Draw String(axis(4).x,axis(4).y),"-Y",Rgb(0,0,0)
            Line (axis(5).x,axis(5).y)-(axis(6).x,axis(6).y),Rgb(100,100,100)
            Draw String(axis(5).x,axis(5).y),"-Z",Rgb(0,0,0)
            Draw String(axis(6).x,axis(6).y),"+Z",Rgb(0,0,0)
            'rotate points and draw
            For z As Integer=1 To Ub
                zeds(z)=vertex(piv,e(z),ang,dilation,col(z))
            Next z
            Screenunlock
            Sleep 1,1
        Loop Until i=Chr(27)
        imagedestroy im
        Sleep
        _END
        

[jdebord]

Nice work dodicat, as usual :)

Did you try the DISLIN library ? (http://www.dislin.de/). There is a dislin.bi file in the FreeBASIC distribution, but I could not find an example.

I have installed the library but not tested with FB yet.

[srvaldez]

the dislin examples are easy enough to translate from c, I played with it some years ago but forgot the details, but in trying the examples again I found that there have been some changes, you will have to translate the new header with swig.

[dodicat]

Thanks for the link jdebord.
Interesting stuff.
I see that you have been at this stuff before srvaldez :
http://www.freebasic.net/forum/viewtopic.php?p=102670

I looked for the lib file, but no luck yet.

@ dafhi

Before I took the dog out for a walk I noticed a reply from you, when I got back with the dog, it had vanished (not the dog)
I was wanting to investigate your optimized combsort.

I don't think that I imagined the reply, otherwise how would I know about the optimized combsort?

[dafhi]

I sometimes see something that needs fixing right after I post.

Anyway, here. Neat project.

Code: Select all

#Macro zInsertionSort()

	i = 1
	For j = 2 To lEnd
		If zeds(i) > zeds(j) Then
			sng_swap = zeds(j)
			uint_swap = cols(j)
			vect_swap = e(j)
			k = i
			zeds(j) = zeds(i)
			e(j)=e(i)
			cols(j)=cols(i)
			For i = j - 2 To 1 Step -1
				If zeds(i) <= sng_swap Then Exit For
				zeds(k) = zeds(i)
				e(k) = e(i)
				cols(k)=cols(i)
				k = i
			Next
			zeds(k) = sng_swap
			e(k) = vect_swap
			cols(k) = uint_swap
		End If
		i = j
	Next
#EndMacro
Private Sub CombSort4(e() As vector3d,zeds() As Single,cols() As UInteger)
Dim As Integer Gap = UBound(zeds)-1, lEnd = Gap
Dim As Integer i,j,k
Dim As Single sng_swap
Dim As UInteger uint_swap
Dim As vector3d vect_swap

	' combsort 4 by dafhi

	Do While Gap > 1
		For i = 1 To lEnd - Gap
			j = i + Gap
			If zeds(i) > zeds(j) Then
				Swap zeds(i), zeds(j)
				Swap e(i),e(j)
				Swap cols(i),cols(j)
			EndIf
		Next
		Gap *= 10
		Gap \= 13
	Loop
	zInsertionSort()
End Sub

[TJF]

Did you try the DISLIN library ? (http://www.dislin.de/). There is a dislin.bi file in the FreeBASIC distribution, but I could not find an example.

I made a package with up-to-date headers and some examples translated from C source. Have a look at

• DISLIN library (10.2 headers and examples)

I didn't test under windows. It'd be nice to get some feedback ;)

[srvaldez]

TJF
I am using the static libraries for mingw and the names are dismg.a and dismg_d.a respectively.

[TJF]

TJF
I am using the static libraries for mingw and the names are dismg.a and dismg_d.a respectively.

Thanks for feedback!

Yes, the naming of this library is a problem, also on LINUX. There the files are named 'dislin-10.2.a' and 'dislin_d-10.2.a' (also without 'lib...' at the beginning), so the linker cannot find them. My solution is using dynamic links, ie named 'libdislin.a', pointing to 'dislin-10.2.a'.

[srvaldez]

hi dodicat
if you want to give dislin a try, I made an import lib from dilnc.dll and put both in a zip file and uploaded it to mediafire http://www.mediafire.com/?9f7uk7v927e9r2w
download the dislin headers and examples from TJF site and change the #inclib entry to #inclib "dislnc"

[dodicat]

hi dodicat
if you want to give dislin a try, I made an import lib from dilnc.dll and put both in a zip file and uploaded it to mediafire http://www.mediafire.com/?9f7uk7v927e9r2w
download the dislin headers and examples from TJF site and change the #inclib entry to #inclib "dislnc"

Thanks srvaldez.

You can also pull the dll's from dl_10_vc.zip (the visual c download from the link given by TJF.
I transferred the dislin.bi from the fb includes and called #inclib "dislnc", then put dislnc.dll and dislnc_d.dll (the quad one) in the same folder.
You of course don't need the dislnc_d.dll for non quad work, in fact it gives errors.
So with TJF's examples and the dislnc.dll (yours or the one pulled from the zip file) and the modified dislin.bi youve got a runner.
So like Cairo, I've just stuck the dll's in a folder without having to install stuff, or actually do any work.
Bit of a cheat really.
The little help text files in the zip package are quite good at explaining calls to the fuctions buried in the binary dlls.

@jdebord
I wonder if Panoramic could use the dll's in this case?
Pointers to arrays, rather than strings.

[rolliebollocks]

Very nice dodicat!

[dodicat]

Very nice dodicat!

Thanks Rollie~
And there's me thinking that you were on holiday, or emigrated.
Have you enrolled yet?

@ others
The disln_d.dll seems to be little different from disln.dll, and not in fact quad precision.
I suppose it was wishful thinking, or just another of my trips into Fool's Paradise.

Fonts are trivial here, and there are loads of routines in dislin.bi to play around with.
It says on the window title "click mb2 to continue", not sure what that is all about.
Here's some trivial fonts.

Code: Select all

#INCLUDE "dislin.bi"

DIM AS ZSTRING*30 sines(11)
sines(0)="Some sines:"
sines(11)="Close window to finish"
for z as integer=1 to 10 
 sines(z)=str(z) & " = " & str(sin(z))
 next z

metafl ("cons")
scrmod ("revers")
setpag("da3p")
disini ()

hwfont()

height(150)

messag (@sines(0), 10, 10)

height(40)
for z as integer=1 to 10 
messag (@sines(z), 20, 200+80*z)
next z
height(100)
messag (@sines(11), 500, 1500)
disfin ()
shell("PAUSE")

[TJF]

The disln_d.dll seems to be little different from disln.dll, and not in fact quad precision.
I suppose it was wishful thinking, or just another of my trips into Fool's Paradise.

My understanding: dislin and dislin_d are not different. Both have the same SUBs/FUNCTIONs, but in dislin_d instead of SINGLE precision DOUBLE is used for variables and arrays.

[jdebord]

The DLL version seems to work well, except for the demo programs Curve.bas, InterpolationMethods.bas and PolarPlots.bas for which I get an error message like this one:

Code: Select all

PolarPlots.o:fake:(.text+0x343): undefined reference to `intrgb'

@dodicat:

1) I think that mb2 means "mouse button 2"

2) It may be difficult to use the DLL with Panoramic. Panoramic cannot pass pointers to arrays, and it uses only double precision reals. If I understand correctly, the DLL is for single precision.

A possibility would be to call the static library (double precision version, dismg_d.a) from TCC and to let Panoramic generate the TCC function. I will try that.

[TJF]

All examples are tested under LINUX. They compile and execute well.

The DLL version seems to work well, except for the demo programs Curve.bas, InterpolationMethods.bas and PolarPlots.bas for which I get an error message like this one:

Code: Select all

PolarPlots.o:fake:(.text+0x343): undefined reference to `intrgb'

Here the compiling went well, but while linking the `intrgb' function couldn't be found in the input library (.a file).

-> Tell the linker where to find the new input library (version 10.2), ie by prepending #LIBPATH "path/to/disln.a" to the code. For a test you may put the .a files into the project folder and use

• #LIBPATH "."