Context. The study of linear waves and instabilities is necessary to understand the physical evolution of an atmosphere, and can provide physical interpretation of the complex flows found in simulations performed using Global Circulation Models (GCM). In particular, the acceleration of superrotating flow at the equator of hot Jupiters has mostly been studied under several simplifying assumptions, the relaxing of which may impact final results. Aims. We develop and benchmark a publicly available algorithm to identify the eigenmodes of an atmosphere around any initial steady state. We also solve for linear steady states. Methods. We linearise the hydrodynamical equations of a planetary atmosphere in a steady state with arbitrary velocities and thermal profile. We then discretise the linearised equations on an appropriate staggered grid, and solve for eigenvectors and linear steady solutions with the use of a parallel library for linear algebra: ScaLAPACK. We also implement a posteriori calculation of an energy equation in order to obtain more information on the underlying physics of the mode. Results. Our code is benchmarked against classical wave and instability test cases in multiple geometries. The steady linear circulation calculations also reproduce expected results for the atmosphere of hot Jupiters. We finally show the robustness of our energy equation, and its power to obtain physical insight into the modes. Conclusions. We have developed and benchmarked a code for the study of linear processes in planetary atmospheres, with an arbitrary steady state. The calculation of an a posteriori energy equation provides both increased robustness and physical meaning to the obtained eigenmodes. This code can be applied to various problems, and notably to further study the initial spin up of superrotation of GCM simulations of hot Jupiter.