Creating (scientific) plots via gnuplot

User projects written in or related to FreeBASIC.
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Postby dodicat » Jun 06, 2019 11:04

Thanks badidea.
I'll give it a try with Julia fractals later, I will have to incorporate a function somehow, but it should be interesting.
How did you load opengl into your Linux system (graphics card driver update perhaps?)
I don't have it.
I tried yum install opengl, which, if it had worked, I would have bought today's lottery ticket.
jj2007
What do you mean "message loop"?
I am using fbgfx to get the opengl driver to run opengl.
example: Screenres 1024,768,32,,2 (The 2 at the end calls the opengl driver)
message loop for a direct winapi application you mean??
jj2007
Posts: 1319
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Postby jj2007 » Jun 06, 2019 13:24

dodicat wrote:jj2007
What do you mean "message loop"?

In a typical C application, it looks like this:

Code: Select all

  while (GetMessage(&msg, NULL, 0, 0)) {
   TranslateMessage(&msg);
   DispatchMessage(&msg);
   }

GetMessage waits for an event. If there is no event (i.e. the user did not hit any key, did not move the mouse, etc), GetMessage simple waits, and cpu usage for the application is exactly 0%. This is the normal behaviour of a Windows application - it does nothing unless you ask for it. Since there are typically about 50 processes up and running, it would be quite a nuisance if their cpu usage was different from 0%.

Google for petzold getmessage to see more detail.
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Postby dodicat » Jun 06, 2019 16:03

Yea jj2007 thanks, that is for direct winapi, (createwindowex(...)) e'.t.c. with a standard winproc to use TranslateMessage( @msg ).
But I have used the freebasic fbgfx screen for opengl.
badidea
Posts: 1782
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Creating (scientific) plots via gnuplot

Postby badidea » Jun 06, 2019 18:11

dodicat wrote:How did you load opengl into your Linux system (graphics card driver update perhaps?)
I don't have it. I tried yum install opengl, which, if it had worked, I would have bought today's lottery ticket.

What does fbc say when you try to compile? ld: cannot find -lGL?
Here it only works with 64-bit fbc. Library name seems to be "libGL"

Code: Select all

badidea@laptop:/media/badidea/data$ ldd dodifractal
   linux-vdso.so.1 =>  (0x00007ffdc83b7000)
   libGL.so.1 => /usr/lib/x86_64-linux-gnu/mesa/libGL.so.1 (0x00007fed03dbd000) <-------------
   libX11.so.6 => /usr/lib/x86_64-linux-gnu/libX11.so.6 (0x00007fed03a83000)
   libXext.so.6 => /usr/lib/x86_64-linux-gnu/libXext.so.6 (0x00007fed03871000)
   libXpm.so.4 => /usr/lib/x86_64-linux-gnu/libXpm.so.4 (0x00007fed0365f000)
   libXrandr.so.2 => /usr/lib/x86_64-linux-gnu/libXrandr.so.2 (0x00007fed03454000)
   libXrender.so.1 => /usr/lib/x86_64-linux-gnu/libXrender.so.1 (0x00007fed0324a000)
   libtinfo.so.5 => /lib/x86_64-linux-gnu/libtinfo.so.5 (0x00007fed03021000)
   libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007fed02d18000)
   libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007fed02b14000)
   libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007fed028f7000)
   libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007fed0252d000)
   libz.so.1 => /lib/x86_64-linux-gnu/libz.so.1 (0x00007fed02313000)
   libexpat.so.1 => /lib/x86_64-linux-gnu/libexpat.so.1 (0x00007fed020ea000)
   libxcb-dri3.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-dri3.so.0 (0x00007fed01ee7000)
   libxcb-present.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-present.so.0 (0x00007fed01ce4000)
   libxcb-sync.so.1 => /usr/lib/x86_64-linux-gnu/libxcb-sync.so.1 (0x00007fed01add000)
   libxshmfence.so.1 => /usr/lib/x86_64-linux-gnu/libxshmfence.so.1 (0x00007fed018da000)
   libglapi.so.0 => /usr/lib/x86_64-linux-gnu/libglapi.so.0 (0x00007fed016a9000)
   libXdamage.so.1 => /usr/lib/x86_64-linux-gnu/libXdamage.so.1 (0x00007fed014a6000)
   libXfixes.so.3 => /usr/lib/x86_64-linux-gnu/libXfixes.so.3 (0x00007fed012a0000)
   libX11-xcb.so.1 => /usr/lib/x86_64-linux-gnu/libX11-xcb.so.1 (0x00007fed0109e000)
   libxcb-glx.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-glx.so.0 (0x00007fed00e85000)
   libxcb-dri2.so.0 => /usr/lib/x86_64-linux-gnu/libxcb-dri2.so.0 (0x00007fed00c80000)
   libxcb.so.1 => /usr/lib/x86_64-linux-gnu/libxcb.so.1 (0x00007fed00a5e000)
   libXxf86vm.so.1 => /usr/lib/x86_64-linux-gnu/libXxf86vm.so.1 (0x00007fed00858000)
   libdrm.so.2 => /usr/lib/x86_64-linux-gnu/libdrm.so.2 (0x00007fed00646000)
   /lib64/ld-linux-x86-64.so.2 (0x00007fed04031000)
   libXau.so.6 => /usr/lib/x86_64-linux-gnu/libXau.so.6 (0x00007fed00442000)
   libXdmcp.so.6 => /usr/lib/x86_64-linux-gnu/libXdmcp.so.6 (0x00007fed0023c000)

