Creating (scientific) plots via gnuplot

User projects written in or related to FreeBASIC.
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Post by jj2007 »

marpon wrote:in fact it is not mysterious at all, you simply do not have the same configuration...

the old dislin.bi did not used namespace , it declared only the sub/functions and sometimes with alias for color , name ... to avoid duplicates
marpon wrote:to use the dislin static lib
download the package for c/c++ 32 or 64 mingw version
I feel sorry for FreeBasic. You really don't realise the problem, right? There is a package, somewhere, and its authors decide out of the blue to change the architecture. There is an old version, somewhere, there is a new version, elsewhere - you do NOT provide a link, you just write "download the package", and then you continue with loads of instructions that certainly make sense for you and your machine...

And all this is spread all over the place, now in a thread about the almost completely unrelated gnuplot, in a few weeks it will be completely outdated because somebody on $%#@ or elsewhere has posted a much better version.

Don't take it personally, marpon. You are a good programmer, a friendly guy, and so are most people in this forum. But FreeBasic is doomed. "Versions" of packages, of IDEs, of manuals, of Wikis, of the language itself all over the place. Code examples that work only if environment variables (dinosaurs of the DOS age) are set "correctly". New "features" that make it resemble more and more that big ugly mess that C++ has become. There are two guiding principles here: 1. Trial and Error, 2. If you are stuck, open a thread on the FB forum. Why would anybody choose a language like FreeBasic? Because it has "BASIC" in its name?
marpon
Posts: 342
Joined: Dec 28, 2012 13:31
Location: Paris - France

Re: Creating (scientific) plots via gnuplot

Post by marpon »

hi jj2007
Don't take it personally, marpon.
not at all, but please, don't you think your are over reacting?

