notebook.community

Edit and run



In [2]:

    
import healpix_util as hu
import astropy as ap
import math, sys, time
import numpy as np
from astropy.io import fits
from astropy.table import Table
import astropy.io.ascii as ascii
from astropy.constants import c
import matplotlib.pyplot as plt
import math
import scipy.special as sp
from scipy import integrate
from astropy.stats import bayesian_blocks
from astropy.stats import histogram
from astropy.visualization import hist
%matplotlib inline
#%precision 4
#plt.style.use('ggplot')



In [16]:

    
# Open file
f = open('./output/DD_1ff_simpler.dat', 'r')
dddata = open("./output/DD_1_R_alpha.dat",'w')
dddata.write("deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\t thetadd \t DD \n")

# Read and ignore header lines
header1 = f.readline()
#print header1

# Loop over lines and extract variables of interest
for line in f:
    line = line.strip()
    columns = line.split()
    deltaZ = float(columns[6])
    deltatheta = float(columns[7])
    Zdeltatheta = float(columns[8])
    R = float(columns[9])
    alpha = float(columns[10])
    thetadd = float(columns[11])
    DD = float(columns[12])
    if Zdeltatheta<0.02:
        dddata.write("%f\t" %deltaZ)
        dddata.write("%f\t" %deltatheta)
        dddata.write("%f\t" %Zdeltatheta)
        dddata.write("%f\t" %R)
        dddata.write("%f\t" %alpha)
        dddata.write("%f \t %f \n"%(thetadd,DD))
    else:
        continue
f.close()
dddata.close()









    



Z1	 RA1	 DEC1	 Z2	 RA2	 DEC2	 deltaZ	 deltatheta	 Zdeltatheta	 R	 alpha	 thetadd 	 DD



In [3]:

    
hmqdat1=ascii.read("./output/DD_1_R_alpha.dat")
hmqdat1
#Z=sum(hmqdat1['Z'])/len(hmqdat1['Z'])
#print Z









    Out[3]:




<Table length=1305226>

deltaZ deltatheta Zdeltatheta R alpha thetadd DD
float64 float64 float64 float64 float64 float64 float64
0.013 0.051491 0.012259 0.017869 0.756076 0.052368 1.0
0.013 0.040922 0.009743 0.016246 0.64315 0.041069 1.0
0.008 0.062346 0.014843 0.016862 1.07647 1.878369 1.0
0.009 0.068156 0.016227 0.018556 1.064399 1.893666 1.0
0.01 0.067917 0.01617 0.019012 1.016931 1.907749 1.0
0.001 0.051615 0.012289 0.012329 1.489599 0.992046 1.0
0.007 0.083981 0.019994 0.021184 1.234034 1.219094 1.0
0.017 0.082431 0.019625 0.025964 0.856957 0.604861 1.0
0.018 0.070245 0.016724 0.02457 0.74867 1.063174 1.0
0.001 0.040053 0.009536 0.009588 1.466311 2.505325 1.0
... ... ... ... ... ... ...
0.0 0.079949 0.019034 0.019034 1.570796 2.53098 1.0
0.0 0.053744 0.012796 0.012796 1.570796 0.988719 2.0
0.0 0.075223 0.017909 0.017909 1.570796 2.53098 2.0
0.0 0.046556 0.011084 0.011084 1.570796 0.538925 1.0
0.0 0.068345 0.016272 0.016272 1.570796 1.952672 1.0
0.0 0.041424 0.009862 0.009862 1.570796 2.466464 1.0
0.0 0.033226 0.007911 0.007911 1.570796 0.988719 1.0
0.0 0.061433 0.014626 0.014626 1.570796 2.551203 1.0
0.0 0.067922 0.016171 0.016171 1.570796 2.551203 1.0
0.0 0.061536 0.014651 0.014651 1.570796 0.988719 1.0



In [4]:

    
hmqdat1.sort('R')
hmqdat1









    Out[4]:




<Table length=1305226>

