chaos_circuit_high_res.py

#!/usr/bin/python3
#high resolution attractor rendering of simple chaotic RC oscillator with ngspice
#https://www.researchgate.net/publication/277939627_Simple_Two-Transistor_Single-Supply_Resistor-Capacitor_Chaotic_Oscillator

import os
from scipy import  *
import matplotlib.pyplot as plt
import multiprocessing


def do_sim(R4,n):
    #random initial condition
    random.seed(n)
    v3,v4,v5,v6 = random.uniform(low=0,high=2,size=4)
    outfilename = 'out_%i.dat' %n
    circuitfilename = 'circuit_%i' %n
    #prepare the circuit file
    f=open(circuitfilename,'w')
    f.write("""*** RC Chaos Generator ***

* Supply
vcc 1 0 5.0V

* Circuit
R1 1 2 5k
R2 4 7 15k
R3 2 6 30k
R4 6 0 %f
R6 2 3 5k
R7 3 4 10k
R8 4 5 10k
C1 3 0 1nF ic=%f
C3 4 0 1nF ic=%f
C4 5 0 1nF ic=%f
C2 6 0 360pF ic=%f
Q1 2 5 0 BC547C
Q2 7 6 0 BC547C

.model BC547C NPN(Is=32f Xti=3 Eg=1.11 Vaf=100 Bf=720
+ Xtb=1.5 Br=6.313 Isc=138f Nc=2.053 Rb=1.3k
+ Ikf=500m Ise=230f Ne=1.583
+ Ikr=.187 Rc=.4216 Cjc=5.25p Mjc=.3147 Vjc=.5697 Fc=.5 Cje=11.5p
+ Mje=.3333 Vje=.5 Tr=10n Tf=410.7p Itf=2.204 Xtf=66.65 Vtf=10)

* Other model
.model BC547CV2   NPN(Is=7.049f Xti=3 Eg=1.11 Vaf=24.76 Bf=543.1 Ise=78.17f
+               Ne=1.679 Ikf=94.96m Nk=.5381 Xtb=1.5 Br=1 Isc=27.51f Nc=1.775
+               Ikr=3.321 Rc=.9706 Cjc=5.25p Mjc=.3147 Vjc=.5697 Fc=.5
+               Cje=11.5p Mje=.6715 Vje=.5 Tr=10n Tf=410.7p Itf=1.12 Xtf=26.19
+               Vtf=10)
*       PHILIPS     pid=bc547c  case=TO92
*       91-07-31 dsq

* Simulations
.control
tran 5ns 20ms uic

wrdata %s  v(2) v(3) v(4) v(5) v(6)
exit
.endc
.end
    """ %(R4,v3,v4,v5,v6,outfilename))
    f.close()
    #do simulation
    os.system('ngspice %s' %circuitfilename)
    x = genfromtxt(outfilename)
    os.remove('%s' %outfilename)
    os.remove('%s' %circuitfilename)
    return x


def do_map(V,_):
    #throw away first transient
    t=V[:,0]
    t0=where(t>0.003)[0][0]
    X=V[t0:,7] #1 3 5 7 9
    Y=V[t0:,3]
    Z=V[t0:,5]
    Q=V[t0:,9]

    #rescale to figure same for all runs, rel to 1st run found limits
    X = X - xmin
    X = borderx + X * (xdim-2*borderx-1)/(xmax-xmin)
    Y = Y - ymin
    Y = bordery + Y  * (ydim-2*bordery-1)/(ymax-ymin)

    Z = (Z - zmin) / (zmax-zmin) #z-axis to 0..1
    Q = (Q - zmin) / (zmax-zmin)

    X.clip(0,xdim-1,out=X)
    Y.clip(0,ydim-1,out=Y)
    Z.clip(0,1,out=Z)
    Q.clip(0,1,out=Q)

    #change here how you map the values of the simulation to the image
    #three channels are red, green, blue
    #4th one (transparency) could be added too (but haven't tried)
    m = zeros((ydim,xdim,3))
    x0,y0,z0 = X[0],Y[0],Z[0]
    for x,y,z,q in zip(X.round().astype(int),Y.round().astype(int),Z,Q):
        m[y,x,0] += x
        m[y,x,1] += y
        s = sqrt(pow(x-x0,2)+ pow(y-y0,2) + (z-z0)**2)
        x0,y0,z0 = x,y,z
        m[y,x,2] += s

        #m[y,x,3]=

    return m


if __name__=='__main__':
    random.seed(42)
    R4 = 44.0e3 #for this value of resistor R4, the circuit is in chaotic regime
    numruns = 4
    numthreads= 4 #8
    borderx = 10
    bordery = 10
    xdim = 2000 #1584
    ydim = 2000 #396
    dpi=150
    szx=xdim/dpi
    szy=ydim/dpi
    M = zeros((ydim,xdim,3))

    #do this once to get the extent of the attractor
    V = do_sim(R4,0)
    #throw away first transient
    #throw away first transient
    t=V[:,0]
    t0=where(t>0.003)[0][0]
    X=V[t0:,7]
    Y=V[t0:,3]
    Z=V[t0:,5]
    Q=V[t0:,9]
    xmin,xmax = X.min(),X.max()
    ymin,ymax = Y.min(),Y.max()
    zmin,zmax = Z.min(),Z.max()
    qmin,qmax = Q.min(),Q.max()

    for k in range(numruns):
        print('run {}:'.format(k))
        print('simulating...')

        #simulate
        tasks = [(R4,random.randint(1000000,9999999)) for _ in range(numthreads)]
        pool = multiprocessing.Pool(numthreads)
        taskoutput = [pool.apply_async(do_sim, params) for params in tasks]
        resultslist = [r.get() for r in taskoutput]
        pool.close()
        pool.join()

        #map
        tasks = [(V,0) for V in resultslist]
        pool = multiprocessing.Pool(numthreads)
        taskoutput = [pool.apply_async(do_map, params) for params in tasks]
        resultslist = [r.get() for r in taskoutput]
        pool.close()
        pool.join()

        #add all results
        for m in resultslist:
            M += m

    #nonlinear mapping to enhance the trajectories
    M[:,:,0] /= M[:,:,0].max()
    M[:,:,1] /= M[:,:,1].max()
    M[:,:,2] /= M[:,:,2].max()
    M[:,:,0] = pow(M[:,:,0],1/3)
    M[:,:,1] = pow(M[:,:,1],1/3)
    M[:,:,2] = pow(M[:,:,2],1/3)

    #M[:,:,3] /= M[:,:,3].max()
    #M[:,:,3] = pow(M[:,:,3],1/3)

    print('plotting...')
    fig = plt.figure(frameon=False,figsize=(szx,szy),dpi=dpi)
    ax = fig.add_axes([0, 0, 1, 1])
    ax.axis('off')
    ax.imshow(M)
    plt.savefig('M.png') #,transparant=True
    plt.close('all')
    print('Done!')