Kinematic model for a 5axis mill with universal (nutating) head

That's great to hear. Currently I have a few other things to work on, like my Vismach model. So I'm all good to wait until you've come to a resolution you're happy with.

If I was putting a control on the machine you pictured I would set up the post processor to put out A and B axis commands (tool tilt) and let the post processor convert those position commands to rotary axis positions.
That may be what your plan is but just thought I would add my thoughts as early as possible.
I would also implement electronic gearing. This would minimize hesitation when moving an axis that needs the second axis to move to keep the tool oriented properly. I saw this on a 5 axis BC head with a Fanuc 15MB control. When electronic gearing was turned on the BC motion was noticeably crisper and smoother.
You can also implement 5 axis tool length offsets if the kinematics is in the machine control instead of the post processor. The 5 axis Cinci (with LinuxCNC) I had was able to use post processed NC programs using any tool length/pivot length combination by putting adjustment values in the tool length register. The important value in 5 axis positioning is the pivot point of the rotary axes. The control needs to position the XYZ of the pivot point using the rotation values of the tool angle(s). By implementing tool length and pivot length into the machine control kinematics you can also have TCP (Tool Control Point) motion.
When compensating a rotary axis that can rotate several rotations you also need to make sure your position compensation adds to ZERO for one complete rotation or you will have compensation addition or subtraction error for rotations more than one complete axis rotation.

I have just updated the config folder in the opening post.

Implemented:
G68.2                  : defines twp using euler-angles, pitch-roll-yaw, 2-vectors, 3 points, optionally with offset in XYZ and rotation in XY
G68.3                  : defines twp from current tool orientation, optionally with offset in XYZ and rotation in XY
G68.4                  : same as G68.2 but as an incremental definition from an active TWP plane
G69                      : cancels the current twp (resets all parameters, moves to G54 and sets Identity kinematics)
G53.1 (P)              : spindle orientation without tcp, switches to G59 and activates tool kinematics
G53.3 (P XYZ)      : same as G53.1 but  with simultaneous move the the XYZ coords on the twp plane 
G53.6 (P)             : same as G53.1 but spindle orientation with tcp

- Spindle is  C primary, A secondary or B secondary as defined in the [TWP] section of the ini file
- All G53.x commands will respect axis limits as set in the ini file for the respective primary and secondary spindle joints.
- The P word sets the orientation strategy: 0(default)=shortest distance, 1=positive rotation only, 2= negative rotation only (this applies to the primary rotary, the secondary move the shortest distance)
 

Some more details for which I did not have time yesterday:

- The tool kinematic mode now operates with the spindle joint angles calculated by the remap instead of the actual joint position. The operator can now jog the spindle rotation without rotating the whole coordinate system as described in the documentation which makes it a lot more fool proof.

- All angles (ie primary, secondary and virtual pre-rotation) calculated to orient the tool to the requested TWP are run through the respective tool-kinematic model in the remap to verify that the values are indeed correct.

Sir, why does this py program take up so many cpu processes? Is there any way to optimize it?

I'm afraid, I cannot help you there but I'm sure the matrix multiplications and atan2 functions are computationally not the most lightweight.