Eq0l:E1*cos(d1)*cos(p)+E2*cos(d2)*sin(p);

Eq0r:Ea*cos(d);

Eq1l:-E1*cos(d1)*sin(p)+E2*cos(d2)*cos(p);

Eq1r:Eb*sin(d);

Eq2l:E1*sin(d1)*cos(p)+E2*sin(d2)*sin(p);

Eq2r:Ea*sin(d);

Eq3l:E1*sin(d1)*sin(p)-E2*sin(d2)*cos(p);

Eq3r:Eb*cos(d);

Sol1l:trigsimp(ratsimp(Eq0l*Eq3l-Eq1l*Eq3l));

Sol1r:trigsimp(ratsimp(Eq0r*Eq3r-Eq1r*Eq3r));

Sol1:trigsimp(ratsimp(Sol1l-Sol1r));

Sol2l:trigsimp(ratsimp(Eq0l/Eq2l+Eq1l/Eq3l));

Sol2r:trigsimp(ratsimp(Eq0r/Eq2r+Eq1r/Eq3r))

Sol2:trigsimp(ratsimp(Sol2l-Sol2r))

Q(S,chi,psi)=Q:S*cos(2*chi)*cos(2*psi);

U(S,chi,psi)=U:S*cos(2*chi)*sin(2*psi);

V(S,chi,psi)=V:S*sin(2*chi);

Jp:trigsimp(jacobian([Q,U,V], [S,chi,psi]));

hQ:Jp.[1,0,0];