deltaZ deltatheta Zdeltatheta R alpha thetadd DD
float64 float64 float64 float64 float64 float64 float64
0.0 5e-06 1e-06 1e-06 1.570796 0.538925 4.0
0.0 7e-06 2e-06 2e-06 1.570796 0.538925 4.0
0.0 1.2e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.1e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.3e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.7e-05 4e-06 4e-06 1.570796 0.538925 9.0
0.0 2.1e-05 5e-06 5e-06 1.570796 0.538925 4.0
0.0 2.2e-05 5e-06 5e-06 1.570796 0.538925 4.0
0.0 2e-05 5e-06 5e-06 1.570796 0.538925 169.0
0.0 2.2e-05 5e-06 5e-06 1.570796 0.538925 4.0
... ... ... ... ... ... ...
0.02 0.083981 0.019995 0.02828 0.785262 0.891568 1.0
0.02 0.083978 0.019994 0.02828 0.785239 2.044048 1.0
0.02 0.083987 0.019996 0.028281 0.785292 0.501326 1.0
0.02 0.083984 0.019995 0.028281 0.785277 0.792125 1.0
0.02 0.083988 0.019996 0.028282 0.785301 2.064619 2.0
0.02 0.083989 0.019996 0.028282 0.785309 0.111224 1.0
0.02 0.083999 0.019999 0.028283 0.785368 0.599163 1.0
0.02 0.083998 0.019998 0.028283 0.785359 0.334718 1.0
0.02 0.084 0.019999 0.028284 0.785372 0.909261 1.0
0.02 0.084 0.019999 0.028284 0.785372 0.704685 2.0



In [22]:

    
#theta=np.arange(0,np.pi,0.001)
#correldat = np.polynomial.legendre.legval(np.cos(theta),(2*ell+1)*np.absolute(pixcl))/(4*math.pi)
plt.figure()
plt.plot(hmqdat1['R'],hmqdat1['DD'])
plt.xlabel('R')
plt.ylabel('DD')
plt.savefig('DDvsR.pdf')
plt.savefig('DDvsR.jpeg')
plt.show()



In [ ]:

    
from scipy.interpolate import spline

Rnew = np.linspace(hmqdat1['R'].min(),hmqdat1['R'].max(),300) #300 represents number of points to make between T.min and T.max

DD_smooth = spline(hmqdat1['R'],hmqdat1['DD'],Rnew)

plt.plot(Rnew,DD_smooth)
plt.show()



In [ ]:



In [14]:

    
hmqdat1[0]['thetadd']
#range(0,len(hmqdat1['thetadd']))









    Out[14]:





0.53892499999999999



In [10]:

    
help(range)









    



Help on built-in function range in module __builtin__:

range(...)
    range(stop) -> list of integers
    range(start, stop[, step]) -> list of integers
    
    Return a list containing an arithmetic progression of integers.
    range(i, j) returns [i, i+1, i+2, ..., j-1]; start (!) defaults to 0.
    When step is given, it specifies the increment (or decrement).
    For example, range(4) returns [0, 1, 2, 3].  The end point is omitted!
    These are exactly the valid indices for a list of 4 elements.



In [12]:

    
len(hmqdat1['thetadd'])









    Out[12]:





1305226



In [16]:

    
hmqdat1.sort('thetadd')



In [17]:

    
hmqdat1['thetadd']









    Out[17]:




<Column name='thetadd' dtype='float64' length=1305226>

0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
0.01812
...
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586
3.060586



In [21]:

    
ddtdata = open("./output/theta_DD.dat",'w')
ddtdata.write("thetadd \t DD \n")
DDsum=hmqdat1[0]['DD']
for i in range(1,len(hmqdat1['thetadd'])-1):
    #print DDsum
    if hmqdat1[i-1]['thetadd']==hmqdat1[i]['thetadd']:
        DDsum=DDsum+hmqdat1[i]['DD']
    else:
        ddtdata.write("%f \t %f \n"%(hmqdat1[i-1]['thetadd'],DDsum))
        DDsum=hmqdat1[i-1]['DD']



In [22]:

    
thedd=ascii.read("./output/theta_DD.dat")



