data!
system:sage

{{{id=1|
#auto
def dist(v, b, left=-1.0r, right=1.0r):
    """
    We divide the interval between left (default: 0) and 
    right (default: pi) up into b bins.
   
    For each number in v (which must left and right), 
    we find which bin it lies in and add this to a counter.
    This function then returns the bins and the number of
    elements of v that lie in each one. 

    ALGORITHM: To find the index of the bin that a given 
    number x lies in, we multiply x by b/length and take the 
    floor. 
    """
    length = right - left
    normalize = float(b/length)
    vals = {}
    d = dict([(i,0) for i in range(b)])
    for x in v:
        n = int(normalize*(float(x)-left))
        d[n] += 1
    return d, len(v)

def frequency_histogram(d, b, num=5000, left=-1.0r, right=1.0r):
    s = Graphics()
    left = float(left); right = float(right)
    length = right - left
    w = length/b
    k = 0
    for i, n in d.iteritems():
        k += n
        # ith bin has n objects in it. 
        s += polygon([(w*i+left,0), (w*(i+1)+left,0), \
                     (w*(i+1)+left, n/(num*w)), (w*i+left, n/(num*w))],\ 
                     rgbcolor=(0,0,0.5))
    return s

def semicircle0():
    alpha = 2.0r / float(pi)
    def f(x):
        return alpha * math.sqrt(1-x^2)
    return plot(f, -1,1,rgbcolor=(0.5,0,0))
semicircle = semicircle0()
}}}

