 xyz=fixed <--- frame' <--- frame'' <--- 123=body
           mphi        mtheta       mpsi

mphi={{Cos[phi],-Sin[phi],0},{Sin[phi],Cos[phi],0},{0,0,1}}
mtheta={{1,0,0},{0,Cos[theta],-Sin[theta]},{0,Sin[theta],Cos[theta]}}
mpsi={{Cos[psi],-Sin[psi],0},{Sin[psi],Cos[psi],0},{0,0,1}}
mphi.mtheta.mpsi
Simplify[%]
Out[5]= {{Cos[phi] Cos[psi] - Cos[theta] Sin[phi] Sin[psi], 
> -(Cos[psi] Cos[theta] Sin[phi]) - Cos[phi] Sin[psi], Sin[phi] Sin[theta]}, 
2> {Cos[psi] Sin[phi] + Cos[phi] Cos[theta] Sin[psi], 
> Cos[phi] Cos[psi] Cos[theta] - Sin[phi] Sin[psi], -(Cos[phi] Sin[theta])}, 
3> {Sin[psi] Sin[theta], Cos[psi] Sin[theta],  Cos[theta]}}

{0,0,dpsi}+Inverse[mpsi].{dtheta,0,0}+Inverse[mpsi].Inverse[mtheta].{0,0,dphi}
Simplify[%]
Out[7]= {dtheta Cos[psi] + dphi Sin[psi] Sin[theta], 
>    -(dtheta Sin[psi]) + dphi Cos[psi] Sin[theta], dpsi + dphi Cos[theta]}
w=%
L={{I1,0,0},{0,I1,0},{0,0,I3}}.w
L.w/2
ke=Collect[%,{I1,I3},Simplify]
                                    2             2       2           2
         I3 (dpsi + dphi Cos[theta])    I1 (dtheta  + dphi  Sin[theta] )
Out[11]= ---------------------------- + --------------------------------
                      2                                2

ppsi=D[ke,dpsi]
Out[13]= I3 (dpsi + dphi Cos[theta])

pphi=D[ke,dphi]
                                                                   2
Out[14]= I3 Cos[theta] (dpsi + dphi Cos[theta]) + dphi I1 Sin[theta]

= ppsi Cos[theta] + dphi I1 Sin[theta]^2

ke2= ppsi^2/(2 I3) + I1 dtheta^2/2 + (pphi-ppsi Cos[theta])^2/(2 I1 Sin[theta]^2)
Simplify[ke2-ke] 
Out[16]= 0

Veff = (Pphi-Ppsi c)^2/(1-c^2)
D[%,c] 
Solve[%==0,c]
                Pphi         Ppsi
Out[19]= {{c -> ----}, {c -> ----}}
                Ppsi         Pphi

Solve[{Pphi==Ppsi Ppsi/Pphi + dphi I1 (1- (Ppsi/Pphi)^2), 
Ppsi== I3 (dpsi + dphi (Ppsi/Pphi))},{dpsi,dphi}]

                   (I1 - I3) Ppsi          Pphi
Out[20]= {{dpsi -> --------------, dphi -> ----}}
                       I1 I3                I1

mphi.mtheta.mpsi.L
Simplify[%]
Out[22]= {dtheta I1 Cos[phi] + (dpsi I3 - dphi (I1 - I3) Cos[theta]) Sin[phi] 
 
>      Sin[theta], dtheta I1 Sin[phi] - 
 
>     Cos[phi] (dpsi I3 - dphi (I1 - I3) Cos[theta]) Sin[theta], 
 
     2 dpsi I3 Cos[theta] + dphi (I1 + I3 + (-I1 + I3) Cos[2 theta])
>    ---------------------------------------------------------------}
                                    2

% /. {dphi->Pphi/I1, dpsi->Cos[theta](I1/I3-1)Pphi/I1,dtheta->0}
Simplify[%]                                                            

Out[24]= {0, 0, Pphi}
conclude: c->Ppsi/Pphi solution has Lfixed directed along z axis

w /. {dphi->Pphi/I1, dpsi->Cos[theta](I1/I3-1)Pphi/I1,dtheta->0}
Simplify[%]                                                            

          Pphi Sin[psi] Sin[theta]  Pphi Cos[psi] Sin[theta]  Pphi Cos[theta]
Out[27]= {------------------------, ------------------------, ---------------}
                     I1                        I1                   I3

mphi.mtheta.mpsi.L
Collect[%[[3]],{I1,I3},Simplify]
                                                                    2
Out[28]= I3 Cos[theta] (dpsi + dphi Cos[theta]) + dphi I1 Sin[theta]

see that this is pphi

Note: Trace of matrix =1 + 2 Cos[angle]
Tr[mphi.mtheta.mpsi]
Simplify[%]
                                phi + psi 2
Out[31]= Cos[phi + psi] + 2 Cos[---------]  Cos[theta]
                                    2
