function [Vpx,Vpy,Vpz,Vsvx,Vsvy,Vsvz,Vshx,Vshy,Vshz]=groupTi03(C,fvp1,fvp2,fsv1,fsv2,fsh1,fsh2,sita0,sita,phi,mi)
%% sita0 是vti介质的倾角,sita是入射角,phi方位角,C为vti介质的刚度矩阵,mi是密度
c11=C(1,1);c13=C(1,3);c33=C(3,3);c44=C(4,4);c66=C(6,6);
G=sin(sita)*cos(phi)*cos(sita0)+cos(sita)*sin(sita0);
F=G^2+sin(sita)^2*sin(phi)^2;
E=-sin(sita)*cos(phi)*sin(sita0)+cos(sita)*cos(sita0);
D=((c11-c44)*F-(c33-c44)*E^2)^2+4*(c13+c44)^2*F*E^2;
vp=(((c44+c11)*F+(c33+c44)*E^2+D^(1/2))/(2*mi))^(1/2);
vsv=(((c44+c11)*F+(c33+c44)*E^2-D^(1/2))/(2*mi))^(1/2);
vsh=((c66*F+c44*E^2)/mi)^(1/2);
vp1=vpa(subs(fvp1));
vp2=vpa(subs(fvp2));
sv1=vpa(subs(fsv1));
sv2=vpa(subs(fsv2));
sh1=vpa(subs(fsh1));
sh2=vpa(subs(fsh2));
Vpx=((vp*sin(sita)+cos(sita)*vp1)*cos(phi)-sin(phi)*vp2/sin(sita));
Vpy=((vp*sin(sita)+cos(sita)*vp1)*sin(phi)+cos(phi)*vp2/sin(sita));
Vpz=(vp*cos(sita)-sin(sita)*vp1);
Vsvx=((vsv*sin(sita)+cos(sita)*sv1)*cos(phi)-sin(phi)*sv2/sin(sita));
Vsvy=((vsv*sin(sita)+cos(sita)*sv1)*sin(phi)+cos(phi)*sv2/sin(sita));
Vsvz=(vsv*cos(sita)-sin(sita)*sv1);
Vshx=((vsh*sin(sita)+cos(sita)*sh1)*cos(phi)-sin(phi)*sh2/sin(sita));
Vshy=((vsh*sin(sita)+cos(sita)*sh1)*sin(phi)+cos(phi)*sh2/sin(sita));
Vshz=(vsh*cos(sita)-sin(sita)*sh1);
% Vpx=vpa(subs((vp*sin(sita)+cos(sita)*fvp1)*cos(phi)-sin(phi)*fvp2/sin(sita)));
% Vpy=vpa(subs((vp*sin(sita)+cos(sita)*fvp1)*sin(phi)+cos(phi)*fvp2/sin(sita)));
% Vpz=vpa(subs(vp*cos(sita)-sin(sita)*fvp1));
% Vsvx=vpa(subs((vsv*sin(sita)+cos(sita)*fsv1)*cos(phi)-sin(phi)*fsv2/sin(sita)));
% Vsvy=vpa(subs((vsv*sin(sita)+cos(sita)*fsv1)*sin(phi)+cos(phi)*fsv2/sin(sita)));
% Vsvz=vpa(subs(vsv*cos(sita)-sin(sita)*fsv1));
% Vshx=vpa(subs((vsh*sin(sita)+cos(sita)*fsh1)*cos(phi)-sin(phi)*fsh2/sin(sita)));
% Vshy=vpa(subs((vsh*sin(sita)+cos(sita)*fsh1)*sin(phi)+cos(phi)*fsh2/sin(sita)));
% Vshz=vpa(subs(vsh*cos(sita)-sin(sita)*fsh1));