%********************************************************************************%
%                                                                                %
%  function tcsolver                                                             %
%                                                                                %
%  function calculates the right hand side of the ODE                            %
%                                                                                %
%********************************************************************************%
%---------------------------------------------------------------------------------
function dy = tcsolver(z,y)

%---------------------------------------------------------------------------------
global gamma theta beta b_const omega;

%---------------------------------------------------------------------------------
%  z is the coordinateje, along which the plane wave propagates
%  y is 4-D vector: y=[vx; vy; tx; ty]
%  F is the right hand side

%---------------------------------------------------------------------------------
% calculation of elastic parameters at point z
%
fi    = b_const*z;

a44   = beta^2*(1+2*gamma*sin(theta)^2*cos(fi)^2);
a45   = 2*beta^2*gamma*sin(theta)^2*sin(fi)*cos(fi);
a55   = beta^2*(1+2*gamma*sin(theta)^2*sin(fi)^2);

A = [ 0   0   1   0;
      0   0   0   1;
     a55 a45  0   0;
     a45 a44  0   0];

B = A^(-1);

right_hand_side = -i*omega*A^(-1)*y;

%--------------------------------------------------------------------------------
% prava strana

dy =	[right_hand_side(1);right_hand_side(2); right_hand_side(3); right_hand_side(4)];