Edit: Maybe it is uses libgl1-mesa-glx because I cannot find libGL in the repository.
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Postby dodicat » Jun 06, 2019 21:58

The mesa-libGL is installed (I checked with yum), but I cannot load it, keep getting cannot find . . .
Never mind, thanks for your tips, I'll try some other time.
Here is the gl 3d plot again, optimised for a greater framerate (using screenevent ).

Code: Select all

 
#include "GL/gl.bi"
#include "fbgfx.bi"

Sub setupgl
    Dim As Integer xres,yres
Screeninfo xres,yres
glDisable (GL_DEPTH_TEST)
glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
glEnable (GL_BLEND)
glEnable (GL_LINE_SMOOTH)
glOrtho 0, xres, yres, 0,-1, 1
glClearColor 0,0,.5,1
End Sub

Sub drawstring(x As Long,y As Long,txt As String,size As Long,c As Ulong,b As Ulong)
     #define GL_RGBA_ 6408
     #define GL_BGRA_ 32993
    Dim As Long xx=128*8,yy=16*size
    Static As Any Ptr i
    Static As gluint texture,s
    If s=0 Then glGenTextures(1, @texture):
    i=Imagecreate(128*4,16,b):'8  16
    glBindTexture( GL_TEXTURE_2D, texture ):s=1
     Draw String i,(0,0),txt,c
    glTexImage2d( GL_TEXTURE_2D, 0, GL_RGBA_, 128*4,16, 0, GL_BGRA_, GL_UNSIGNED_BYTE, i+32 )
    glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST )
    glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST )
    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL)
    glEnable( GL_TEXTURE_2D )
    glcolor3ub(Cast(Ubyte Ptr,@c)[2],Cast(Ubyte Ptr,@c)[1],Cast(Ubyte Ptr,@c)[0])
    glbegin gl_quads
    glTexCoord2f(0,0)
    glvertex2f(x,y)
    glTexCoord2f(0,1)
    glvertex2f(x,y+yy)
    glTexCoord2f(1,1)
    glvertex2f(x+xx,y+yy)
    glTexCoord2f(1,0)
    glvertex2f(x+xx,y)
    glend
    gldisable( GL_TEXTURE_2D )
    Imagedestroy i
End Sub


Type V3
    As Single x,y,z
End Type

Operator -(v1 As v3,v2 As v3) As v3 'v1-v2
Return Type(v1.x-v2.x,v1.y-v2.y,v1.z-v2.z)
End Operator

Operator ^ (Byref v1 As v3,Byref v2 As v3) As v3 'cross product
Return Type(v1.y*v2.z-v2.y*v1.z,-(v1.x*v2.z-v2.x*v1.z),v1.x*v2.y-v2.x*v1.y)
End Operator

Type float As V3

Type box
    As v3 p(1 To 4)
    As Ulong c    'colour
    As Single z
End Type

Type angle3D             'FLOATS for angles
    As Single sx,sy,sz
    As Single cx,cy,cz
    Declare Static Function construct(As Single,As Single,As Single) As Angle3D
End Type