In [23]:

    
plt.figure()
plt.plot(thedd['thetadd'],thedd['DD'])
plt.xlabel('thetadd')
plt.ylabel('DD')
plt.savefig('DDvstheta.pdf')
plt.savefig('DDvstheta.jpeg')
plt.show()



In [24]:

    
hist(thedd['thetadd'])









    Out[24]:





(array([ 261.,  579.,  939.,  980.,  870.,  191.,  194.,  362.,  323.,  123.]),
 array([ 0.01812 ,  0.268416,  0.518712,  0.769008,  1.019304,  1.2696  ,
         1.519896,  1.770192,  2.020488,  2.270784,  2.52108 ]),
 <a list of 10 Patch objects>)



In [26]:

    
max(hmqdat1['thetadd'])









    Out[26]:





3.0605859999999998



In [27]:

    
len(hmqdat1['thetadd'])









    Out[27]:





1305226



In [30]:

    
hmqdat1[1305224]['thetadd']









    Out[30]:





3.0605859999999998



In [40]:

    
bl=histogram(hmqdat1['thetadd'],'blocks')



In [37]:

    
help(histogram)









    



Help on function histogram in module astropy.stats.histogram:

histogram(a, bins=10, range=None, weights=None, **kwargs)
    Enhanced histogram function, providing adaptive binnings
    
    This is a histogram function that enables the use of more sophisticated
    algorithms for determining bins.  Aside from the ``bins`` argument allowing
    a string specified how bins are computed, the parameters are the same
    as ``numpy.histogram()``.
    
    Parameters
    ----------
    a : array_like
        array of data to be histogrammed
    
    bins : int or list or str (optional)
        If bins is a string, then it must be one of:
    
        - 'blocks' : use bayesian blocks for dynamic bin widths
    
        - 'knuth' : use Knuth's rule to determine bins
    
        - 'scott' : use Scott's rule to determine bins
    
        - 'freedman' : use the Freedman-Diaconis rule to determine bins
    
    range : tuple or None (optional)
        the minimum and maximum range for the histogram.  If not specified,
        it will be (x.min(), x.max())
    
    weights : array_like, optional
        Not Implemented
    
    other keyword arguments are described in numpy.histogram().
    
    Returns
    -------
    hist : array
        The values of the histogram. See ``normed`` and ``weights`` for a
        description of the possible semantics.
    bin_edges : array of dtype float
        Return the bin edges ``(length(hist)+1)``.
    
    See Also
    --------
    numpy.histogram



In [39]:

    
hist(hmqdat1['thetadd'],'blocks')









    Out[39]:





(array([ 306.,   60.,  215., ...,  409.,  447.,  680.]),
 array([ 0.01812  ,  0.0211935,  0.025516 , ...,  2.980921 ,  2.988078 ,
         3.060586 ]),
 <a list of 3966 Patch objects>)



In [43]:

    
plt.plot(bl[0],bl[1])









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-43-8316ea71b986> in <module>()
----> 1 plt.plot(bl[0],bl[1])

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/pyplot.pyc in plot(*args, **kwargs)
   3316                       mplDeprecation)
   3317     try:
-> 3318         ret = ax.plot(*args, **kwargs)
   3319     finally:
   3320         ax._hold = washold

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/__init__.pyc in inner(ax, *args, **kwargs)
   1890                     warnings.warn(msg % (label_namer, func.__name__),
   1891                                   RuntimeWarning, stacklevel=2)
-> 1892             return func(ax, *args, **kwargs)
   1893         pre_doc = inner.__doc__
   1894         if pre_doc is None:

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.pyc in plot(self, *args, **kwargs)
   1404         kwargs = cbook.normalize_kwargs(kwargs, _alias_map)
   1405 
-> 1406         for line in self._get_lines(*args, **kwargs):
   1407             self.add_line(line)
   1408             lines.append(line)

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _grab_next_args(self, *args, **kwargs)
    405                 return
    406             if len(remaining) <= 3:
--> 407                 for seg in self._plot_args(remaining, kwargs):
    408                     yield seg
    409                 return

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _plot_args(self, tup, kwargs)
    383             x, y = index_of(tup[-1])
    384 
--> 385         x, y = self._xy_from_xy(x, y)
    386 
    387         if self.command == 'plot':

/Users/rohin/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in _xy_from_xy(self, x, y)
    242         if x.shape[0] != y.shape[0]:
    243             raise ValueError("x and y must have same first dimension, but "
--> 244                              "have shapes {} and {}".format(x.shape, y.shape))
    245         if x.ndim > 2 or y.ndim > 2:
    246             raise ValueError("x and y can be no greater than 2-D, but have "

ValueError: x and y must have same first dimension, but have shapes (3966,) and (3967,)



In [45]:

    
hmqdat1.sort('R')
hmqdat1









    Out[45]:




<Table length=1305226>

deltaZ deltatheta Zdeltatheta R alpha thetadd DD
float64 float64 float64 float64 float64 float64 float64
0.0 5e-06 1e-06 1e-06 1.570796 0.538925 4.0
0.0 7e-06 2e-06 2e-06 1.570796 0.538925 4.0
0.0 1.3e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.1e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.2e-05 3e-06 3e-06 1.570796 0.538925 4.0
0.0 1.7e-05 4e-06 4e-06 1.570796 0.538925 9.0
0.0 2.2e-05 5e-06 5e-06 1.570796 0.538925 4.0
0.0 2e-05 5e-06 5e-06 1.570796 0.538925 169.0
0.0 2.1e-05 5e-06 5e-06 1.570796 0.538925 4.0
0.0 2.2e-05 5e-06 5e-06 1.570796 0.538925 4.0
... ... ... ... ... ... ...
0.02 0.083981 0.019995 0.02828 0.785262 0.891568 1.0
0.02 0.083978 0.019994 0.02828 0.785239 2.044048 1.0
0.02 0.083984 0.019995 0.028281 0.785277 0.792125 1.0
0.02 0.083987 0.019996 0.028281 0.785292 0.501326 1.0
0.02 0.083989 0.019996 0.028282 0.785309 0.111224 1.0
0.02 0.083988 0.019996 0.028282 0.785301 2.064619 2.0
0.02 0.083999 0.019999 0.028283 0.785368 0.599163 1.0
0.02 0.083998 0.019998 0.028283 0.785359 0.334718 1.0
0.02 0.084 0.019999 0.028284 0.785372 0.704685 2.0
0.02 0.084 0.019999 0.028284 0.785372 0.909261 1.0



In [50]:

    
DD_R=hist(hmqdat1['R'],bins='blocks')
plt.savefig("DD_based_onRhist.pdf")
plt.savefig("DD_based_onRhist.jpeg")



In [47]:

    
DD_R_1=histogram(hmqdat1['R'],bins='blocks')



In [48]:

    
help(hist)









    



Help on function hist in module astropy.visualization.hist:

hist(x, bins=10, ax=None, **kwargs)
    Enhanced histogram function
    
    This is a histogram function that enables the use of more sophisticated
    algorithms for determining bins.  Aside from the ``bins`` argument allowing
    a string specified how bins are computed, the parameters are the same
    as pylab.hist().
    
    This function was ported from astroML: http://astroML.org/
    
    Parameters
    ----------
    x : array_like
        array of data to be histogrammed
    
    bins : int or list or str (optional)
        If bins is a string, then it must be one of:
    
        - 'blocks' : use bayesian blocks for dynamic bin widths
    
        - 'knuth' : use Knuth's rule to determine bins
    
        - 'scott' : use Scott's rule to determine bins
    
        - 'freedman' : use the Freedman-diaconis rule to determine bins
    
    ax : Axes instance (optional)
        specify the Axes on which to draw the histogram.  If not specified,
        then the current active axes will be used.
    
    **kwargs :
        other keyword arguments are described in ``plt.hist()``.
    
    Notes
    -----
    Return values are the same as for ``plt.hist()``
    
    See Also
    --------
    astropy.stats.histogram



In [49]:

    
DD_R_1









    Out[49]:





(array([   949,     52,    210,    605,   2806,    524,   5276,    387,
          7940,  11464,  15021,  19270,  23822,  28337,  34104,  40473,
         47122,  27490,  25806,  62260,  70046,  66430,  11750,  57827,
         30518,  98591,  76127,  32001, 129991,  13201,  17653,  20958,
         23873,  25857,  26699,  30252,  29896,  28590,  28502,  26619,
         25662,  22361,  19863,  16113,  12192,   7394,   2342]),
 array([  1.00000000e-06,   9.99500000e-04,   1.00150000e-03,
          1.02650000e-03,   1.16250000e-03,   1.99950000e-03,
          2.06350000e-03,   2.99950000e-03,   3.03650000e-03,
          3.99850000e-03,   4.99950000e-03,   5.99950000e-03,
          6.99750000e-03,   7.99950000e-03,   8.99850000e-03,
          9.99750000e-03,   1.09975000e-02,   1.19995000e-02,
          1.25205000e-02,   1.29895000e-02,   1.39985000e-02,
          1.49995000e-02,   1.58485000e-02,   1.59925000e-02,
          1.66555000e-02,   1.69955000e-02,   1.80045000e-02,
          1.87135000e-02,   1.90025000e-02,   2.00995000e-02,
          2.02225000e-02,   2.03965000e-02,   2.06155000e-02,
          2.08785000e-02,   2.11845000e-02,   2.15175000e-02,
          2.19255000e-02,   2.23645000e-02,   2.28235000e-02,
          2.33225000e-02,   2.38395000e-02,   2.44175000e-02,
          2.49905000e-02,   2.56105000e-02,   2.62465000e-02,
          2.69065000e-02,   2.75855000e-02,   2.82840000e-02]))



In [51]:

    
DD_R=hist(hmqdat1['R'],bins=1000)
plt.savefig("./images/DD_based_onRhist1000.pdf")
plt.savefig("./images/DD_based_onRhist1000.jpeg")



In [52]:

    
DD_R=hist(hmqdat1['R'],bins='knuth')
plt.savefig("./images/DD_based_onRhistknuth.pdf")
plt.savefig("./images/DD_based_onRhistknuth.jpeg")



In [53]:

    
DD_R









    Out[53]:





(array([   119.,    103.,    108.,    126.,    107.,    120.,    131.,
           131.,    674.,    484.,    381.,    429.,    405.,    429.,
           441.,    419.,    856.,    729.,    688.,    635.,    710.,
           716.,    727.,    709.,   1067.,   1020.,   1054.,    995.,
          1009.,   1045.,   1054.,   1042.,   1460.,   1359.,   1415.,
          1419.,   1414.,   1457.,   1477.,   1389.,   1876.,   1779.,
          1819.,   1897.,   1915.,   1861.,   1854.,   1914.,   2410.,
          2404.,   2333.,   2414.,   2380.,   2416.,   2357.,   2464.,
          2894.,   2850.,   2958.,   2925.,   2982.,   2981.,   3023.,
          3016.,   3511.,   3507.,   3374.,   3568.,   3566.,   3532.,
          3527.,   3529.,   4129.,   4160.,   4123.,   4309.,   4306.,
          4312.,   4309.,   4271.,   4813.,   4963.,   4950.,   5022.,
          5079.,   5104.,   5063.,   5101.,   5719.,   5817.,   5727.,
          5841.,   5897.,   5858.,   5825.,   5990.,   6349.,   6617.,
          6572.,   6607.,   6765.,   6840.,   6820.,   6848.,   7499.,
          7512.,   7699.,   7640.,   7634.,   7625.,   7837.,   7761.,
          8327.,   8521.,   8731.,   8597.,   8789.,   8938.,   8657.,
          8836.,   9447.,   9696.,   9631.,   9861.,   9799.,   9816.,
          9643.,  10188.,  10542.,  10789.,  11014.,  10844.,  10924.,
         10898.,  11184.,  11215.,  11662.,  11931.,  12176.,  12152.,
         12168.,  12148.,  12358.,  12278.,  12877.,  13122.,  13337.,
         13436.,  13462.,  13574.,  13706.,  13822.,  14133.,  14768.,
         14611.,  14628.,  14965.,  14541.,  14909.,  14984.,  14890.,
         13780.,  12615.,  12517.,  11941.,  11517.,  11434.,  10865.,
         10614.,  10537.,   9865.,  10217.,   9716.,   9122.,   9196.,
          9293.,   8387.,   8497.,   8524.,   8133.,   7588.,   7860.,
          7797.,   7360.,   7092.,   7030.,   7186.,   6486.,   6374.,
          6366.,   6482.,   5831.,   5480.,   5487.,   5599.,   5585.,
          4801.,   4788.,   4954.,   4881.,   4574.,   4034.,   3871.,
          4077.,   3899.,   3631.,   3127.,   3196.,   3232.,   3090.,
          2861.,   2249.,   2261.,   2245.,   2317.,   2373.,   1354.,
          1311.,   1268.,   1415.,   1404.,    833.,    431.,    416.,
           422.,    404.,    407.]),
 array([  1.00000000e-06,   1.25594714e-04,   2.50189427e-04,
          3.74784141e-04,   4.99378855e-04,   6.23973568e-04,
          7.48568282e-04,   8.73162996e-04,   9.97757709e-04,
          1.12235242e-03,   1.24694714e-03,   1.37154185e-03,
          1.49613656e-03,   1.62073128e-03,   1.74532599e-03,
          1.86992070e-03,   1.99451542e-03,   2.11911013e-03,
          2.24370485e-03,   2.36829956e-03,   2.49289427e-03,
          2.61748899e-03,   2.74208370e-03,   2.86667841e-03,
          2.99127313e-03,   3.11586784e-03,   3.24046256e-03,
          3.36505727e-03,   3.48965198e-03,   3.61424670e-03,
          3.73884141e-03,   3.86343612e-03,   3.98803084e-03,
          4.11262555e-03,   4.23722026e-03,   4.36181498e-03,
          4.48640969e-03,   4.61100441e-03,   4.73559912e-03,
          4.86019383e-03,   4.98478855e-03,   5.10938326e-03,
          5.23397797e-03,   5.35857269e-03,   5.48316740e-03,
          5.60776211e-03,   5.73235683e-03,   5.85695154e-03,
          5.98154626e-03,   6.10614097e-03,   6.23073568e-03,
          6.35533040e-03,   6.47992511e-03,   6.60451982e-03,
          6.72911454e-03,   6.85370925e-03,   6.97830396e-03,
          7.10289868e-03,   7.22749339e-03,   7.35208811e-03,
          7.47668282e-03,   7.60127753e-03,   7.72587225e-03,
          7.85046696e-03,   7.97506167e-03,   8.09965639e-03,
          8.22425110e-03,   8.34884581e-03,   8.47344053e-03,
          8.59803524e-03,   8.72262996e-03,   8.84722467e-03,
          8.97181938e-03,   9.09641410e-03,   9.22100881e-03,
          9.34560352e-03,   9.47019824e-03,   9.59479295e-03,
          9.71938767e-03,   9.84398238e-03,   9.96857709e-03,
          1.00931718e-02,   1.02177665e-02,   1.03423612e-02,
          1.04669559e-02,   1.05915507e-02,   1.07161454e-02,
          1.08407401e-02,   1.09653348e-02,   1.10899295e-02,
          1.12145242e-02,   1.13391189e-02,   1.14637137e-02,
          1.15883084e-02,   1.17129031e-02,   1.18374978e-02,
          1.19620925e-02,   1.20866872e-02,   1.22112819e-02,
          1.23358767e-02,   1.24604714e-02,   1.25850661e-02,
          1.27096608e-02,   1.28342555e-02,   1.29588502e-02,
          1.30834449e-02,   1.32080396e-02,   1.33326344e-02,
          1.34572291e-02,   1.35818238e-02,   1.37064185e-02,
          1.38310132e-02,   1.39556079e-02,   1.40802026e-02,
          1.42047974e-02,   1.43293921e-02,   1.44539868e-02,
          1.45785815e-02,   1.47031762e-02,   1.48277709e-02,
          1.49523656e-02,   1.50769604e-02,   1.52015551e-02,
          1.53261498e-02,   1.54507445e-02,   1.55753392e-02,
          1.56999339e-02,   1.58245286e-02,   1.59491233e-02,
          1.60737181e-02,   1.61983128e-02,   1.63229075e-02,
          1.64475022e-02,   1.65720969e-02,   1.66966916e-02,
          1.68212863e-02,   1.69458811e-02,   1.70704758e-02,
          1.71950705e-02,   1.73196652e-02,   1.74442599e-02,
          1.75688546e-02,   1.76934493e-02,   1.78180441e-02,
          1.79426388e-02,   1.80672335e-02,   1.81918282e-02,
          1.83164229e-02,   1.84410176e-02,   1.85656123e-02,
          1.86902070e-02,   1.88148018e-02,   1.89393965e-02,
          1.90639912e-02,   1.91885859e-02,   1.93131806e-02,
          1.94377753e-02,   1.95623700e-02,   1.96869648e-02,
          1.98115595e-02,   1.99361542e-02,   2.00607489e-02,
          2.01853436e-02,   2.03099383e-02,   2.04345330e-02,
          2.05591278e-02,   2.06837225e-02,   2.08083172e-02,
          2.09329119e-02,   2.10575066e-02,   2.11821013e-02,
          2.13066960e-02,   2.14312907e-02,   2.15558855e-02,
          2.16804802e-02,   2.18050749e-02,   2.19296696e-02,
          2.20542643e-02,   2.21788590e-02,   2.23034537e-02,
          2.24280485e-02,   2.25526432e-02,   2.26772379e-02,
          2.28018326e-02,   2.29264273e-02,   2.30510220e-02,
          2.31756167e-02,   2.33002115e-02,   2.34248062e-02,
          2.35494009e-02,   2.36739956e-02,   2.37985903e-02,
          2.39231850e-02,   2.40477797e-02,   2.41723744e-02,
          2.42969692e-02,   2.44215639e-02,   2.45461586e-02,
          2.46707533e-02,   2.47953480e-02,   2.49199427e-02,
          2.50445374e-02,   2.51691322e-02,   2.52937269e-02,
          2.54183216e-02,   2.55429163e-02,   2.56675110e-02,
          2.57921057e-02,   2.59167004e-02,   2.60412952e-02,
          2.61658899e-02,   2.62904846e-02,   2.64150793e-02,
          2.65396740e-02,   2.66642687e-02,   2.67888634e-02,
          2.69134581e-02,   2.70380529e-02,   2.71626476e-02,
          2.72872423e-02,   2.74118370e-02,   2.75364317e-02,
          2.76610264e-02,   2.77856211e-02,   2.79102159e-02,
          2.80348106e-02,   2.81594053e-02,   2.82840000e-02]),
 <a list of 227 Patch objects>)



In [54]:

    
help(np.correlate)









    



Help on function correlate in module numpy.core.numeric:

correlate(a, v, mode='valid')
    Cross-correlation of two 1-dimensional sequences.
    
    This function computes the correlation as generally defined in signal
    processing texts::
    
        c_{av}[k] = sum_n a[n+k] * conj(v[n])
    
    with a and v sequences being zero-padded where necessary and conj being
    the conjugate.
    
    Parameters
    ----------
    a, v : array_like
        Input sequences.
    mode : {'valid', 'same', 'full'}, optional
        Refer to the `convolve` docstring.  Note that the default
        is 'valid', unlike `convolve`, which uses 'full'.
    old_behavior : bool
        `old_behavior` was removed in NumPy 1.10. If you need the old
        behavior, use `multiarray.correlate`.
    
    Returns
    -------
    out : ndarray
        Discrete cross-correlation of `a` and `v`.
    
    See Also
    --------
    convolve : Discrete, linear convolution of two one-dimensional sequences.
    multiarray.correlate : Old, no conjugate, version of correlate.
    
    Notes
    -----
    The definition of correlation above is not unique and sometimes correlation
    may be defined differently. Another common definition is::
    
        c'_{av}[k] = sum_n a[n] conj(v[n+k])
    
    which is related to ``c_{av}[k]`` by ``c'_{av}[k] = c_{av}[-k]``.
    
    Examples
    --------
    >>> np.correlate([1, 2, 3], [0, 1, 0.5])
    array([ 3.5])
    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "same")
    array([ 2. ,  3.5,  3. ])
    >>> np.correlate([1, 2, 3], [0, 1, 0.5], "full")
    array([ 0.5,  2. ,  3.5,  3. ,  0. ])
    
    Using complex sequences:
    
    >>> np.correlate([1+1j, 2, 3-1j], [0, 1, 0.5j], 'full')
    array([ 0.5-0.5j,  1.0+0.j ,  1.5-1.5j,  3.0-1.j ,  0.0+0.j ])
    
    Note that you get the time reversed, complex conjugated result
    when the two input sequences change places, i.e.,
    ``c_{va}[k] = c^{*}_{av}[-k]``:
    
    >>> np.correlate([0, 1, 0.5j], [1+1j, 2, 3-1j], 'full')
    array([ 0.0+0.j ,  3.0+1.j ,  1.5+1.5j,  1.0+0.j ,  0.5+0.5j])