the official distribution for dislin is here https://www.mps.mpg.de/dislin/distributions
that direction was already given on previous posts of that topic
by Carlos Herrera » May 20, 2019 16:40
Go to the dislin page (https://www.mps.mpg.de/dislin/win-64-bit). You will find there the official distribution for Freebasic.
dodicat was also saying what configuration he is using
by dodicat » May 19, 2019 21:56
Here is dislin for 32 bit win.
https://www.mediafire.com/file/ptg394tg ... n.zip/file
Forum as the name says is to discuss , and discussion is speaking and listening
so please calm down

and next time (smile) better not pollute the topic gnuplot with dislin... for me too.
Sorry badidea
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Post by jj2007 »

marpon wrote:the official distribution for dislin is here https://www.mps.mpg.de/dislin/distributions
that direction was already given on previous posts of that topic
Right, it's where you find

Code: Select all

dl_11_fb.zip
FreeBASIC 9534 KB 15-Mar-2019
This is what I am using, this is the one that refuses to work with the latest fashionable namespace syntax...
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Post by dodicat »

The thing is gnuplot is huge (~85 MB in windows) and can be accesses via shell. But dislin is one dll file (1.387 MB) with excellent plotting capabilities, and can be used totally by a freebasic.bi binding file.
(Or the .a file by marpon also via a .bi)
Off topic yes, I suppose.
Anyways, I have no dislin whatsoever in my windows 64 bit fb 1.06, and the dislin.bi in my 32 bit distribution has no mention of namespace.
There are a few different dowloadable fb 1.06.
I am using the .zip win packages which are not depending on minGW for gcc.
I feel off topic again, so I should post a dislin plot to counterbalance, but I am taking the dogs for walkies instead.
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Creating (scientific) plots via gnuplot

Post by badidea »

I read that in gnuplot 5.2 arrays have been added. This way I could create one script file with the data included.
E.g., script file:

Code: Select all

...
set title 'U/I plot'
...
array I[10] = [120, 125, 122, 119, 134, 138, 141, 159, 162, 178]
array U1[10] = [250, 210, 180, 155, 140, 125, 110, 100, 90, 80]
array U2[10] = [255, 215, 185, 160, 145, 130, 115, 105, 95, 85]
...
plot \
I : U1 title 'U1 [V]' with linespoints, \
I : U2 title 'U2 [V]' with linespoints
...
Unfortunately, I have gnuplot 5.0 installed. So I cannot test that currently.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: Creating (scientific) plots via gnuplot

Post by srvaldez »

@badidea
yes, plotting arrays looks interesting

Code: Select all

dim as string plot_command
const max=20

plot_command=chr(34)+"array A["+str(max)+"]=["
for i as long=1 to max
	plot_command+=str(i)
	if i<max then
		plot_command+=", "
	end if
next
plot_command+="];"

plot_command+="array B["+str(max)+"]=["
for i as long=1 to max
	plot_command+=str(log(i))
	if i<max then
		plot_command+=", "
	end if
next
plot_command+="];"
plot_command+="plot sample [i=1:"+str(max)+"] '+' using (A[i]):(B[i]) with linespoints"
plot_command+=chr(34)

shell("gnuplot -p -e " + plot_command)
marcov
Posts: 3454
Joined: Jun 16, 2005 9:45
Location: Netherlands
Contact:

Re: Creating (scientific) plots via gnuplot

Post by marcov »

Carlos Herrera wrote: If an old and buggy fbgfx form is replaced by the gnuplot engine, we will have something
that matches TAChart from Lazarus. After all, Maxima, Octave, nextnano, ... use gnuplot engine,
so its not a shame.
TAChart doesn't use gnuplot or anything external.
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Creating (scientific) plots via gnuplot

Post by badidea »

Slowly making progress on learning gnuplot. Getting some nice graphs now:
Image
Dashed lines is data fitted by gnuplot.
Code cleanup needed before I can post that...

I also found MathGL. I wonder if this could be used from freeBASIC.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Post by dodicat »

Mathgl has tons of big dll's, they all look very c.

I tried to do a Julia fractal with gnu plot, but after about 200 points (an array) it choked and said command line (script) was too long.
Here is an experimental 3D surface plot of functions in x and y.
The input function is right at the bottom.
The framerate is poor (especially on Linux) but gentle keying does the job.
I have omitted the grid lines but have a pattern and some pseudo shading to show up the surface.

Code: Select all


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

'============ 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 Double df,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 fill_polygon(a() As Long, Byval c As Ulong)
            'translation of a c snippet by Angad
            'source of c code: http://code-heaven.blogspot.it/2009/10/simple-c-program-for-scan-line-polygon.html
            Dim As Long      i, j, k, dy, dx, x, y, temp
            Dim As Long      xi(0 To Ubound(a, 1))
            Dim As Single    slope(0 To Ubound(a, 1))
            
            'join first and last vertex
            a(Ubound(a, 1), 0) = a(0, 0)
            a(Ubound(a, 1), 1) = a(0, 1)
            
            For i = 0 To Ubound(a, 1) - 1
                
                dy = a(i+1, 1) - a(i, 1)
                dx = a(i+1, 0) - a(i, 0)
                
                If (dy = 0) Then slope(i) = 1.0
                If (dx = 0) Then slope(i) = 0.0
                
                If (dy <> 0) Andalso (dx <> 0) Then slope(i) = dx / dy
            Next i
            
            For y = 0 To yres - 1
                k = 0
                ' using FB's short-cut operators (which C doesn't have!)
                For i = 0 To Ubound(a, 1) - 1
                    If (a(i, 1) <= y Andalso a(i+1, 1) > y) Orelse _
                    (a(i, 1) > y Andalso a(i+1, 1) <= y) Then
                    xi(k) = Clng(a(i, 0) + slope(i) * (y - a(i, 1)))
                    k += 1
                End If
            Next i
            
            For j = 0 To k - 2
                'Arrange x-intersections in order
                For i = 0 To k - 2
                    If (xi(i) > xi(i + 1)) Then
                        temp = xi(i)
                        xi(i) = xi(i + 1)
                        xi(i + 1) = temp
                    End If
                Next i
            Next j
            'line filling
            For i = 0 To k - 2 Step 2
                Line (xi(i), y)-(xi(i + 1) + 1, y), c
            Next i
        Next y
    End Sub
    
    Sub convert(p() As v3,a() As Long)
        Redim  a((Ubound(p)-Lbound(p))+1,1)
        Dim As Long ctr
        For n As Long=Lbound(p) To Ubound(p)
            a(ctr,0)=p(n).x
            a(ctr,1)=p(n).y
            ctr+=1
        Next n
    End Sub
    
    Sub drawboxes(b() As box)
        Redim As Long a()
        For n As Long=Lbound(b) To Ubound(b)
            convert(b(n).p(),a())
            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)
            fill_polygon(a(),Rgb(f*rd,f*gr,f*bl))
        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
        Screenset 1,0
        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     
            
            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))
            
            Cls
            Draw String(50,50),"Framerate "&fps
            draw string (50,80),"keys q and w to rotate round vertical (y) axis" 
            Draw String(50,110),"Use the arrow keys for x and z axis"
            draw string(50,140), "Mouse wheel to magnify"
            Draw String(50,170),"Right mouse click to reset"
            drawboxes(rot1())
            Flip
            
            Sleep regulate(80,fps),1
        Loop Until key=Chr(27)
        return 0
    End function
    
    end display()
    Sleep
    
    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
        
        'return sin(x)*cos(y)*3  'egg box
        Return  (20*Cos(x/3)+.5*Sin(y/4)*Log(x)*(x^.3-15*x^.2))/12 ' << --------------- INPUT function -----------
    End Function
    
   
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Post by jj2007 »