Declare Function InputFunction(x As Double,y As Double) As Double

Screenres 1024,768,32,,2
Width 1024\8,768\16 'max dos font size
setupgl
'============ globals =============
Const pi=4*Atn(1)
Redim Shared As box b()
Redim Shared As box rot1()
Dim Shared As Angle3D A3d
Dim Shared As V3 CC       'grid centre
Dim Shared As fb.event e
'Dim Shared As Double x,y  'for inputfunction()
Dim Shared As Single MinX
Dim Shared As Single MaxX
Dim Shared As Single MinY
Dim Shared As Single MaxY
MinX=-3:MaxX=3:Miny=-3:MaxY=3
Dim Shared As Integer xres,yres
Screeninfo xres,yres
'================================== functions ================

       Sub QsortZ(array() As box,begin As Long,Finish As Long)
         Dim As Long i=begin,j=finish
          Dim As box x =array(((I+J)\2))
        While I <= J
            While array(I).z > X.z:I+=1:Wend
            While array(J).z < X.z:J-=1:Wend
        If I<=J Then Swap array(I),array(J): I+=1:J-=1
        Wend
            If J >begin Then QsortZ(array(),begin,J)
            If I <Finish Then QsortZ(array(),I,Finish)
        End Sub
       
        Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
            Return   Type (Sin(x),Sin(y),Sin(z), _
            Cos(x),Cos(y),Cos(z))
        End Function
       
        Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=Type(1,1,1)) As V3
            Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
            Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
            (scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
            (scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z)
        End Function
       
        Function perspective(p As V3,eyepoint As V3) As V3
            Dim As Single   w=1+(p.z/eyepoint.z)
            If w=0 Then w=1e-6
            Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
            (p.y-eyepoint.y)/w+eyepoint.y,_
            (p.z-eyepoint.z)/w+eyepoint.z)
        End Function
       
        Function dot(v1 As v3,v2 As v3) As Single
            Dim As Single d1=Sqr(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z)
            Dim As Single d2=Sqr(v2.x*v2.x + v2.y*v2.y + v2.z*v2.z)
            Dim As Single v1x=v1.x/d1,v1y=v1.y/d1,v1z=v1.z/d1 'normalize
            Dim As Single v2x=v2.x/d2,v2y=v2.y/d2,v2z=v2.z/d2 'normalize
            Return v1x*v2x+v1y*v2y+v1z*v2z  'dot product
        End Function
       
        Function map(a As Single,b As Single,x As Single,c As Single,d As Single) As Single
            Return ((d)-(c))*((x)-(a))/((b)-(a))+(c)
        End Function
       
        Function setgrid(sx As Single,bx As Single,sy As Single,by As Single,st As Single,p() As box,fn As Function(x As Double,y As Double=0) As Double) As v3
            ''485 515   335 365
            #define U Ubound(p)
            Redim p(0)
            Dim As Single cx,cy,ctr
            Static As Single q=15
            Var sttx=st*(MaxX-MinX)/(bx-sx)
            Var stty=st*(MaxY-MinY)/(by-sy)
            For y As Single=sy To by+st/2 Step st
                For x As Single=sx To bx+st/2 Step st
                    ctr+=1
                    cx+=x
                    cy+=y
                    Redim Preserve p(1 To U+1)
                    Var lx=map(sx,bx,x,MinX,MaxX)'+500
                    Var ly=map(sy,by,y,MinY,MaxY)'+350
                   
                    'temp adjust to use limits for .z
                    p(u).p(1)=Type<v3>(lx,ly,          fn(p(u).p(1).x,p(u).p(1).y))
                    p(u).p(2)=Type<v3>(lx+sttx,ly,     fn(p(u).p(2).x,p(u).p(2).y))
                    p(u).p(3)=Type<v3>(lx+sttx,ly+stty,fn(p(u).p(3).x,p(u).p(3).y))
                    p(u).p(4)=Type<v3>(lx,ly+stty,     fn(p(u).p(4).x,p(u).p(4).y))
                    're set
                    p(u).p(1).x=x:     p(u).p(1).y=y
                    p(u).p(2).x=x+st:  p(u).p(2).y=y
                    p(u).p(3).x=x+st:  p(u).p(3).y=y+st
                    p(u).p(4).x=x:     p(u).p(4).y=y+st
                    p(u).c=Rgb(x*q, x*q xor y*q,y*q)
                Next x
            Next y
            Return Type(cx/ctr,cy/ctr)'centre
        End Function
       
       
    Sub drawboxes(b() As box)
        Redim As Long a()
        Redim As V3 aa()
        For n As Long=Lbound(b) To Ubound(b)
       
            Var rd=Cast(Ubyte Ptr,@b(n).c)[2]
            Var gr=Cast(Ubyte Ptr,@b(n).c)[1]
            Var bl=Cast(Ubyte Ptr,@b(n).c)[0]
            Dim As v3 screencentre=(xres\2,yres\2)
            Var v1=b(n).p(2)-b(n).p(1)
            Var v2=b(n).p(3)-b(n).p(2)
            Var norm=v1^v2 'cross product
            Var dt=dot(norm,Type(1,.5,0))
            Var f=map(-1,1,dt,.2,1)
             glbegin gl_quads
        glcolor3ub f*rd,f*gr,f*bl
        For m As Long=4 To 1 Step -1
        glvertex2f b(n).p(m).x,b(n).p(m).y
       Next m
       glend
        Next
    End Sub
   
    Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
        Static As Double timervalue,_lastsleeptime,t3,frames
        frames+=1
        If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
        Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
        If sleeptime<1 Then sleeptime=1
        _lastsleeptime=sleeptime
        timervalue=Timer
        Return sleeptime
    End Function
   
    Sub setup(x1 As Single,x2 As Single,y1 As Single,y2 As Single,meshsize As Single)
        CC= setgrid(x1,x2,y1,y2,meshsize,b(),@InputFunction)'create grid, CC is the centre
        Redim rot1(Lbound(b) To Ubound(b))                   'working array
        A3d=angle3D.construct(0,-pi/2,0)
        Var dx=x2-x1,dy=y2-y1
        Var s=20-map(0,30,dx,0,10)
        For n As Long=Lbound(b) To Ubound(b)
            For m As Long=1 To 4
                rot1(n).p(m)=rotate(CC,B(n).p(m),A3D,Type(s,s,s)) 'align boxes horizontally based
                rot1(n).c=B(n).c
                B(n).p(m)=rot1(n).p(m)
            Next m
        Next n
    End Sub
   
    Function display() As Long
        #define resetwheel(w,fl) fl=w
        #define wheel(w,f) w-f
        Static As float ang=(0,-pi/7,pi/2)  'default
        Static As Long fps
        Static As String key
        Static As Long mx,my,mw,mb,rflag
        Static As Single sc=1
       
        Const k=40
        Var f=map(0,40,k,0,.5)
       
        Do
           
            setup(485,485+k,335,335+k,f)
            Getmouse mx,my,mw,mb
             
            If mb=2 Then 'reset
                ang.z=pi/2:ang.y=-pi/7:ang.x=0
                resetwheel(mw,rflag)
            End If
           
            mw=wheel(mw,rflag)
            If mx>0 Then sc=2+(mw/10)'scaler
            key=Inkey
            If key=Chr(255)+"K" Then ang.z-=.05     'left
            If key=Chr(255)+"M" Then ang.z+=.05     'right
            If key=Chr(255)+"P" Then ang.y-=.05     'down
            If key=Chr(255)+"H" Then ang.y+=.05     'up
            If key="q" Then ang.x+=.05     
            If key="w" Then ang.x-=.05
            '=================== do nothing until an event ===========
            If Screenevent(@e) Then
            A3D=Angle3D.construct(ang.x,ang.y,ang.z)      'set the rotate trigs
           
            For n As Long=Lbound(b) To Ubound(b)
                For m As Long=1 To 4
                    rot1(n).p(m) =rotate(CC,B(n).p(m),A3D,Type(sc,sc,sc))
                    rot1(n).p(m) =perspective(rot1(n).p(m),Type(cc.x,cc.y,400*sc))'eyepoint
                    If mb=1 Then rot1(n).p(m).x-=cc.x-mx: rot1(n).p(m).y-=cc.y-my'follow the mouse
                Next m
               
                rot1(n).z=(rot1(n).p(1).z+rot1(n).p(3).z)/2
            Next n
       
            qsortz(rot1(),Lbound(rot1), Ubound(rot1))
           
            glClear(GL_COLOR_BUFFER_BIT)

            drawstring (50,80,"keys q and w to rotate round vertical (y) axis",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            DrawString(50,110,"Use the arrow keys for x and z axis",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            drawstring(50,140, "Mouse wheel to magnify",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            DrawString(50,170,"Right mouse click to reset",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            drawboxes(rot1())
            End If 'screenevent
         DrawString(50,50,"Framerate "&fps,1.5,Rgb(200,200,200),Rgb(0,0,255\2))' always show fps
            Flip
           
            Sleep regulate(80,fps),1
        Loop Until key=Chr(27)
        Return 0
    End Function
   
    Function InputFunction(x As Double,y As Double) As Double
        'set the x/y domains
        MinX=-pi*4
        MaxX=pi*4
        MinY=-pi*4
        MaxY=pi*4
        If MaxX<MinX Then Swap MaxX,MinX
        If MaxY<MinY Then Swap MaxY,MinY
        ' << --------------- INPUT functions ----------->>
        'return  (((2*x)^2+(2*y)^2))^.5 -20
        'return sin(x)+cos(y)
        'Return  (2*Cos(x/3)+.5*Sin(y/4)*(abs(x)^.3-25*abs(x)^.2))
         Return  3*Sin(5*((x/10)^2+(y/10)^2))
    End Function
   
   
    End display()
    Sleep
   
   
   
 
srvaldez
Posts: 2234
Joined: Sep 25, 2005 21:54

Re: Creating (scientific) plots via gnuplot

Postby srvaldez » Jun 06, 2019 22:23

impressive dodicat :-)
however, if you zoom in with the scroll button, there places where the image flickers a lot, like it flips from one size to the other without deciding what size to settle on.
btw, 60 fps on my 10 year ol Mac :-)
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

Re: Creating (scientific) plots via gnuplot

Postby deltarho[1859] » Jun 06, 2019 22:47

@dodicat

Just for the fun of it I compiled and run using srvaldez's fbc 1.06/gcc 9.1 in 32-bit and 64-bit. I am getting no flicker but my machine may have a 'stronger' CPU, holding down the 'q' key only sees the CPU peak at about 5.8%.
jj2007
Posts: 1319
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Postby jj2007 » Jun 07, 2019 2:07

Truly impressive, dodicat!

I tried to debug it using a MessageBox(), but #include "Windows.bi" causes loads of errors. And a simple print "hello" does not output anything; boh...

Here is a version that does not use the cpu excessively. Search for fpsmax, it should be pretty clear what it does.

Code: Select all

 
#include "GL/gl.bi"
#include "fbgfx.bi"

Sub setupgl
    Dim As Integer xres,yres
Screeninfo xres,yres
glDisable (GL_DEPTH_TEST)
glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA)
glEnable (GL_BLEND)
glEnable (GL_LINE_SMOOTH)
glOrtho 0, xres, yres, 0,-1, 1
glClearColor 0,0,.5,1
End Sub

Sub drawstring(x As Long,y As Long,txt As String,size As Long,c As Ulong,b As Ulong)
     #define GL_RGBA_ 6408
     #define GL_BGRA_ 32993
    Dim As Long xx=128*8,yy=16*size
    Static As Any Ptr i
    Static As gluint texture,s
    If s=0 Then glGenTextures(1, @texture):
    i=Imagecreate(128*4,16,b):'8  16
    glBindTexture( GL_TEXTURE_2D, texture ):s=1
     Draw String i,(0,0),txt,c
    glTexImage2d( GL_TEXTURE_2D, 0, GL_RGBA_, 128*4,16, 0, GL_BGRA_, GL_UNSIGNED_BYTE, i+32 )
    glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST )
    glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST )
    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL)
    glEnable( GL_TEXTURE_2D )
    glcolor3ub(Cast(Ubyte Ptr,@c)[2],Cast(Ubyte Ptr,@c)[1],Cast(Ubyte Ptr,@c)[0])
    glbegin gl_quads
    glTexCoord2f(0,0)
    glvertex2f(x,y)
    glTexCoord2f(0,1)
    glvertex2f(x,y+yy)
    glTexCoord2f(1,1)
    glvertex2f(x+xx,y+yy)
    glTexCoord2f(1,0)
    glvertex2f(x+xx,y)
    glend
    gldisable( GL_TEXTURE_2D )
    Imagedestroy i
