For seismic modeling in frequency-domain it is priority to have an accurate solver for wave equation, a.k.a Helmholtz equation, in many values of the frequency. Due to numerical dispersion and "pollution" of the numerical solutions with relatively high frequencies, to design stable solvers for Helmholtz equation it has become a great computational challenge. We study a flexible quadrature scheme with hexagonal grids, which is based on radial basis functions (RBF) interpolation. We benefit of RBF for finding optimal weights for several stencils with hexagonal symmetry. Using these stencils improves the accuracy and mitigates the numerical dispersion. Numerical experiments are presented to show our results. Moreover, we present implementations with perfectly matched layers (PML) and we use this method to compute imaging conditions in single scattering models.