That looks fantastic, dodicat!

I expected loads of compiler errors but nope, it works out of the box. One minor point: cpu usage goes up 100% for one core, so apparently there is no proper message loop. Inserting a sleep(1000) between Do and Setup cools it down a little bit, but that is not really the right way to do it, of course.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Post by dodicat »

Thanks jj2007.
To ease the cpu workload I have tried opengl rendering (very easy to convert from fbgfx)
It runs a bit better than fbgfx on this box.

Code: Select all

#include "gl/gl.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
setupgl

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 Double df,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     
            
            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,50,"Framerate "&fps,1.5,rgb(200,200,200),rgb(0,0,255\2))
            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())
            Flip
            
            Sleep regulate(80,fps),1
        Loop Until key=Chr(27)
        return 0
    End function
    
    end display()
    Sleep
    
    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
        
        'return sin(x)*cos(y)*3  'egg box
        Return  (20*Cos(x/3)+.5*Sin(y/4)*Log(x)*(x^.3-15*x^.2))/12 ' << --------------- INPUT function -----------
    End Function
    
   
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Creating (scientific) plots via gnuplot

Post by jj2007 »

dodicat wrote:To ease the cpu workload I have tried opengl rendering (very easy to convert from fbgfx)
It runs a bit better than fbgfx on this box.
Indeed a lot faster, but to get cpu usage down to an acceptable level, you need something like Sleep regulate(5,fps),1 on a Core i5. There is no proper message loop. If there was one, cpu usage would be 0% unless user presses a key.
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Creating (scientific) plots via gnuplot

Post by badidea »

Not sure what this 'dodifractal' is doing in my gnuplot topic, but it looks cool and the opengl version runs very smooth on linux was well.
I only had to change the include to "GL/gl.bi"
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Creating (scientific) plots via gnuplot

Post by dodicat »

badidea
Regarding the fractals.
I was wondering if gnu plot could plot pixels (dots) each with a different colour.
srvaldez demonstrated arrays for gnu plot.
So I obtained three arrays x(),y() and c() which is 32 bit colour.
But before I could figure out the commands gnu plot threw an error with "command line too long" at about > 500 array points.
So I didn't continue.
Although I have been off topic, I have still been testing out win gnu plot.
Only, as yet, I have not produced anything worth posting.
Here are the three fractal arrays which I tried to send to gnu plot via script (text file).

Code: Select all

Type complex
    As Single re,im
