I like Frank's idea.

You'll have 2 components to the libration: The sun moving north and south of its average position, and the sun moving east and west of its average position.

The north/south movement will only apply to a planet whose axis of rotation differs from its orbital plane.

Your planet's rotation will simply be linear. If the planet is in a circular orbit, its true longitude will also be linear. If you simulate a full orbit and output your data, you can figure out how far the true longitude deviates from linear.

Here's an image I made. I used Gravity Simulator to create a planet around a star with eccentricity = 0.5. I then computed its period and used the screen shot feature to take 12 screen shots every T=1/12 orbital period, with Trails = off. I stacked the images together, along with a complete orbit with Trails = on. Then I made a two-tone circle and pasted a copy of it on each of the 12 positions, rotating it 30 degrees (360/12) each time. This helps illustrate the libration.

In reality, a planet with an eccentricity this high would probably be locked into a spin-coupled resonance like Mercury, but this exaggerated effect helps illustrate the libration.

For the North South libration, give your planet inclination that equals its axial tilt. Then you could output the position vectors over a complete orbit and compute with the z-position how many degrees above or below the ecliptic the planet is. That will tell you how far the star has drifted from its average position in its sky.

Maybe you could use some sine regression software to come up with a formula that fits your data.