{{{id=2|
#auto
from math import asin, log, sqrt

def redline(xmin,xmax):
    return line([(xmin,1/2),(xmax,1/2)], rgbcolor=(1,0,0))

def Xab(a,b):
    bb = (asin(b)/2r + b*sqrt(1r-b^2r)/2r)
    aa = (asin(a)/2r + a*sqrt(1r-a^2r)/2r)
    def X(T):
        return (asin(T)/2r + T*sqrt(1r-T^2r)/2r - aa)/(bb - aa)
    return X

import bisect

class SatoTate:
    def __init__(self, E, n):
        self._E = E    
        self._n = n
        self.init_aplist(n)

    def init_aplist(self, n):
        t = cputime()
        v = self._E.aplist(n, python_ints=True)
        print 'computed aplist ', cputime(t)
        P = prime_range(n)
        self._aplist = v
        two = float(2)
        t = cputime()
        self._normalized_aplist = [float(v[i])/(two*math.sqrt(P[i])) for i in range(len(v))]
        print 'time to normalize ', cputime(t)

    def __repr__(self): 
        return "Sato-Tate data for %s using primes up to %s"%(self._E, self._n)
 
    def normalized_aplist(self, n):  # returns a copy
        k = prime_pi(n)
        v = self._normalized_aplist
        if k > len(v):
             raise ValueError, "call init_aplist"
        return v[:prime_pi(n)]

    def sorted_aplist(self, n):
        v = self.normalized_aplist(n)
        v.sort()
        return v

    def YCab(self, Cmax, a=-1, b=1):
        v = self.sorted_aplist(Cmax)
    
        denom = bisect.bisect_right(v, float(b)) - bisect.bisect_left(v, float(a))
        try:
           normalize = float(1)/denom
        except:
           def Y(T):
               return 1.0r
           return Y
        start_pos = bisect.bisect_left(v, float(a))

        def Y(T):
            # find position that T would go in if it were inserted
            # in the sorted list v.
            n = bisect.bisect_right(v, float(T)) - start_pos
            return n * normalize
        return Y


    def xyplot(self, C, a=-1, b=1):
        """
        Return the quantile-quantile plot for given a,b, up to C.
        """
        Y = self.YCab(C,a=a,b=b)
        X = Xab(a=a,b=b)
        pX = plot(X, a, b, rgbcolor=(1,0,0))
        pY = plot(Y, a, b, rgbcolor=(0,0,1))
        return pX + pY

    def qqplot(self, C, a=-1, b=1):
        """
        Return the quantile-quantile plot for given a,b, up to C.
        """
        Y = self.YCab(C,a=a,b=b)
        X = Xab(a=a,b=b)
        pl = parametric_plot((X, Y), a,b)
        ll = line([(0,0), (1.1,1.1)], rgbcolor=(1,0,0))
        return pl+ll
 
    def Delta(self, C, a, b, max_points=300, L_norm=2):
        """
        Delta_{a}^{b} function:

        INPUT: C - cutoff
             a,b - evaluate over the interval (a,b)
             max_points - number of points used in numerical integral         
             L_norm --the integer n=2 or n=oo. 
                      Compute the L_n norm.   For n finite this
                      is the integral of the difference to the power n.
                      For n = +oo, this is the L_oo norm, which is the max
                      of the absolute value of the difference (where the max
                      is evaluated at max_points equidistributed points).
        """
        key = (C,a,b,max_points, L_norm)
        try:
           return self._delta[key]
        except AttributeError:
           self._delta = {}
        except KeyError:
           pass
        X = Xab(a,b)
        Y = self.YCab(C,a,b)     # This takes all the time.

        if L_norm == oo:
            val = max([abs(X(T)-Y(T)) for T in srange(a,b,float(b-a)/max_points)])
            err = 0
        else:
            n = int(L_norm)   # usually n = 2.
            def h(T):
                return (X(T) - Y(T))^n
            val, err = integral_numerical(h, a, b, max_points=max_points, 
              algorithm='qag', rule=1,eps_abs=1e-10, eps_rel=1e-10)
            val = float(val^(1.0/n)) # compute L_n norm.
            err = float(err^(1.0/n))
    
        self._delta[key] = (val, err)
        return val, err

    def plot_Delta(self, Cmax=None, **kwds):
        if Cmax is None: Cmax = self._n
        v2 = self.compute_Delta(Cmax, L_norm=2,  **kwds)
        voo = self.compute_Delta(Cmax, L_norm=oo, **kwds)
        return line(v2, rgbcolor=(0,0,1)) + line(voo, rgbcolor=(0.8,0,0))

    def compute_Delta(self, Cmax, plot_points=30, a=-1, b=1, 
                            max_points=300, L_norm=2, verbose=False):
        a,b = (float(a), float(b))
        def f(C):
           val, err = self.Delta(C, a, b, max_points=max_points, L_norm=L_norm)
           return val
        return [(x,f(x)) for x in range(100, Cmax, int(Cmax/plot_points))]

    def theta(self, C, a=-1, b=1, max_points=300, L_norm=2):
        val, err = self.Delta(C, a, b, max_points=max_points, L_norm=L_norm)
        return -log(val)/log(C), val, err

    def theta_interval(self, C, a=-1, b=1, max_points=300, L_norm=2):
        val, err = self.Delta(C, a, b, max_points=max_points, L_norm=L_norm)
        return -log(val-abs(err))/log(C), -log(val+abs(err))/log(C)

    def compute_theta(self, Cmax, plot_points=30, a=-1, b=1, 
                 max_points=300, L_norm=2, verbose=False):
        a,b = (float(a), float(b))
        def f(C):
           z = self.theta(C, a, b, max_points=max_points, L_norm=L_norm)
           if verbose: print C, z
           return z[0]
        return [(x,f(x)) for x in range(100, Cmax, int(Cmax/plot_points))]

    def compute_theta_interval(self, Cmax, plot_points=30, a=-1, b=1, 
                max_points=300, L_norm=2, verbose=False):
        a,b = (float(a), float(b))
        vmin = []; vmax = []
        for C in range(100, Cmax, int(Cmax/plot_points)):
            zmin,zmax = self.theta_interval(C, a, b, max_points=max_points, L_norm=L_norm)
            vmin.append((C, zmin))
            vmax.append((C, zmax))
            if verbose: print C, zmin, zmax
        return vmin, vmax

    def plot_theta_interval(self, Cmax=None, clr=(0,0,1), *args, **kwds):
        if Cmax is None: Cmax = self._n
        vmin, vmax = self.compute_theta_interval(Cmax, *args, **kwds)
        v = self.compute_theta(Cmax, *args, **kwds)
        grey = (0.7,0.7,0.7)
        return line(vmin,rgbcolor=grey)+line(vmax,rgbcolor=grey) \
                 + line(v,rgbcolor=clr) + redline(0, Cmax) # + point(v,rgbcolor=(0,0,0), pointsize=3) 

    def histogram(self, Cmax, num_bins):
        v = self.normalized_aplist(Cmax)
        d, total_number_of_points = dist(v, num_bins)
        return frequency_histogram(d, num_bins, total_number_of_points) + semicircle

    def x_times_Delta(self, x):
        return x*self.Delta(x, -1,1, max_points=500)[0]

def showtheta(P):
    P.show(xmin=0, ymin=0, ymax=0.6, figsize=[9,4])
}}}