End Type
Operator +(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re+n2.re,n1.im+n2.im)
End Operator

Operator *(n1 As complex,n2 As complex) As complex
Return Type<complex>(n1.re*n2.re - n1.im*n2.im,n1.im*n2.re + n1.re*n2.im)
End Operator

Sub Julia (Scaler As Single,x() As long,y() as long,c() as ulong)
    #define absolute2(p) p.re*p.re+p.im*p.im
    Dim As Integer xres ,yres,i,ctr
    Screeninfo xres,yres
    Redim x(5*xres*yres),y(5*xres*yres),c(5*xres*yres)
    Dim Z As Complex
    Dim  As Complex initial=Type(.138341,.649857)
    Dim As Single ConstX = -.090010013,ConstY = -.00002514951
    Dim As Long x1,y1
    Dim As Single StartX = -1/Scaler+ConstX, EndX = 1/Scaler+ConstX
    Dim As Single  StartY = -1/Scaler+ConstY,EndY=1/Scaler+ConstY
    Dim As Single IncrementReal=.006245/(Scaler*8),IncrementImag=.01/(Scaler*8)
    For real As Single  = StartX To EndX Step IncrementReal
        For imag As Single = StartY To EndY Step IncrementImag
            ' If Inkey <> "" Then End
            z=Type(real,imag)
            i=0
            Do
                i+=1 
                i=i and 255
                z=z*z+initial  '(Julia iterator)
            Loop Until absolute2(z)>5
            x1=(real-ConstX)*xres*Scaler
            If x1>=0 Then Exit For
            x1=Abs(x1)
            y1=Abs((imag-ConstY)*yres*Scaler)
            If Point(x1,y1)=0 Then
                x(ctr)=x1
                y(ctr)=y1
                c(ctr)=i
                ctr+=1
            End If
            Pset(x1,y1),7
        Next imag
    Next real
    Redim Preserve x(ctr)
     Redim Preserve y(ctr)
      Redim Preserve c(ctr)
End Sub

Sub set(clr() as ulong,lim as long=255)
   ReDim As ULong tmp(101*101)
for i as integer = 0 to lim
   ScreenRes 100,100,8,,-1
   Dim As Any Ptr im=ImageCreate(100,100,,8)
   Paint im,(5,5),i
   BSave("tmp.bmp",im)
   ScreenRes 100,100,32,,-1
   BLoad("tmp.bmp",@tmp(0))
    clr(i)=tmp(50)
   ImageDestroy im:im=0
Next i
Kill "tmp.bmp"
End Sub

sub Make32bitcolor(c()as ulong,clr() as ulong) 'change 8 bit colour to 32 bit colour
    for n as long=lbound(c) to ubound(c)
        var m=(clr(c(n)))
        c(n)=m
    next n
    end sub

Sub CreateJulia(wide As Long,high As Long,num As Single,x() As long,y() as long,c() as long)
    If Screenres(wide,high,8,,-1)=0 Then Julia num,x(),y(),c()
    Screen 0, 0, 0, &h80000000
End Sub

Sub ShowJulia(x() as long,y() as long,c() as ulong)
    Print "done, array sizes =  ";Ubound(x)
    Print "Press a key"
    Sleep
    Screen 20,32
    Windowtitle "Julia points:"
    For n As Long=Lbound(x) To Ubound(x)
        Pset(x(n),y(n)),c(n)
    Next 
    Sleep
    Screen 0, 0, 0, &h80000000
End Sub




print "Please wait, creating arrays . . ."
redim as long x(),y()
redim as ulong c()
dim as ulong clr(255)
set(clr())'get the tranfer array, 8 to 32 bits
createjulia(800,600,200,x(),y(),c())
make32bitcolor(c(),clr())
showjulia x(),y(),c()

sleep
  
badidea
Posts: 2586
Joined: May 24, 2007 22:10
Location: The Netherlands

Re: Creating (scientific) plots via gnuplot

Post by badidea »

