Pushing point-spread function reconstruction to the next level: application to SPHERE/ZIMPOL