In [55]:

    
def AutoCorrelation(x):
    #x = np.asarray(x)
    y = x-x.mean()
    result = np.correlate(y, y, mode='full')
    result = result[len(result)//2:]
    result /= result[0]
    return result



In [65]:

    
dat=DD_R[0]



In [66]:

    
dd=AutoCorrelation(dat)



In [70]:

    
plt.plot(DD_R[1][1:len(DD_R[1])],dd)









    Out[70]:





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



In [63]:

    
len(DD_R[1])









    Out[63]:





228



In [64]:

    
len(dat)









    Out[64]:





227



In [72]:

    
print time.time()-134432.3









    



1492880139.73



In [ ]:

deltaZ	deltatheta	Zdeltatheta	R	alpha	thetadd	DD
float64	float64	float64	float64	float64	float64	float64
0.013	0.051491	0.012259	0.017869	0.756076	0.052368	1.0
0.013	0.040922	0.009743	0.016246	0.64315	0.041069	1.0
0.008	0.062346	0.014843	0.016862	1.07647	1.878369	1.0
0.009	0.068156	0.016227	0.018556	1.064399	1.893666	1.0
0.01	0.067917	0.01617	0.019012	1.016931	1.907749	1.0
0.001	0.051615	0.012289	0.012329	1.489599	0.992046	1.0
0.007	0.083981	0.019994	0.021184	1.234034	1.219094	1.0
0.017	0.082431	0.019625	0.025964	0.856957	0.604861	1.0
0.018	0.070245	0.016724	0.02457	0.74867	1.063174	1.0
0.001	0.040053	0.009536	0.009588	1.466311	2.505325	1.0
...	...	...	...	...	...	...
0.0	0.079949	0.019034	0.019034	1.570796	2.53098	1.0
0.0	0.053744	0.012796	0.012796	1.570796	0.988719	2.0
0.0	0.075223	0.017909	0.017909	1.570796	2.53098	2.0
0.0	0.046556	0.011084	0.011084	1.570796	0.538925	1.0
0.0	0.068345	0.016272	0.016272	1.570796	1.952672	1.0
0.0	0.041424	0.009862	0.009862	1.570796	2.466464	1.0
0.0	0.033226	0.007911	0.007911	1.570796	0.988719	1.0
0.0	0.061433	0.014626	0.014626	1.570796	2.551203	1.0
0.0	0.067922	0.016171	0.016171	1.570796	2.551203	1.0
0.0	0.061536	0.014651	0.014651	1.570796	0.988719	1.0