End Sub


Type V3
    As Single x,y,z
End Type

Operator -(v1 As v3,v2 As v3) As v3 'v1-v2
Return Type(v1.x-v2.x,v1.y-v2.y,v1.z-v2.z)
End Operator

Operator ^ (Byref v1 As v3,Byref v2 As v3) As v3 'cross product
Return Type(v1.y*v2.z-v2.y*v1.z,-(v1.x*v2.z-v2.x*v1.z),v1.x*v2.y-v2.x*v1.y)
End Operator

Type float As V3

Type box
    As v3 p(1 To 4)
    As Ulong c    'colour
    As Single z
End Type

Type angle3D             'FLOATS for angles
    As Single sx,sy,sz
    As Single cx,cy,cz
    Declare Static Function construct(As Single,As Single,As Single) As Angle3D
End Type

Declare Function InputFunction(x As Double,y As Double) As Double

Screenres 1024,768,32,,2
Width 1024\8,768\16 'max dos font size
setupgl
'============ globals =============
Const pi=4*Atn(1)
Redim Shared As box b()
Redim Shared As box rot1()
Dim Shared As Angle3D A3d
Dim Shared As V3 CC       'grid centre
Dim Shared As fb.event e
'Dim Shared As Double x,y  'for inputfunction()
Dim Shared As Single MinX
Dim Shared As Single MaxX
Dim Shared As Single MinY
Dim Shared As Single MaxY
MinX=-3:MaxX=3:Miny=-3:MaxY=3
Dim Shared As Integer xres,yres
Screeninfo xres,yres
'================================== functions ================

       Sub QsortZ(array() As box,begin As Long,Finish As Long)
         Dim As Long i=begin,j=finish
          Dim As box x =array(((I+J)\2))
        While I <= J
            While array(I).z > X.z:I+=1:Wend
            While array(J).z < X.z:J-=1:Wend
        If I<=J Then Swap array(I),array(J): I+=1:J-=1
        Wend
            If J >begin Then QsortZ(array(),begin,J)
            If I <Finish Then QsortZ(array(),I,Finish)
        End Sub
       
        Function Angle3D.construct(x As Single,y As Single,z As Single) As Angle3D
            Return   Type (Sin(x),Sin(y),Sin(z), _
            Cos(x),Cos(y),Cos(z))
        End Function
       
        Function Rotate(c As V3,p As V3,a As Angle3D,scale As float=Type(1,1,1)) As V3
            Dim As Single dx=p.x-c.x,dy=p.y-c.y,dz=p.z-c.z
            Return Type<V3>((scale.x)*((a.cy*a.cz)*dx+(-a.cx*a.sz+a.sx*a.sy*a.cz)*dy+(a.sx*a.sz+a.cx*a.sy*a.cz)*dz)+c.x,_
            (scale.y)*((a.cy*a.sz)*dx+(a.cx*a.cz+a.sx*a.sy*a.sz)*dy+(-a.sx*a.cz+a.cx*a.sy*a.sz)*dz)+c.y,_
            (scale.z)*((-a.sy)*dx+(a.sx*a.cy)*dy+(a.cx*a.cy)*dz)+c.z)
        End Function
       
        Function perspective(p As V3,eyepoint As V3) As V3
            Dim As Single   w=1+(p.z/eyepoint.z)
            If w=0 Then w=1e-6
            Return Type<V3>((p.x-eyepoint.x)/w+eyepoint.x,_
            (p.y-eyepoint.y)/w+eyepoint.y,_
            (p.z-eyepoint.z)/w+eyepoint.z)
        End Function
       
        Function dot(v1 As v3,v2 As v3) As Single
            Dim As Single d1=Sqr(v1.x*v1.x + v1.y*v1.y + v1.z*v1.z)
            Dim As Single d2=Sqr(v2.x*v2.x + v2.y*v2.y + v2.z*v2.z)
            Dim As Single v1x=v1.x/d1,v1y=v1.y/d1,v1z=v1.z/d1 'normalize
            Dim As Single v2x=v2.x/d2,v2y=v2.y/d2,v2z=v2.z/d2 'normalize
            Return v1x*v2x+v1y*v2y+v1z*v2z  'dot product
        End Function
       
        Function map(a As Single,b As Single,x As Single,c As Single,d As Single) As Single
            Return ((d)-(c))*((x)-(a))/((b)-(a))+(c)
        End Function
       
        Function setgrid(sx As Single,bx As Single,sy As Single,by As Single,st As Single,p() As box,fn As Function(x As Double,y As Double=0) As Double) As v3
            ''485 515   335 365
            #define U Ubound(p)
            Redim p(0)
            Dim As Single cx,cy,ctr
            Static As Single q=15
            Var sttx=st*(MaxX-MinX)/(bx-sx)
            Var stty=st*(MaxY-MinY)/(by-sy)
            For y As Single=sy To by+st/2 Step st
                For x As Single=sx To bx+st/2 Step st
                    ctr+=1
                    cx+=x
                    cy+=y
                    Redim Preserve p(1 To U+1)
                    Var lx=map(sx,bx,x,MinX,MaxX)'+500
                    Var ly=map(sy,by,y,MinY,MaxY)'+350
                   
                    'temp adjust to use limits for .z
                    p(u).p(1)=Type<v3>(lx,ly,          fn(p(u).p(1).x,p(u).p(1).y))
                    p(u).p(2)=Type<v3>(lx+sttx,ly,     fn(p(u).p(2).x,p(u).p(2).y))
                    p(u).p(3)=Type<v3>(lx+sttx,ly+stty,fn(p(u).p(3).x,p(u).p(3).y))
                    p(u).p(4)=Type<v3>(lx,ly+stty,     fn(p(u).p(4).x,p(u).p(4).y))
                    're set
                    p(u).p(1).x=x:     p(u).p(1).y=y
                    p(u).p(2).x=x+st:  p(u).p(2).y=y
                    p(u).p(3).x=x+st:  p(u).p(3).y=y+st
                    p(u).p(4).x=x:     p(u).p(4).y=y+st
                    p(u).c=Rgb(x*q, x*q xor y*q,y*q)
                Next x
            Next y
            Return Type(cx/ctr,cy/ctr)'centre
        End Function
       
       
    Sub drawboxes(b() As box)
        Redim As Long a()
        Redim As V3 aa()
        For n As Long=Lbound(b) To Ubound(b)
       
            Var rd=Cast(Ubyte Ptr,@b(n).c)[2]
            Var gr=Cast(Ubyte Ptr,@b(n).c)[1]
            Var bl=Cast(Ubyte Ptr,@b(n).c)[0]
            Dim As v3 screencentre=(xres\2,yres\2)
            Var v1=b(n).p(2)-b(n).p(1)
            Var v2=b(n).p(3)-b(n).p(2)
            Var norm=v1^v2 'cross product
            Var dt=dot(norm,Type(1,.5,0))
            Var f=map(-1,1,dt,.2,1)
             glbegin gl_quads
        glcolor3ub f*rd,f*gr,f*bl
        For m As Long=4 To 1 Step -1
        glvertex2f b(n).p(m).x,b(n).p(m).y
       Next m
       glend
        Next
    End Sub
   
    Function Regulate(Byval MyFps As Long,Byref fps As Long) As Long
        Static As Double timervalue,_lastsleeptime,t3,frames
        frames+=1
        If (Timer-t3)>=1 Then t3=Timer:fps=frames:frames=0
        Var sleeptime=_lastsleeptime+((1/myfps)-Timer+timervalue)*1000
        If sleeptime<1 Then sleeptime=1
        _lastsleeptime=sleeptime
        timervalue=Timer
        Return sleeptime
    End Function
   
    Sub setup(x1 As Single,x2 As Single,y1 As Single,y2 As Single,meshsize As Single)
        CC= setgrid(x1,x2,y1,y2,meshsize,b(),@InputFunction)'create grid, CC is the centre
        Redim rot1(Lbound(b) To Ubound(b))                   'working array
        A3d=angle3D.construct(0,-pi/2,0)
        Var dx=x2-x1,dy=y2-y1
        Var s=20-map(0,30,dx,0,10)
        For n As Long=Lbound(b) To Ubound(b)
            For m As Long=1 To 4
                rot1(n).p(m)=rotate(CC,B(n).p(m),A3D,Type(s,s,s)) 'align boxes horizontally based
                rot1(n).c=B(n).c
                B(n).p(m)=rot1(n).p(m)
            Next m
        Next n
    End Sub
   
    Function display() As Long
        #define resetwheel(w,fl) fl=w
        #define wheel(w,f) w-f
        Static As float ang=(0,-pi/7,pi/2)  'default
        Static As Long fps, fpsmax=50
        Static As String key
        Static As Long mx,my,mw,mb,rflag
        Static As Single sc=1
       
        Const k=40
        Var f=map(0,40,k,0,.5)
       
        Do
           
            setup(485,485+k,335,335+k,f)
            Getmouse mx,my,mw,mb
             
            If mb=2 Then 'reset
                ang.z=pi/2:ang.y=-pi/7:ang.x=0
                resetwheel(mw,rflag)
            End If
           
            mw=wheel(mw,rflag)
            If mx>0 Then sc=2+(mw/10)'scaler
            key=Inkey
            If key=Chr(255)+"K" Then ang.z-=.05:fpsmax=99     'left
            If key=Chr(255)+"M" Then ang.z+=.05:fpsmax=99     'right
            If key=Chr(255)+"P" Then ang.y-=.05:fpsmax=99     'down
            If key=Chr(255)+"H" Then ang.y+=.05:fpsmax=99     'up
            If key="q" Then ang.x+=.05:fpsmax=99     
            If key="w" Then ang.x-=.05:fpsmax=99
            '=================== do nothing until an event ===========
            If Screenevent(@e) Then
            A3D=Angle3D.construct(ang.x,ang.y,ang.z)      'set the rotate trigs
           
            For n As Long=Lbound(b) To Ubound(b)
                For m As Long=1 To 4
                    rot1(n).p(m) =rotate(CC,B(n).p(m),A3D,Type(sc,sc,sc))
                    rot1(n).p(m) =perspective(rot1(n).p(m),Type(cc.x,cc.y,400*sc))'eyepoint
                    If mb=1 Then rot1(n).p(m).x-=cc.x-mx: rot1(n).p(m).y-=cc.y-my'follow the mouse
                Next m
               
                rot1(n).z=(rot1(n).p(1).z+rot1(n).p(3).z)/2
            Next n
       
            qsortz(rot1(),Lbound(rot1), Ubound(rot1))
           
            glClear(GL_COLOR_BUFFER_BIT)

            drawstring (50,80,"keys q and w to rotate round vertical (y) axis",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            DrawString(50,110,"Use the arrow keys for x and z axis",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            drawstring(50,140, "Mouse wheel to magnify",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            DrawString(50,170,"Right mouse click to reset",1.5,Rgb(200,200,200),Rgb(0,0,255\2))
            drawboxes(rot1())
            End If 'screenevent
         DrawString(50,50,"Framerate "&fps,1.5,Rgb(200,200,200),Rgb(0,0,255\2))' always show fps
            Flip
            if fpsmax>1 Then fpsmax=fpsmax-1
            Sleep regulate(fpsmax,fps),1
        Loop Until key=Chr(27)
        Return 0
    End Function
   
    Function InputFunction(x As Double,y As Double) As Double
        'set the x/y domains
        MinX=-pi*4
        MaxX=pi*4
        MinY=-pi*4
        MaxY=pi*4
        If MaxX<MinX Then Swap MaxX,MinX
        If MaxY<MinY Then Swap MaxY,MinY
        ' << --------------- INPUT functions ----------->>
        'return  (((2*x)^2+(2*y)^2))^.5 -20
        'return sin(x)+cos(y)
        'Return  (2*Cos(x/3)+.5*Sin(y/4)*(abs(x)^.3-25*abs(x)^.2))
         Return  3*Sin(5*((x/10)^2+(y/10)^2))
    End Function
   
   
    End display()
    Sleep
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Postby dodicat » Jun 07, 2019 12:03

jj2007
The word float has conflict with windows.
You can #undef float or change float in the code to another name (_float maybe), only about 4 instances to change.
That single conflict results in a pile of errors ending in too many errors - exiting.
Including windows checks thousands of lines of code during compile.
I even had problems with lastsleeptime in my regulator. (It is used in winnt.bi), who would have thought?

Thanks for testing srvaldez and deltarho[].
You guys are way into fb mac and gcc 9(.1)
I have no problem zooming in here, with 64 bit -Wc -03 I maintain 80 fps while no screenevent.
(using the official build)

Return to “Projects”

Who is online

Users browsing this forum: No registered users and 1 guest