This code was presented and used in V. V Makarov & C. T Berghea, PASP, 2026, "Gaia parallax bias via spherical harmonics: A Python tool and discussion of possible causes"(https://ui.adsabs.harvard.edu/abs/2026arXiv260812619M/abstract), where our method is presented in detail and also includes a description of the code. From our abstract: "Parallaxes in Gaia DR3 are known to suffer from a complex set of sky-correlated and magnitude-dependent offsets or biases at the level of a few tens of $\mu$as. Estimated from a sample of one million distant quasars and AGNs from the CRF catalog, the average offset is negative, but the actual distribution of this important parameter shows significant variations on the sky. We propose a practical method to evaluate the parallax correction as a function of sky position and, optionally, of G magnitude using a spherical harmonic series" The main code is varpi3.py. The function which computes the parallax correction is fittingly called "parallax_correction()". It takes as input sky coordinates (ra, dec), fitted SSH coefficients (from the included file sshmean.csv) and optionally Gaia gmag. The input parameters can be single values or Numpy arrays. Here are some examples using the tool: # load the SSH coefficients try: ssh_df = pd.read_csv("sshmean.csv", header=None) ssh = ssh_df.values except FileNotFoundError: print(f"Error: SSH coefficient file not found ") exit() # Example for one coordinate with magnitude correction = parallax_correction(90.0, 75.0, ssh, 11.2) print(f"{correction}") # [-12.384747546510695] # Same coordinates without magnitude correction = parallax_correction(90.0, 75.0, ssh) print(f"{correction}") # [-8.65546462826444] # Example with an array of three coordinates with magnitudes correction = parallax_correction(np.array([220.22, 0. , 359.99, 90. ]), np.array([-89.99, 0. , -45. , 75. ]), ssh, np.array([19.7 , 24. , 19.45, 11.2 ])) print(f"{correction}") # [-14.648153061666195, 0.735334820590475, -13.475237214258028, -12.384747546510695] The code used to produce Figure 4 in our paper is presented as a Jupyter notebook gaiaplx.ipynb. This figure was generated with our code and it shows the sky distribution of Gaia CRF3 parallax bias reconstructed with 81 spherical harmonic functions to degree 8 and is also included here for convenience (SkyPlot.png). Interestingly, the plot shows a curious relation beween the parallax correction and the quasar density dipole reported in recent publications. From our paper abstract: "The directions of the smallest and largest parallax offsets are (l,b) ~ (220 deg,+43 deg) and (l,b) ~ (45 deg, -45 deg), which are close to the orientation of the quasar density dipole reported in recent publications. Motivated by this curious coincidence, we review possible physical effects resulting in a negative bias of measured parallaxes, including an anisotropic universe with a positive curvature and an orbital aberration component" The sky coordinates are loaded from the file testcoo.csv.