Thank you Jo?l to reply. The thing is a bit more complex; maybe with the Moon the problem is more clear being already computed.
In the image below you can see the position angle of the bright limb (b.l.), where the Sun's vector is at: (360 - 252.452 = 107.548 E) good agreeding, imho. Good agreeding is also the declination of Sun, which it is important here for the calc.
pos.angle.jpg
BUT, for what concern the position angle of the axis, the things are a bit differents, since it should be visualizable like in the image below, as the basculating of the body's axis about a position which can be "simulated" in the image through the reference frame axis (or the converse, accordingly with the relative point of view drawn from the image). Note that I'm unaware whether such axial toggling is available through SPICE, since I've never experimenticed with them (my fault, indeed), but if one follows, chase, or syncronize the planetary motion never sees such axial displacement, that is: the geographic pole which "slide" about from where is the body axis.
pos.angle_2.jpg
For the Moon I got such motion by adapting parts of the XEphem source (thanks to Fridger to have pointed me there) with SOFA ANSI C source. Through the "celutil.lua", I got the ecliptic coords. of Sol and Moon and did retrive the "vector" latitude:
local BH = math.rad(moon_lat) * math.rad(dist_moon_earth / dist_sol_earth)then is calculated the longitudinal component of the optical libration:
local w = math.rad(sol_long) - omega;note that from the Sol's longitude mighst subtracted the nutation in longitude and the nutation in obliquity, but they are precalculated by Celestia and passed to the "celutil.lua"'s ecliptic coords. function. In fact, if I calculate and do insert them within the equation, the result differ from the check against other ephemerides programs because their re-iteration. "omega" is the argument of the latitude of the Moon.
From there I get the Sun's declination:
local B0 = math.deg(math.asin(-math.sin(w) * math.cos(BH) * math.sin(I) - math.sin(BH) * math.cos(I)));where "I" is the Moon's tilt axis respect to the equator (1.54242 deg.). The axial components are got through:
local pax = math.sin(I) * math.sin(omega);
local pay = math.sin(I) * math.cos(omega) * cosEPS - math.cos(I) * sinEPS;
local paw = math.atan2(pax, pay);where sinEPS, cosEPS are costants: the sine and cosine of the J2000.0 mean obliquity of the ecliptic (IAU 1976):
local sinEPS = 0.3977771559319137;
local cosEPS = 0.9174820620691818;
sinEPS/cosEPS of course can't be avoided since the calcolus is a "standalone" calcolus. Because the axis is toggling, atan2 put in the correct range (+/- < of 25 degrees for the Moon). Finally, its position is:
local pa = math.deg(math.asin(math.sqrt(pax^2 + pay^2) * math.cos(math.rad(moon_long) - paw)
/ math.cos(math.rad(B0))));Such value agreed with the tests at least until 2099, because then the other ephemerides programs which I use does cease to calc. BUT: if I change "omega" and "I" with that of the planets ("omega" from SOFA ANSI C source code, "I" from Wikipedia), I lost the agreeding after 6 months; and note that there isn't nothing in the equations relevant to the Moon only, being irrilevant the optical libration's component "w" once its parameters are given. Thus, I was in search of another method to get the position angle of the axis, at least to check out from my part whether I've done a trivial set-up error with the planets and henceforth to avoid definitely such "omega" fundamental, being, by the way, got without using trigonometric series like Celestia's engine does, and therefore is lesser accurate (I'm somewhat sure to have well-converted arcseconds to degrees, at least...).

EDIT LATER:
"...but if one follows, chase, or syncronize the planetary motion never sees such axial displacement, that is: the geographic pole which "slide" about from where is the body axis." Read: the Moon's axis should moves like the velocity vector, making a circle and the Moon with it.
Never at rest.
Massimo