def calculateFDunc(self):
        #Calculates the uncertainty of the FFT according to:
        #   - J. M. Fornies-Marquina, J. Letosa, M. Garcia-Garcia, J. M. Artacho, "Error Propagation for the transformation of time domain into frequency domain", IEEE Trans. Magn, Vol. 33, No. 2, March 1997, pp. 1456-1459
        #return asarray _tdData
        #Assumes tha the amplitude of each time sample is statistically independent from the amplitude of the other time
        #samples

        # Calculates uncertainty of the real and imaginary part of the FFT and ther covariance
        unc_E_real = []
        unc_E_imag = []
        cov = []
        for f in self.getfreqs():
            unc_E_real.append(py.sum((py.cos(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            unc_E_imag.append(py.sum((py.sin(2*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX())**2))
            cov.append(-0.5*sum(py.sin(4*py.pi*f*self._tdData.getTimes())*self._tdData.getUncEX()**2))
        
        unc_E_real = py.sqrt(py.asarray(unc_E_real))
        unc_E_imag = py.sqrt(py.asarray(unc_E_imag))
        cov = py.asarray(cov)
        
        # Calculates the uncertainty of the modulus and phase of the FFT
        unc_E_abs = py.sqrt((self.getFReal()**2*unc_E_real**2+self.getFImag()**2*unc_E_imag**2+2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**2)
        unc_E_ph = py.sqrt((self.getFImag()**2*unc_E_real**2+self.getFReal()**2*unc_E_imag**2-2*self.getFReal()*self.getFImag()*cov)/self.getFAbs()**4)
        
        t=py.column_stack((self.getfreqs(),unc_E_real,unc_E_imag,unc_E_abs,unc_E_ph))
        return self.getcroppedData(t)  

def get_first_init(x0,epsilon,N):
    x_new = pl.copy(x0)
    print('getting the first initial condition')
    print('fiducial initial: '+str(x0))
    # multi particle array layout [nth particle v, (n-1)th particle v , ..., 0th v, nth particle x, x, ... , 0th particle x]
    
    # we will use a change of coordinates to get the location of the particle relative to x0. First
    # we just find some random point a distace epsilon from the origin.
    # need 2DN random angles 
    angle_arr = pl.array([])
    purturbs = pl.array([])

    # This is just an n-sphere 
    for i in range(2*N):
        angle_arr = pl.append(angle_arr,random.random()*2.0*pl.pi)
        cur_purt = epsilon
        for a,b in enumerate(angle_arr[:-1]):
            cur_purt *= pl.sin(b)
        if i == (2*N-1):
            cur_purt = pl.sin(angle_arr[i])
        else:
            cur_purt = pl.cos(angle_arr[i])

        purturbs = pl.append(purturbs,cur_purt)

    print('sqrt of sum of squars should be epsilon -> is it? --> ' +str(pl.sqrt(pl.dot(purturbs,purturbs))))
    print('len(purturbs) == 2N ? ' +str(len(purturbs)==(2*N)))

    return x_new+purturbs

def findcurve(psi1,psi2,n=3,nn_fit=4,nn_out=100):
    '''
    Function to find the elastica curve for start and end orientations
    psi1 and psi2. It finds the best curve across all directions from start
    and end, i.e. the direction independent elastica curve.
    
    Inputs
    ------------
    psi1,psi2: start and end orientations.
    n:     degree of estimation polynomial.
    nn:    number of points on the curve.
             - nn_fit: for fittin purposes
             - nn_out: for the output
    
    Outputs
    ------------
    Returns a tuple (s,psi). 
    s:   points on the curve.
    psi: curvature of the curve as a function of s.
    E:   curvature energy of the curve
    '''
    # 
    
    # define the starting conditions
    a0 = pl.zeros(n+1) 
    
    # Set a high energy: 
    E_best = 10000  
    
    # and predfine output curve
    s       = pl.linspace(0,1,nn_out) # points on the curve
    psi_out = pl.zeros(nn_out)        # curvature at points in curve
    
    
    # across all the start and end directions find the curve with the lowest energy    
    for dpsi1 in (-pl.pi,0,pl.pi):
        for dpsi2 in (-pl.pi,0,pl.pi):
            # For the starting variables,
            # the first two polygon variables can be estimated from the Sharon paper derivation
            # For different starting variables the solution can be hard to find            
            a0[-2] = 4*(   pl.arcsin(- (pl.sin(psi1+dpsi1)+ pl.sin(psi2+dpsi2))/4)    -(psi1+dpsi1+psi2+dpsi2)/2       )
            a0[-1] = 2*a0[-2]/pl.cos( (psi1+dpsi1+psi2+dpsi2)/2 + a0[-2]/4  )               
            
            # find the best variables to minimize the elastica energy
            fit = fsolve(errors,a0,args=(psi1+dpsi1,psi2+dpsi2,nn_fit))
    
            # find the curve and its derivative for the fitted variables
            a    = fit[:-1]
            psi  = Psi(a,s,psi1+dpsi1,psi2+dpsi2)
            dpsi = dPsi(a,s,psi1+dpsi1,psi2+dpsi2)
    
            # find the energy of this curve
            E = sum(dpsi**2)*s[1]
            
            # check against the lowest energy
            if E_best > E:
                E_best = E
                psi_out[:] = pl.copy(psi)    
    
    return (s,psi_out,E_best)

    def get_touchdown(self, t, y, params):
        """
        Compute the touchdown position of the leg w.r.t. CoM velocity
        
        :args:
            t (float): time
            y (6x float): state of the CoM
            params (4x float): leg parameter: stiffness, l0, alpha, beta
        
        :returns:
            [xFoot, yFoot, zFoot] the position of the leg tip
        """
        k, l0, alpha, beta = params
        vx, vz = y[3], y[5]
        
        a_v_com = -arctan2(vz, vx) # correct with our coordinate system
        #for debugging
        #print "v_com_angle:", a_v_com * 180. / pi
                
        xf = y[0] + l0 * cos(alpha) * cos(beta + a_v_com)
        yf = y[1] - l0 * sin(alpha)
        zf = y[2] - l0 * cos(alpha) * sin(beta + a_v_com)

        #for debugging
        #print "foot: %2.3f,%2.3f,%2.3f," % ( xf,yf, zf)
        
        return array([xf, yf, zf])

def tests():
	x = pylab.arange(0.0, 2*pylab.pi, 0.01)
	list_y = (pylab.sin(x),pylab.sin(2*x))
	plot_and_format((x,), (list_y[0],))
	exportplot('test/test1_one_curve.png')
	pylab.clf()
	list_x = (x, x)
	plot_and_format(list_x, list_y)
	exportplot('test/test2_two_curves.png')
	pylab.clf()
	list_format = ('k-', 'r--')
	plot_and_format(list_x, list_y, list_format=list_format)
	exportplot('test/test3_two_curves_formatting.png')
	pylab.clf()
	plot_and_format(list_x, list_y, list_format=list_format, xlabel='hello x axis')
	exportplot('test/test4_two_curves_formatting_xlab.png')
	pylab.clf()
	plot_and_format(list_x, list_y, list_format=list_format, legend=['sin($x$)', 'sin($2x$)'])
	exportplot('test/test5_two_curves_formatting_legend.png')
	pylab.clf()
	plot_and_format(list_x, list_y, list_format=list_format, xticks={'ticks': [0, pylab.pi, 2*pylab.pi], 'labels':['0', '$\pi$', '$2\pi$']})
	exportplot('test/test6_two_curves_formatting_xticks.png')
	pylab.clf()
	plot_and_format(list_x, list_y, list_format=list_format, xticks={'ticks': [0, pylab.pi, 2*pylab.pi], 'labels':['0', '$\pi$', '$2\pi$']}, xlim=[0,2*pylab.pi])
	exportplot('test/test7_two_curves_formatting_xticks_xlim.png')
	pylab.clf()

def specgram_demo():
   '''
   the demo in matplotlib. But calls
   interactive.specgram
   '''
   from pylab import arange, sin, where, logical_and, randn, pi

   dt = 0.0005
   t = arange(0.0, 20.0, dt)
   s1 = sin(2*pi*100*t)
   s2 = 2*sin(2*pi*400*t)

   # create a transient "chirp"
   mask = where(logical_and(t>10, t<12), 1.0, 0.0)
   s2 = s2 * mask

   # add some noise into the mix
   nse = 0.01*randn(len(t))

   x = s1 + s2 + nse # the signal
   NFFT = 1024       # the length of the windowing segments
   Fs = int(1.0/dt)  # the sampling frequency

   from ifigure.interactive import figure, specgram, nsec, plot, isec, clog, hold

   figure()
   hold(True)
   nsec(2)
   isec(0)
   plot(t, x)
   isec(1)
   specgram(x, NFFT=NFFT, Fs=Fs, noverlap=900)
   clog()

def hillshade(data,scale=10.0,azdeg=165.0,altdeg=45.0):
  ''' 
    This code thanks to Ran Novitsky Nof
  http://rnovitsky.blogspot.co.uk/2010/04/using-hillshade-image-as-intensity.html
  Repeated here to make my cyclopean uk_map code prettier.

  convert data to hillshade based on matplotlib.colors.LightSource class.
    input:
         data - a 2-d array of data
         scale - scaling value of the data. higher number = lower gradient
         azdeg - where the light comes from: 0 south ; 90 east ; 180 north ;
                      270 west
         altdeg - where the light comes from: 0 horison ; 90 zenith
    output: a 2-d array of normalized hilshade
  '''
  
  from pylab import pi, gradient, arctan, hypot, arctan2, sin, cos
  # convert alt, az to radians
  az = azdeg*pi/180.0
  alt = altdeg*pi/180.0
  # gradient in x and y directions
  dx, dy = gradient(data/float(scale))
  slope = 0.5*pi - arctan(hypot(dx, dy))
  aspect = arctan2(dx, dy)
  intensity = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(-az - aspect - 0.5*pi)
  intensity = (intensity - intensity.min())/(intensity.max() - intensity.min())
  return intensity

def main():
    # Create the grid
    x = arange(-100, 101)
    y = arange(-100, 101)

    # Create the meshgrid
    Y, X = meshgrid(x, y)
    A = 1
    B = 2
    V = 6*pi / 201
    W = 4*pi / 201
    F = A*sin(V*X) + B*cos(W*Y)
    Fx = V*A*cos(V*X)
    Fy = W*B*-sin(W*Y)

    # Show the images
    show_image(F)
    show_image(Fx)
    show_image(Fy)

    # Create the grid for the quivers
    xs = arange(-100, 101, 10)
    ys = arange(-100, 101, 10)

    # Here we determine the direction of the quivers
    Ys, Xs = meshgrid(ys, xs)
    FFx = V*A*cos(V*Xs)
    FFy = W*B*-sin(W*Ys)

    # Draw the quivers and the image
    clf()
    imshow(F, cmap=cm.gray, extent=(-100, 100, -100, 100))
    quiver(ys, xs, -FFy, FFx, color='red')
    show()

def f_mdl(mdl, phi, maxord):
    """
    given a periodic function "mdl" consisting of n data points (ranging from  [0,2pi)),
    a fourier model of order maxord is computed for all phases in phi

    :args:
        mdl (1-by-k array): datapoints of the reference model, ranging from
           [0, 2pi). Length in datapoints is arbitrary.
        phi (1-by-n array): phases at which to evaluate the model
        maxord (int): order of the Fourier model to compute

    :returns:
       mdl_val (1-by-n array): value of the Fourier model obtained from mdl for
           the given phases phi
    
    """
    spec_fy = fft.fft(mdl)
    as_, ac_ = spec_fy.imag, spec_fy.real
    sigout = zeros(len(phi))
    for order in range(maxord):
        sigout -= sin(order * phi) * as_[order]
        sigout += cos(order * phi) * ac_[order]    
        sigout += cos(order * phi) * ac_[-order]
        sigout += sin(order * phi) * as_[-order]
    sigout /= len(mdl)
    return sigout

def hillshade(data, scale=10.0, azdeg=165.0, altdeg=45.0):
    '''
    Convert data to hillshade based on matplotlib.colors.LightSource class.

    Args:
        data - a 2-d array of data
        scale - scaling value of the data. higher number = lower gradient
        azdeg - where the light comes from: 0 south ; 90 east ; 180 north ;
                        270 west
        altdeg - where the light comes from: 0 horison ; 90 zenith

    Returns:
        a 2-d array of normalized hilshade
    '''
    # convert alt, az to radians
    az = azdeg*pi/180.0
    alt = altdeg*pi/180.0
    # gradient in x and y directions
    dx, dy = gradient(data/float(scale))
    slope = 0.5*pi - arctan(hypot(dx, dy))
    aspect = arctan2(dx, dy)
    az = -az - aspect - 0.5*pi
    intensity = sin(alt)*sin(slope) + cos(alt)*cos(slope)*cos(az)
    mi, ma = intensity.min(), intensity.max()
    intensity = (intensity - mi)/(ma - mi)
    return intensity

def haversine (latlong1, latlong2, r):

    deltaLatlong = latlong1 - latlong2
    
    dLat = deltaLatlong[0]
    dLon = deltaLatlong[1]

    lat1 = latlong1[0]
    lat2 = latlong2[0]

    a = (sin (dLat/2) * sin (dLat/2) +
         sin (dLon/2) * sin (dLon/2) * cos (lat1) * cos (lat2))
    c = 2 * arctan2 (sqrt (a), sqrt (1-a))
    d = r * c

    # initial bearing
    y = sin (dLon) * cos (lat2)
    x = (cos (lat1)*sin (lat2) -
         sin (lat1)*cos (lat2)*cos (dLon))
    b1 = arctan2 (y, x);

    # final bearing
    dLon = -dLon
    dLat = -dLat
    tmp = lat1
    lat1 = lat2
    lat2 = tmp
    y = sin (dLon) * cos (lat2)
    x = (cos (lat1) * sin (lat2) - 
         sin (lat1) * cos (lat2) * cos (dLon))
    b2 = arctan2 (y, x)
    b2 = mod ((b2 + pi), 2*pi)

    return (d, b1, b2)

def sortAnglesCW(t1,t2):

    """
    Sort angles so that t2>t1 in a clockwise sense
    idea from `StackOverflow <http://stackoverflow.com/questions/242404/sort-four-points-in-clockwise-order>`_
    more description: `SoftSurfer <http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm>`_

    If the signed area of the triangle formed between the points on a unit circle with angles t1 and t2
    and the origin is positive, the angles are sorted counterclockwise. Otherwise, the angles
    are sorted in a counter-clockwise manner.  Here we want the angles to be sorted CCW, so
    if area is negative, swap angles
    
    Area obtained from the cross product of a vector from origin 
    to 1 and a vector to point 2, so use right hand rule to get 
    sign of cross product with unit length
    """

    while (cos(t1)*sin(t2)-sin(t1)*cos(t2)>0):
        ##Swap angles
        temp=t1;
        t1=t2;
        t2=temp;
    #Make t1 between 0 and 2pi
    while (t1<0 or t1> 2.0*pi):
        if t1>2.0*pi:
            t1=t1-2*pi
        else:
            t1=t1+2*pi
    #Want t2 to be less than t1, but no more than 2*pi less
    while (t2<t1 and t1-t2>2*pi):
        t2=t2+2*pi
    while (t2>t1):
        t2=t2-2*pi
    return (t1,t2)

def haversine(location1, location2=None):  # calculates great circle distance
    __doc__ = """Returns the great circle distance of the given
    coordinates.
    
    INPUT:  location1 = ((lat1, lon1), ..., n(lat1, lon1))
           *location2 = ((lat2, lon2), ..., n(lat2, lon2))
           *if location2 is not given a square matrix of distances
             for location1 will be put out
    OUTPUT: distance in km
            (dist1  ...  ndist
              :            : 
             ndist1 ...  ndist)
            shape will depend on the input
    METHOD: a = sin(dLat / 2) * sin(dLat / 2) + 
                sin(dLon / 2) * sin(dLon / 2) * 
                cos(lat1) * cos(lat2)
            c = 2 * arctan2(sqrt(a), sqrt(1 - a))
            d = R * c
            
            where R is the earth's radius (6371 km)
            and d is the distance in km"""
    
    from itertools import product, combinations
    from pylab import   deg2rad, sin, cos, arctan2, \
                        meshgrid, sqrt, array, arange
    
    if location2: 
        location1 = array(location1, ndmin=2)
        location2 = array(location2, ndmin=2)
    elif location2 is None:
        location1 = array(location1, ndmin=2)
        location2 = location1.copy()
    
    # get all combinations using indicies
    ind1 = arange(location1.shape[0])
    ind2 = arange(location2.shape[0])
    ind  = array(list(product(ind1, ind2)))
    
    # using combination inds to get lats and lons
    lat1, lon1 = location1[ind[:,0]].T
    lat2, lon2 = location2[ind[:,1]].T
    
    # setting up variables for haversine
    R = 6371.
    dLat = deg2rad(lat2 - lat1)
    dLon = deg2rad(lon2 - lon1)
    lat1 = deg2rad(lat1)
    lat2 = deg2rad(lat2)
    
    # haversine formula
    a = sin(dLat / 2) * sin(dLat / 2) + \
        sin(dLon / 2) * sin(dLon / 2) * \
        cos(lat1) * cos(lat2)
    c = 2 * arctan2(sqrt(a), sqrt(1 - a))
    d = R * c
    
    # reshape accodring to the input
    D = d.reshape(location1.shape[0], location2.shape[0])
    
    return D

def measureDistance(lat1, lon1, lat2, lon2):
    R = 6383.137 # Radius of earth at Chajnantor aprox. in KM
    dLat = (lat2 - lat1) * np.pi / 180.
    dLon = (lon2 - lon1) * np.pi / 180.
    a = pl.sin(dLat/2.) * pl.sin(dLat/2.) + pl.cos(lat1 * np.pi / 180.) * pl.cos(lat2 * np.pi / 180.) * pl.sin(dLon/2.) * pl.sin(dLon/2.)
    c = 2. * atan2(pl.sqrt(a), pl.sqrt(1-a))
    d = R * c
    return d * 1000. # meters

 def rotated_rectangle(x0,y0,w,h,rot):
     
     x = np.array([-w/2,w/2,w/2,-w/2,-w/2])
     y = np.array([-h/2,-h/2,h/2,h/2,-h/2])
     
     xrot = x*cos(rot)-y*sin(rot)
     yrot = x*sin(rot)+y*cos(rot) 
     
     return xrot+x0, yrot+y0

def random_euler_angles():
    r1,r2,r3 = pylab.random(3)
    q1 = pylab.sqrt(1.0-r1)*pylab.sin(2.0*pylab.pi*r2)
    q2 = pylab.sqrt(1.0-r1)*pylab.cos(2.0*pylab.pi*r2)
    q3 = pylab.sqrt(r1)*pylab.sin(2.0*pylab.pi*r3)
    q4 = pylab.sqrt(r1)*pylab.cos(2.0*pylab.pi*r3)
    phi = math.atan2(2.0*(q1*q2+q3*q4), 1.0-2.0*(q2**2+q3**2))
    theta = math.asin(2.0*(q1*q3-q4*q2))
    psi = math.atan2(2.0*(q1*q4+q2*q3), 1.0-2.0*(q3**2+q4**2))
    return [phi,theta,psi]

def plot_trace(X, scale=1., angle=0.):
    fig = pl.figure(figsize=(12,4.75))
    
    ax1 = fig.add_subplot(1, 2, 1)
    # plot boundary
    t = pl.arange(0,2*pl.pi,.01)
    ax1.plot(pl.cos(angle)*pl.cos(t) - pl.sin(angle)*pl.sin(t)/scale, pl.cos(angle)*pl.sin(t)/scale + pl.sin(angle)*pl.cos(t), 'k:')
    
    # plot samples
    if isinstance(X, mc.Stochastic):
        tr = [X.trace()[:,0], X.trace()[:,1]]
    else:
        tr = [X[0].trace(), X[1].trace()]

    ax1.plot(tr[0], tr[1], 'ko-')
        
    # decorate plot
    book_graphics.set_font()
    pl.xlabel('$X_1$')
    pl.ylabel('$X_2$', rotation=0)
    pl.axis([-1.1,1.1,-1.1,1.1])
    pl.text(-1,1,'(a)', fontsize=16, va='top', ha='left')

    
    for i in range(2):
        if i == 0:
            ax2 = fig.add_subplot(2, 4, 3+4*i)
            ax2.plot(tr[i], 'k', drawstyle='steps-mid')
        else:
            ax2a = fig.add_subplot(2, 4, 3+4*i, sharex=ax2)
            ax2a.plot(tr[i], 'k', drawstyle='steps-mid')
            pl.xlabel('Sample')
        
        pl.xticks([25,50,75])
        pl.yticks([-.5,0,.5])
        pl.ylabel('$X_%d$'%(i+1), rotation=0)
        pl.axis([-5,105,-1.5,1.5])
        pl.text(-1,1.25,'(%s)'%'bc'[i], fontsize=16, va='top', ha='left')
        
        if i == 0:
            ax3 = fig.add_subplot(2, 4, 4+4*i)
            ax3.acorr(tr[i].reshape(100), color='k')
        else:
            ax3a = fig.add_subplot(2, 4, 4+4*i, sharex=ax3)
            ax3a.acorr(tr[i].reshape(100), color='k')
            pl.xlabel('Autocorrelation')
        
        pl.xticks([-5,0,5])
        pl.yticks([0., .5, 1])
        pl.axis([-12,12,-.1,1.1])
        pl.text(-10,1,'(%s)'%'de'[i], fontsize=16, va='top', ha='left')
        
    pl.setp(ax2.get_xticklabels(), visible=False)
    pl.setp(ax3.get_xticklabels(), visible=False)
    pl.subplots_adjust(wspace=.55, hspace=.1, bottom=.14,left=.13)

def CoordsOrbScroll(theta,geo,shaveOn=True, just_involutes = False, Ndict = {}):
    
    from PDSim.scroll import scroll_geo
    
    shaveDelta=None
    if shaveOn==True:
        shaveDelta = pi/2
    else:
        shaveDelta = 1e-16
    
    (xshave, yshave) = Shave(geo, theta, shaveDelta)
    
    Nphi = Ndict.get('phi',500)
    Narc1 = Ndict.get('arc1',100)
    Nline = Ndict.get('line',100)
    Narc2 = Ndict.get('arc2',100)
    
    phi = np.linspace(geo.phi_ois, geo.phi_oie, Nphi)
    (x_oi,y_oi) = scroll_geo.coords_inv(phi,geo,theta,flag="oi")
    phi = np.linspace(geo.phi_oos, geo.phi_ooe - shaveDelta, Nphi)
    (x_oo,y_oo) = scroll_geo.coords_inv(phi,geo,theta,flag="oo")
    
    xarc1=geo.xa_arc1+geo.ra_arc1*cos(np.linspace(geo.t2_arc1,geo.t1_arc1,Narc1))
    yarc1=geo.ya_arc1+geo.ra_arc1*sin(np.linspace(geo.t2_arc1,geo.t1_arc1,Narc1))
    xline=np.linspace(geo.t1_line,geo.t2_line,Nline)
    yline=geo.m_line*xline+geo.b_line
    xarc2=geo.xa_arc2+geo.ra_arc2*cos(np.linspace(geo.t1_arc2,geo.t2_arc2,Narc2))
    yarc2=geo.ya_arc2+geo.ra_arc2*sin(np.linspace(geo.t1_arc2,geo.t2_arc2,Narc2))
    
    ro = geo.rb*(pi-geo.phi_fi0+geo.phi_fo0)
    om = geo.phi_fie-theta+3.0*pi/2.0
    xarc1_o=-xarc1+ro*cos(om)
    yarc1_o=-yarc1+ro*sin(om)
    xline_o=-xline+ro*cos(om)
    yline_o=-yline+ro*sin(om)
    xarc2_o=-xarc2+ro*cos(om)
    yarc2_o=-yarc2+ro*sin(om)
    
    if just_involutes:
        x=np.r_[x_oo,x_oi[::-1]]
        y=np.r_[y_oo,y_oi[::-1]]
    else:
        if shaveOn:
            x=np.r_[x_oo,xshave,x_oi[::-1],xarc1_o,xline_o,xarc2_o]
            y=np.r_[y_oo,yshave,y_oi[::-1],yarc1_o,yline_o,yarc2_o]
        else:
            x=np.r_[x_oo,x_oi[::-1],xarc1_o,xline_o,xarc2_o]
            y=np.r_[y_oo,y_oi[::-1],yarc1_o,yline_o,yarc2_o]
            
    
    #Output as a column vector
    x=x.reshape(len(x),1)
    y=y.reshape(len(y),1)
    return x,y

    def Rot_x(self,x):
        """ Returns rotation matrix in x direction

        Parameters
        ----------
        x: float
           Angle to rotate the matrix with
        """
        return matrix([[1,      0,      0],
                       [0, cos(x), sin(x)],
                       [0,-sin(x), cos(x)]])

    def Rot_y(self,y):
        """ Returns rotation matrix in y direction

        Parameters
        ----------
        y: float
           Angle to rotate the matrix with
        """
        return matrix([[cos(y), 0, -sin(y)],
                       [     0, 1,       0],
                       [sin(y), 0,  cos(y)]])