sstability.py

from matplotlib import pyplot as plt
from scipy import arange, meshgrid, zeros, cos, sin, log10, linspace, logspace, around, exp, absolute
from scipy.optimize import fmin

from idmlin import IDMLin


def maxStringResponse(idm, n, tau = 0):
    gammainv = humanstring(idm, n, tau)
    wmax = fmin(gammainv, 10)
    # print "wmax =",  wmax
    return -gammainv(wmax)

def human(idm, tau=0):
    f, g, h = idm.f, idm.g, idm.h
    k = g+h
    def gamma(w):
        s = w * 1j
        num = f + g * s
        den = f + k * s + exp(tau * s) * s * s
        return num/den
    return gamma

def humanstring(idm, n, tau = 0):
    gamma = human(idm, tau)
    def dbgammainv(w):
        db = 20 * log10(absolute(gamma(w)))
        return -n*db
    return dbgammainv

def followedrobot(idm):#, ff, fb, gf, gb, hr, taur=0, tau=0):
    tau = 0
    taur = 0
    ff = idm.f/2
    fb = idm.f/2
    gf = idm.g/2
    gb = idm.g/2
    hr = idm.h

    gamma = human(idm, tau)
    def dbgammainv(w):
        s = w * 1j
        num = ff + gf * s
        den = ff + fb + (gf + gb + hr) * s + exp(taur * s) * s * s - (fb + gb * s) * gamma(w)
        dbnum = 20 * log10(absolute(num))
        dbden = 20 * log10(absolute(den))
        return (dbden-dbnum)
    return dbgammainv

def teststability(idm, N):
    mr = N * idm.maxresponse()
    return max(mr, 0)

a = 0.45
t = 1.5
b = 2.0
s = 2.0
v = 33.3

idm = IDMLin(t=t, a=a, b=b, s0=s)

nrange = logspace(1, 2.5, 20)
rerange = logspace(-2, -0.2, 20)

mgrid = zeros((len(nrange), len(rerange)))
mgrid2 = zeros((len(nrange), len(rerange)))

for i, n in enumerate(nrange):
    for j, re in enumerate(rerange):
        idm.go(v, re)
        mgrid[i,j] = teststability(idm, n)
        mgrid2[i,j] = maxStringResponse(idm, n)

fig, axarr = plt.subplots(1,3, sharex=True, sharey=True)

im = axarr[2].imshow(mgrid, cmap=plt.cm.Reds, 
                interpolation='none', 
                aspect='auto',
		vmin=0, vmax=1,
                )
im = axarr[1].imshow(mgrid2, cmap=plt.cm.Reds, 
                interpolation='none', 
                aspect='auto',
		vmin=0, vmax=1,
                )
im = axarr[0].imshow((mgrid2-mgrid)*100, cmap=plt.cm.Reds, 
                interpolation='none', 
                aspect='auto',
		vmin=0, vmax=1,
                )

axarr[0].set_yticks(range(len(nrange)))
axarr[1].set_yticks(range(len(nrange)))
axarr[1].set_xticks(range(len(rerange)))

axarr[0].set_yticklabels(nrange)
axarr[1].set_yticklabels(nrange)
axarr[1].set_xticklabels(around(rerange, 2))

axarr[0].set_ylabel("N")
axarr[1].set_ylabel("N")
axarr[1].set_xlabel("R_e")

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im, cax=cbar_ax)

plt.show()