def solver(AA, N_max, 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 = (J_0,J_1,J_2,S...) where
    x contains all Fourier components of the generating function with |n|<N_max """

    # Find all integer component n_vectors which lie within sphere of radius N_max
    # Here we have assumed that the potential is symmetric x->-x, y->-y, z->-z
    # This can be relaxed by changing symN to 1
    # Additionally due to time reversal symmetry S_n = -S_-n so we only consider
    # "half" of the n-vector-space

    angs = unroll_angles(AA.T[3:].T,np.ones(3))

    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

    xxx = check_each_direction(n_vectors,angs)

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

    n = len(n_vectors)+3
    b = np.zeros(shape=(n, ))
    a = np.zeros(shape=(n,n))

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

    for i in AA:
        a[:3,3:]+=2.*n_vectors.T[:3]*np.cos(np.dot(n_vectors,i[3:]))
        a[3:,3:]+=4.*np.dot(n_vectors,n_vectors.T)*np.outer(np.cos(np.dot(n_vectors,i[3:])),np.cos(np.dot(n_vectors,i[3:])))
        b[:3]+=i[:3]
        b[3:]+=2.*np.dot(n_vectors,i[:3])*np.cos(np.dot(n_vectors,i[3:]))

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

    return np.array(solve(a,b)), n_vectors
Constructs the matrix A and the vector b from a timeseries of toy
    action-angles AA to solve for the vector x = (J_0,J_1,J_2,S...) where
    x contains all Fourier components of the generating function with |n|<N_max