{{{id=44|
# make curve
time e37a = SatoTate(EllipticCurve('37a'), 10^6)
///
computed aplist  5.120321
time to normalize  0.080005
Time: CPU 5.85 s, Wall: 7.80 s
}}}

{{{id=47|
show(e37a.plot_Delta(10^3, plot_points=400, max_points=100), ymax=0.1, ymin=0, figsize=[10,3])
}}}

{{{id=45|
show(e37a.plot_Delta(10^4, plot_points=200, max_points=100), ymax=0.1, ymin=0, figsize=[10,3])
}}}

{{{id=46|
show(e37a.plot_Delta(10^6, plot_points=200, max_points=100), ymax=0.1, ymin=0, figsize=[10,3])
}}}

{{{id=41|
# make curve
time S0 = SatoTate(EllipticCurve('5077a'), 10^6)
///
computed aplist  5.356334
time to normalize  0.072004
Time: CPU 6.06 s, Wall: 8.39 s
}}}

{{{id=42|
show(S0.plot_Delta(10^6, plot_points=100, max_points=100), ymax=0.1, ymin=0, figsize=[10,3])
}}}

{{{id=43|
show(S0.plot_Delta(10^5, plot_points=100, max_points=100, L_norm=oo), ymax=0.1, ymin=0, figsize=[10,3])
}}}

{{{id=49|
# make curve
time S0 = SatoTate(EllipticCurve('37a'), 10^6)
///
computed aplist  5.452341
time to normalize  0.088005
Time: CPU 6.23 s, Wall: 7.71 s
}}}

{{{id=50|
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
Time: CPU 3.36 s, Wall: 3.82 s
CPU time: 8.47 s,  Wall time: 9.97 s
}}}

{{{id=48|
#fit = plot(lambda x: 1/2 - 1/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
}}}

{{{id=51|
showtheta(P0+P1)
}}}