dodicat wrote:I was wondering if gnu plot could plot pixels (dots) each with a different colour.
srvaldez demonstrated arrays for gnu plot.
So I obtained three arrays x(),y() and c() which is 32 bit colour.
But before I could figure out the commands gnu plot threw an error with "command line too long" at about > 500 array points.
So I didn't continue.
Although I have been off topic, I have still been testing out win gnu plot.
Only, as yet, I have not produced anything worth posting.
I have only tried 'line plots' in gnuplot so far. And I write the data to a separate data file, not into the script. I expect that gnuplot can handle more data this way. I even create a separate temporary data file for each plot / line. So my script looks like this now:

Code: Select all

# Gnuplot script file
set title "Source - 8.0 slm H2 - 1300 mm"
set xlabel "Cathode current [A]"
set ylabel "Cathode voltage [A]"
set xrange [80:220]
set yrange [80:240]
set grid
set output "Source - 8.0 slm H2 - 1300 mm.png"
set key outside
set terminal pngcairo size 800,600 font "arial,16" linewidth 3
f0(x) = a0*x**b0;
a0 = 1000
b0 = -1;
fit f0(x) "gpdata00.dat" using 1:2 via a0, b0
f1(x) = a1*x**b1;
a1 = 1000
b1 = -1;
fit f1(x) "gpdata01.dat" using 1:2 via a1, b1
f2(x) = a2*x**b2;
a2 = 1000
b2 = -1;
fit f2(x) "gpdata02.dat" using 1:2 via a2, b2
f3(x) = a3*x**b3;
a3 = 1000
b3 = -1;
fit f3(x) "gpdata03.dat" using 1:2 via a3, b3
f4(x) = a4*x**b4;
a4 = 1000
b4 = -1;
fit f4(x) "gpdata04.dat" using 1:2 via a4, b4
f5(x) = a5*x**b5;
a5 = 1000
b5 = -1;
fit f5(x) "gpdata05.dat" using 1:2 via a5, b5
f6(x) = a6*x**b6;
a6 = 1000
b6 = -1;
fit f6(x) "gpdata06.dat" using 1:2 via a6, b6
f7(x) = a7*x**b7;
a7 = 1000
b7 = -1;
fit f7(x) "gpdata07.dat" using 1:2 via a7, b7
plot \
 "gpdata00.dat" using 1:2 title "B: 0.1 T" with points pointtype 7 dashtype 3 lc 1, \
 f0(x) with lines dashtype 3 lc 1 title '', \
 "gpdata01.dat" using 1:2 title "B: 0.2 T" with points pointtype 7 dashtype 3 lc 2, \
 f1(x) with lines dashtype 3 lc 2 title '', \
 "gpdata02.dat" using 1:2 title "B: 0.4 T" with points pointtype 7 dashtype 3 lc 3, \
 f2(x) with lines dashtype 3 lc 3 title '', \
 "gpdata03.dat" using 1:2 title "B: 0.6 T" with points pointtype 7 dashtype 3 lc 4, \
 f3(x) with lines dashtype 3 lc 4 title '', \
 "gpdata04.dat" using 1:2 title "B: 0.8 T" with points pointtype 7 dashtype 3 lc 5, \
 f4(x) with lines dashtype 3 lc 5 title '', \
 "gpdata05.dat" using 1:2 title "B: 1.0 T" with points pointtype 7 dashtype 3 lc 6, \
 f5(x) with lines dashtype 3 lc 6 title '', \
 "gpdata06.dat" using 1:2 title "B: 1.2 T" with points pointtype 7 dashtype 3 lc 7, \
 f6(x) with lines dashtype 3 lc 7 title '', \
 "gpdata07.dat" using 1:2 title "B: 1.4 T" with points pointtype 7 dashtype 3 lc 8, \
 f7(x) with lines dashtype 3 lc 8 title ''
Or just put the fractal formula in the script, see e.g.: http://lowrank.net/gnuplot/fractal/mandelbrot-e.html
But I would chose a different tool for displaying fractals.
Post Reply