Hfszeeman95 is an updated and extended Fortran 95 version of the Hfszeeman program (Andersson and Jönsson, 2008). Given relativistic atomic state functions generated by the Grasp2018 package (Fischer et al., 2019), Hfszeeman95 together with the accompanying Matlab/GNU Octave program Mithit allows for: (1) the computation and plotting of Zeeman energy splittings of magnetic fine- and hyperfine structure substates as functions of the strength of an external magnetic field, (2) the computation of transition rates between different magnetic fine- and hyperfine structure substates in the presence of an external magnetic field and rates of hyperfine-induced transitions in the field free limit, (3) the synthesization of spectral profiles for transitions obtained from (2). With the new features, Hfszeeman95 and the accompanying Matlab/GNU Octave program Mithit are useful for the analysis of observational spectra and to resolve the complex features due to the splitting of the fine and hyperfine levels. Program summaryProgram Title:Hfszeeman95Program Files doi:http://dx.doi.org/10.17632/rv2vycs7pg.1Licensing provisions: GNU General Public License 3Programming language: Fortran 95, Matlab/GNU OctaveNature of problem: Calculation of transition energies and rates between different magnetic fine- and hyperfine structure substates in the presence of an external magnetic field and rates of hyperfine induced transitions in the field free limit. Synthesization of spectral profiles.Solution method: Wave functions for magnetic fine structure substates in the field free case are given by atomic state functions (ASFs). The ASFs are expansions over configuration state functions (CSFs) |ΓJMJ〉=∑γcγ|γJMJ〉.The ASFs are computed by the Grasp2018 relativistic atomic structure package (Fischer et al., 2018) and are supposed to be available. Wave functions for magnetic fine structure substates in an external magnetic field are expanded in a basis of ASFs |Γ˜MJ〉=∑ΓJdΓJ|ΓJMJ〉.Wave functions for magnetic hyperfine structure substates in an external magnetic field are expanded in a basis of the combined nuclear and atomic system |Γ˜IMF〉=∑ΓJFdΓJF|ΓIJFMF〉,where |ΓIJFMF〉 are coupled nuclear and atomic functions |ΓIJFMF〉=∑MI,MJ〈IJMIMJ|IJFMF〉|IMI〉|ΓJMJ〉.Reduced hyperfine and Zeeman matrix elements, used to construct the total interaction matrix in the given basis, are computed as sums over reduced one-particle matrix elements of orbitals building the CSFs. By diagonalizing the interaction matrix, Zeeman energy splittings of fine- and hyperfine structure substates are obtained together with the expansion coefficients of the basis functions. Transition rates between different magnetic fine- and hyperfine structure substates are computed as sums over reduced transition matrix elements between fine structure states weighted by the expansion coefficients of the basis functions and angular factors. Given energies of the magnetic substates along with transition rates, a synthetic spectrum is obtained by convolving the spectral lines with a Gaussian function with a user defined value of the full width half maximum (FWHM).Additional comments including restrictions and unusual features : 1. The complexity of the cases that can be handled is determined by the Grasp2018 package used for the generation of the atomic state functions. 2. The current programs can be used for the calculations of electric dipole (E1), electric quadrupole (E2), magnetic dipole (M1) and magnetic quadrupole (M2) magnetic-field- and hyperfine-induced transitions, which are caused by the mixing of states due to hyperfine and Zeeman interaction. 3. The present model does not include the necessary non-perturbative treatment of the uncommon case involving near-degeneracies where the radiative width of a fine structure state is of the same order as the hyperfine or magnetic-field perturbation; an effect often termed radiation damping (Indelicato et al., 1989 [1]; Robicheaux et al., 1995 [2]; Johnson et al., 1997 [3]).