def angle_solver(AA, timeseries, N_max, sign, symNx = 2, throw_out_modes=False):
    """ Constructs the matrix A and the vector b from a timeseries of toy
    action-angles AA to solve for the vector x = (theta_0,theta_1,theta_2,omega_1,
    omega_2,omega_3, dSdx..., dSdy..., dSdz...) where x contains all derivatives
    of the Fourier components of the generating function with |n| < N_max """

    # First unroll angles
    angs = unroll_angles(AA.T[3:].T,sign)

    # Same considerations as above
    symNz = 2
    NNx = range(-N_max, N_max+1, symNx)
    NNy = range(-N_max, N_max+1, symNz)
    NNz = range(-N_max, N_max+1, symNz)
    n_vectors = np.array([[i,j,k] for (i,j,k) in product(NNx,NNy,NNz)
                          if(not(i==0 and j==0 and k==0)    # exclude zero vector
                             and (k>0                          # northern hemisphere
                                  or (k==0 and j>0)                 # half of x-y plane
                                  or (k==0 and j==0 and i>0))       # half of x axis
                             and np.sqrt(i*i+j*j+k*k)<=N_max     # inside sphere
                             )])

    if(throw_out_modes):
        n_vectors = np.delete(n_vectors,check_each_direction(n_vectors,angs),axis=0)

    nv = len(n_vectors)
    n = 3*nv+6

    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

    a[:3,:3]=len(AA)*np.identity(3)
    a[:3,3:6]=np.sum(timeseries)*np.identity(3)
    a[3:6,:3]=a[:3,3:6]
    a[3:6,3:6]=np.sum(timeseries*timeseries)*np.identity(3)

    for i,j in zip(angs,timeseries):
        a[6:6+nv,0]+=-2.*np.sin(np.dot(n_vectors,i))
        a[6:6+nv,3]+=-2.*j*np.sin(np.dot(n_vectors,i))
        a[6:6+nv,6:6+nv]+=4.*np.outer(np.sin(np.dot(n_vectors,i)),np.sin(np.dot(n_vectors,i)))

        b[:3]+=i
        b[3:6]+=j*i

        b[6:6+nv]+=-2.*i[0]*np.sin(np.dot(n_vectors,i))
        b[6+nv:6+2*nv]+=-2.*i[1]*np.sin(np.dot(n_vectors,i))
        b[6+2*nv:6+3*nv]+=-2.*i[2]*np.sin(np.dot(n_vectors,i))

    a[6+nv:6+2*nv,1]=a[6:6+nv,0]
    a[6+2*nv:6+3*nv,2]=a[6:6+nv,0]
    a[6+nv:6+2*nv,4]=a[6:6+nv,3]
    a[6+2*nv:6+3*nv,5]=a[6:6+nv,3]
    a[6+nv:6+2*nv,6+nv:6+2*nv]=a[6:6+nv,6:6+nv]
    a[6+2*nv:6+3*nv,6+2*nv:6+3*nv]=a[6:6+nv,6:6+nv]

    a[:6,:]=a[:,:6].T

    return np.array(solve(a,b))
Constructs the matrix A and the vector b from a timeseries of toy
    action-angles AA to solve for the vector x = (theta_0,theta_1,theta_2,omega_1,
    omega_2,omega_3, dSdx..., dSdy..., dSdz...) where x contains all derivatives
    of the Fourier components of the generating function with |n| < N_max