'-----------------------------------------------------------
' Chua's Chaos Circuit. Double Scroll. Translated to FreeBASIC 
'-----------------------------------------------------------
Dim As Integer sx, sy, sz, colour = 0, counter = 2
Screeninfo sx, sy, sz
sx -= 32
sy -= 16
sz = 4
Screenres sx, sy, sz

Window (-10, -10)-(10, 10)

'-----------------------------------------------
Function g(Byval MN1 As Double,_
    Byval M0 As Double,_
    Byval M1 As Double,_
    Byval BKP1 As Double,_
    Byval BKP2 As Double,_
    Byval VOLT1 As Double) As Double
    Return MN1 * VOLT1 + 0.5 * ( M0 - MN1 ) * ( Abs( VOLT1 + BKP1 ) - Abs( VOLT1 - BKP1 ) )_
    + 0.5 * ( M1 - M0 )  * ( Abs( VOLT1 + BKP2 ) - Abs( VOLT1 - BKP2 ) )
End Function

'-----------------------------------------------
Function tvc1(Byval RE As Double,_
    Byval C1 As Double,_
    Byval MN1 As Double,_
    Byval M0 As Double,_
    Byval M1 As Double,_
    Byval BPK1 As Double,_
    Byval BKP2 As Double,_
    Byval VOLT1 As Double,_
    Byval VOLT2 As Double) As Double
    Return ( 1/(RE*C1)) * ( VOLT2 - VOLT1 - RE * g(MN1, M0, M1, BPK1, BKP2, VOLT1 ) )
End Function

'-----------------------------------------------
Function tvc2(Byval RE As Double,_
    Byval C2 As Double,_
    Byval VOLT1 As Double,_
    Byval VOLT2 As Double,_
    Byval INTL As Double) As Double
    Return ( 1 / ( RE * C2 ) ) * ( VOLT1 - VOLT2 + RE * INTL )
End Function

'-----------------------------------------------
Function til(Byval LE As Double, Byval VOLT2 As Double ) As Double
    Return -VOLT2/LE
End Function

'-----------------------------------------------
' variable parameters of the circuit
Const As Double R = 2000 ' was 1200.  should be 2000
Const As Double RL = 2140.  ' make it  2140.

' fixed parameters in the circuit
Const As Double r1 = 220.
Const As Double r3 = 2200.
Const As Double r4 = 22000.
Const As Double r6 = 3300.
Const As Double c1 = 10.e-9
Const As Double c2 = 100.e-9

' equivalent inductance given by gyrator circuit
Const As Double r8 = 1000.
Const As Double r9 = 1000.
Const As Double r10 = 100.
Const As Double c3 = 100.e-9

Const As Double L = ( r10 * r8 * c3 * RL ) / r9

' parameters for the function of chua diode
Const As Double Esat = 8.3
Const As Double mn1 = 1/r1 + 1/r4
Const As Double m0 =  1/r4 - 1/r3
Const As Double m1 = -1/r6 - 1/r3
Const As Double bp1 = r3 * Esat / ( r3 + r1 )
Const As Double bp2 = r6 * Esat / ( r6 + r4 )

'Print
'Print mn1, m0, m1, bp1, bp2
'Print

'-----------------------------------------------
' USING THE RK4 METHOD FOR SOLVING THE SYSTEM OF DIFFERENTIAL EQUATIONS

' real life time of circuit simulation
Dim As Double TT = 0.001   ' was 0.01

' number of temporal points
Dim As Integer n = 1000      ' was 1000

' size of temporal points
Dim As Double h = TT / n 

' initial time
Dim As Double t

'-----------------------------------------------
' initial conditions
Randomize
'Dim As Double v1 = 3, v2 = 1, il = 1
Dim As Double v1 = 2.5 + Rnd, v2 = 0.5 + Rnd / 5, il = 1

'Print v1, v2, il
Pset( v1, V2 ), colour
Dim As Double min_i = 1e10, max_i = -1e10 

'-----------------------------------------------
Dim As Double kvn11, kvn21, kiln1
Dim As Double kvn12, kvn22, kiln2
Dim As Double kvn13, kvn23, kiln3
Dim As Double kvn14, kvn24, kiln4
Dim As Double  vn1 ,  vn2 ,  iln 

' Print ( 1 / ( R * c1 ) ) * (v2 - v1 - R * g( mn1, m0, m1, bp1, bp2, v1 ) )

'-----------------------------------------------------------

Do
    ' finding slopes with the k's
    '  k1
    kvn11 = tvc1(R, c1, mn1, m0, m1, bp1, bp2, v1, v2)
    kvn21 = tvc2(R, c2, v1, v2, il)
    kiln1 = til(L, v2)
    
    '  k2
    kvn12 = tvc1(R, c1, mn1, m0, m1, bp1, bp2, v1+kvn11*h/2., v2+kvn21*h/2.)
    kvn22 = tvc2(R, c2, v1+kvn11*h/2., v2+kvn21*h/2., il+kiln1*h/2.)
    kiln2 = til(L, v2+kvn21*h/2.)
    
    '  k3
    kvn13 = tvc1(R, c1, mn1, m0, m1, bp1, bp2, v1+kvn12*h/2., v2+kvn22*h/2.)
    kvn23 = tvc2(R, c2, v1+kvn12*h/2., v2+kvn22*h/2., il+kiln2*h/2.)
    kiln3 = til(L, v2+kvn22*h/2.)
    
    '  k4
    kvn14 = tvc1( R, c1, mn1, m0, m1, bp1, bp2, v1+kvn13*h, v2+kvn23*h)
    kvn24 = tvc2( R, c2, v1+kvn13*h, v2+kvn23*h, il+kiln3*h )
    kiln4 = til( L, v2 + kvn23 * h )
    
    ' finding new points
    vn1 = v1 + (1./6.) * (kvn11+2*kvn12+2*kvn13+kvn14) * h
    vn2 = v2 + (1./6.) * (kvn21+2*kvn22+2*kvn23+kvn24) * h
    iln = il + (1./6.) * (kiln1+2*kiln2+2*kiln3+kiln4) * h
    
    ' terminal print of new set of points
    ' Print vn1, vn2, iln
    
    ' Pset( v1, v2 - v1/5 ), colour
    Line -( v1, v2 - v1/5 ), colour
    
    if counter > 0 then    ' change colour after counts
        counter -= 1
    else
        colour = 6
    end if
    ' If min_i > iln Then min_i = iln
    ' If max_i < iln Then max_i = iln
    
    ' next iteration preparation
    v1 = vn1
    v2 = vn2
    il = iln
    t = t + h
Loop Until Len(Inkey) ' While t < TT

'-----------------------------------------------------------
Screen 0
' Print " Min i ="; min_i
' Print " Max i ="; max_i
' Sleep