{{{id=54|
# make curve
time S0 = SatoTate(EllipticCurve('11a'), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  6.004376
time to normalize  0.084006
Time: CPU 6.75 s, Wall: 8.61 s
Time: CPU 10.05 s, Wall: 11.99 s
Time: CPU 6.80 s, Wall: 7.73 s
}}}

{{{id=53|
showtheta(P0+P1)
}}}

{{{id=58|
# make curve
time S0 = SatoTate(EllipticCurve('389a'), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.812363
time to normalize  0.088006
Time: CPU 6.58 s, Wall: 8.91 s
Time: CPU 10.13 s, Wall: 14.74 s
Time: CPU 7.09 s, Wall: 9.34 s
}}}

{{{id=57|
showtheta(P0+P1)
}}}

{{{id=56|
# make curve
time S0 = SatoTate(EllipticCurve('5077a'), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.724357
time to normalize  0.084005
Time: CPU 6.44 s, Wall: 8.98 s
Time: CPU 10.18 s, Wall: 14.21 s
Time: CPU 6.98 s, Wall: 9.24 s
}}}

{{{id=60|
showtheta(P0+P1)
}}}

{{{id=66|
# make curve
time S0 = SatoTate(EllipticCurve('5077a'), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.372336
time to normalize  0.076004
Time: CPU 6.09 s, Wall: 7.20 s
Time: CPU 9.86 s, Wall: 12.41 s
Time: CPU 7.06 s, Wall: 8.90 s
}}}

{{{id=86|
fit = plot(lambda x: 1/2 - 1/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
}}}

{{{id=85|
showtheta(P0 + P1 + fit)
}}}

{{{id=84|
%time
time S0 = SatoTate(EllipticCurve([19,234]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.416338
time to normalize  0.076005
Time: CPU 6.09 s, Wall: 7.23 s
Time: CPU 9.72 s, Wall: 11.44 s
Time: CPU 6.84 s, Wall: 8.45 s
CPU time: 22.64 s,  Wall time: 27.12 s
}}}

{{{id=88|
fit = plot(lambda x: 1/2 - 1/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
}}}

{{{id=87|
showtheta(P0 + P1 + fit)
}}}

{{{id=91|
time S0 = SatoTate(EllipticCurve([1,-1,0,-79,289]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.704355
time to normalize  0.116007
Time: CPU 9.35 s, Wall: 9.74 s
Time: CPU 6.41 s, Wall: 6.47 s
CPU time: 22.27 s,  Wall time: 23.47 s
}}}

{{{id=100|
log=math.log
r = 4.0
fit_watkins = plot(lambda x: 1/2 - (2/r)*log(log(x))/log(x), 10, 10^6, rgbcolor=(0,0.5,0.5))
}}}

{{{id=92|
fit = plot(lambda x: 1/2 - (4/3)/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
}}}

{{{id=90|
showtheta(P0 + P1 + fit + fit_watkins)
}}}

{{{id=104|

}}}

{{{id=103|
time S3 = SatoTate(EllipticCurve('5077a'), 10^6)
time P0_3 = S3.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1_3 = S3.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
///
computed aplist  5.212326
time to normalize  0.088006
Time: CPU 9.36 s, Wall: 9.63 s
Time: CPU 6.62 s, Wall: 6.74 s
CPU time: 21.85 s,  Wall time: 22.40 s
}}}

{{{id=102|
r = 4.0
fit_watkins = plot(lambda x: 1/2 - (2/r)*log(log(x))/log(x), 10, 10^6, rgbcolor=(0.2,0.1,0.5))
}}}

{{{id=89|
showtheta(P0_3 + P1_3 + fit_watkins)
}}}

{{{id=83|

}}}

{{{id=65|
time S0 = SatoTate(EllipticCurve([1,-1,0,-79,289]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.708357
time to normalize  0.084005
Time: CPU 9.96 s, Wall: 13.32 s
Time: CPU 6.73 s, Wall: 9.17 s
Time: CPU 0.21 s, Wall: 0.24 s
CPU time: 23.33 s,  Wall time: 31.74 s
}}}

{{{id=59|

}}}

{{{id=69|

}}}

{{{id=68|
time S0 = SatoTate(EllipticCurve([0, 0, 1, -79, 342]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.496343
time to normalize  0.080005
Time: CPU 10.03 s, Wall: 12.37 s
Time: CPU 6.89 s, Wall: 9.98 s
Time: CPU 0.22 s, Wall: 0.25 s
CPU time: 23.34 s,  Wall time: 31.16 s
}}}

{{{id=71|
time S0 = SatoTate(EllipticCurve([1, 1, 0, -2582, 48720]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.60835
time to normalize  0.084005
Time: CPU 9.56 s, Wall: 11.41 s
Time: CPU 6.55 s, Wall: 7.51 s
Time: CPU 0.22 s, Wall: 0.22 s
CPU time: 22.69 s,  Wall time: 26.81 s
}}}

{{{id=82|
time S0 = SatoTate(EllipticCurve([0, 0, 0, -10012, 346900]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.44034
time to normalize  0.068004
Time: CPU 9.69 s, Wall: 11.92 s
Time: CPU 6.94 s, Wall: 8.90 s
Time: CPU 0.19 s, Wall: 0.20 s
CPU time: 22.93 s,  Wall time: 28.44 s
}}}

{{{id=81|
time S0 = SatoTate(EllipticCurve([0, 0, 1, -23737, 960366]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.816364
time to normalize  0.084005
Time: CPU 9.72 s, Wall: 12.06 s
Time: CPU 6.70 s, Wall: 8.50 s
Time: CPU 0.21 s, Wall: 0.23 s
CPU time: 23.30 s,  Wall time: 29.06 s
}}}

{{{id=99|
fit = plot(lambda x: 1/2 - (19/9)/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
showtheta(P0 + P1 + fit)
}}}

{{{id=80|
time S0 = SatoTate(EllipticCurve([1,-1,1,-20067762415575526585033208209338542750930230312178956502,34481611795030556467032985690390720374855944359319180361266008296291939448732243429 ]), 10^6)
time P0 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  6.008374
time to normalize  0.084004
Time: CPU 9.82 s, Wall: 13.54 s
Time: CPU 6.54 s, Wall: 7.59 s
Time: CPU 0.20 s, Wall: 0.43 s
CPU time: 23.27 s,  Wall time: 29.89 s
}}}

{{{id=79|
fit = plot(lambda x: 1/2 - (28/9)/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
showtheta(P0 + P1 + fit)
}}}

{{{id=78|

}}}

{{{id=98|
time S8 = SatoTate(EllipticCurve([0, 0, 1, -23737, 960366]), 10^6)
time P0 = S8.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=2, clr=(0,0,1))
time P1 = S8.plot_theta_interval(10^6, plot_points=300, max_points=80, L_norm=oo, clr=(0,0.5,0))
time showtheta(P0+P1)
///
computed aplist  5.632352
time to normalize  0.080005
Time: CPU 16.38 s, Wall: 45.97 s
Time: CPU 11.85 s, Wall: 18.06 s
Time: CPU 1.67 s, Wall: 3.15 s
CPU time: 36.25 s,  Wall time: 74.68 s
}}}

{{{id=97|
fit = plot(lambda x: 1/2 - 2/math.log(x), 10, 10^6, rgbcolor=(0,0,0))
}}}

{{{id=96|
showtheta(P0 + P1 + fit)
}}}

{{{id=95|

}}}

{{{id=94|

}}}

{{{id=93|

}}}

{{{id=77|

}}}

{{{id=76|

}}}

{{{id=75|

}}}

{{{id=74|

}}}

{{{id=73|

}}}

{{{id=72|

}}}

{{{id=67|

}}}

{{{id=55|

}}}

{{{id=52|

}}}

{{{id=35|
# make curve
time S0 = SatoTate(EllipticCurve('5077a'), 5*10^6)
///
computed aplist  50.23914
time to normalize  1.796112
Time: CPU 59.98 s, Wall: 72.94 s
}}}

{{{id=38|
time P0 = S0.plot_theta_interval(5*10^6, plot_points=100, max_points=100, L_norm=2, clr=(0,0,1))
time P1 = S0.plot_theta_interval(5*10^6, plot_points=100, max_points=100, L_norm=oo, clr=(0,0.5,0))
///
Time: CPU 12.93 s, Wall: 14.71 s
CPU time: 29.03 s,  Wall time: 34.74 s
}}}

{{{id=40|
fit = plot(lambda x: 1/2 - 1/math.log(x), 10, 5*10^6, rgbcolor=(0,0,0))
}}}

{{{id=36|
showtheta(P0+P1+fit)
}}}

{{{id=63|

}}}

{{{id=62|

}}}

{{{id=61|

}}}

{{{id=39|

}}}

{{{id=3|
S0 = SatoTate(EllipticCurve('11a'))
S1 = SatoTate(EllipticCurve('37a'))
S2 = SatoTate(EllipticCurve('389a'))
S3 = SatoTate(EllipticCurve('5077a'))
S4 = SatoTate(EllipticCurve([1,-1,0,-79,289]))
S5 = SatoTate(EllipticCurve([0, 0, 1, -79, 342]))
S6 = SatoTate(EllipticCurve([1, 1, 0, -2582, 48720]))
S7 = SatoTate(EllipticCurve([0, 0, 0, -10012, 346900]))
S8 = SatoTate(EllipticCurve([0, 0, 1, -23737, 960366]))
S28 = SatoTate(EllipticCurve([1,-1,1,\
  -20067762415575526585033208209338542750930230312178956502,\
  34481611795030556467032985690390720374855944359319180361266008296291939448732243429 ]))
}}}

{{{id=8|
%time
P0 = S0.plot_theta_interval(10^5, plot_points=10, max_points=10, L_norm=2)
P0.show(dpi=100)
///
CPU time: 0.18 s,  Wall time: 0.20 s
}}}

{{{id=0|
%time
P1 = S1.plot_theta_interval(10^6,plot_points=15,max_points=50)
P1.show(dpi=200)
}}}

{{{id=4|
%time
P2 = S2.plot_theta_interval(10^6,plot_points=15,max_points=50)
P2.show(dpi=200)
}}}

{{{id=5|
%time
P3 = S3.plot_theta_interval(10^6,plot_points=15,max_points=50)
P3.show(dpi=200)
}}}

{{{id=34|
%time
P8 = S8.plot_theta_interval(10^6,plot_points=15,max_points=50)
P8.show(dpi=200)
}}}

{{{id=6|
%time
P0oo = S0.plot_theta_interval(10^6,plot_points=15,max_points=50, L_norm=oo)
P0oo.show(dpi=200)
}}}

{{{id=7|
%time
P1oo = S1.plot_theta_interval(10^6,plot_points=15,max_points=50, L_norm=oo)
P1oo.show(dpi=200)
}}}

{{{id=9|
%time
P2oo = S2.plot_theta_interval(10^6,plot_points=15,max_points=50, L_norm=oo)
P2oo.show(dpi=200)
}}}

{{{id=10|
%time
P3oo = S3.plot_theta_interval(10^6,plot_points=15,max_points=50, L_norm=oo)
P3oo.show(dpi=200)
}}}

{{{id=11|
%time
P8oo = S8.plot_theta_interval(10^6,plot_points=15,max_points=50, L_norm=oo)
P8oo.show(dpi=200)
}}}