In [2]:

    

# Boilerplate
%matplotlib inline
from pylab import rcParams
rcParams['figure.figsize'] = 18, 12
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
import math
import numpy as np
from ClassExpoSin import *


    

In [9]:

    

r1 = 1.0
r2 = 1.1
psi = 9*math.pi


    

In [10]:

    

def k2_simple_range(r1, r2):
    """The possible-k2 domain falls within this range."""
    return -math.pow(4/math.log(r1/r2)**2, 0.25), math.pow(4/math.log(r1/r2)**2, 0.25)
def check_k2_possible(k2, r1, r2, psi):
    """Check that delta in equation 11 is positive."""
    return 0 <= 2 * (1 - math.cos(k2 * psi)) / k2**4 - math.log(r1 / r2)**2

def tany1r(k2, r1, r2, psi):
    logr1r2 = math.log(r1 / r2)
    cosk2O = math.cos(k2 * psi)
    delta = 2*(1-cosk2O)/k2**4 - logr1r2**2
    if delta < 0: # no feasible trajectories
        return [None, None, None]
    tany1mid = k2/2 * (-logr1r2 / math.tan(k2*psi/2))
    tany1min = k2/2 * (-logr1r2 / math.tan(k2*psi/2) - math.sqrt(delta))
    tany1max = k2/2 * (-logr1r2 / math.tan(k2*psi/2) + math.sqrt(delta))
    return (tany1min, tany1max, tany1mid)

k2s = linspace(0.001, k2_simple_range(r1,r2)[1], 10000)
xlabel('k2')
ylabel('tany1 min, max, center')
title('tany1 range for given k2')
# Plot where no valid tany1 exists (thus k2 impossible)
plot(k2s, [None if check_k2_possible(k2,r1,r2,psi) else 0 for k2 in k2s], 'r')
# Plot tany1 range
plot(k2s, [tany1r(k2,r1,r2,psi)[0] for k2 in k2s], 'k')
plot(k2s, [tany1r(k2,r1,r2,psi)[1] for k2 in k2s], 'k')
plot(k2s, [tany1r(k2,r1,r2,psi)[2] for k2 in k2s], 'k--')


    

    
Out[10]:





[<matplotlib.lines.Line2D at 0x7f051df152d0>]






    

In [12]:

    

k2s = linspace(0.001,k2_simple_range(r1,r2)[1],200)
tany1s = linspace(-15.0,15.0,200)
k2s, tany1s = meshgrid(k2s, tany1s)

@np.vectorize
def dt(k2, tany1):
    try:
        clss = ClassExpoSin(k2, r1, r2, psi)
        mintany1, maxtany1 = clss.tany1Range()
        extinct = 1.0
        if tany1 > maxtany1 * extinct or tany1 < mintany1 * extinct:
            return np.nan
        exps = clss.createExpoSin(tany1)
        return math.log(exps.dT(psi, 1.0))
    except:
        return np.nan

dts = dt(k2s,tany1s)

mdts = np.ma.array(dts, mask=np.isnan(dts))
cmap = cm.hsv
cmap.set_bad('w',1.0)

fig = plt.figure()
ax = fig.gca()
ax.set_title('TOF matrix for k2, tany1')
surf = ax.pcolor(k2s, tany1s, mdts, cmap=cmap, vmin=np.nanmin(dts), vmax=np.nanmax(dts))
plt.show()


    

In [ ]: