%***************************************************************************************%
%                                                                                       %
%  TC-FDY                                                                               %
%                                                                                       %
%  calculation of the S-wave coupling in a medium with twisting anisotropy axis         %
%  by the numerical solution of the ordinary differential equation                      %
%  the incident wave is polarized along the y-axis                                      %
%                                                                                       %
%***************************************************************************************%
% description of variables
%
% a44, a45, a55               elastic parameters of twisting transverse isotropy
% gamma                       parameter of of twisting transverse isotropy
% theta                       angle between the symmetry axis and the axis of rotation
% b_const                     rate of rotation of the symmetry axis
% beta                        S-wave velocity along the symmetry axis in TI
% vx,vy                       particle velocities along the x and y coordinates
% tz,ty                       traction along the x and y coordinates
% omega                       frequency
% z                           axis of rotation
% fmin,fmax,df                minimum frequency, maximum frequency, frequency step
%-----------------------------------------------------------------------------------------
global gamma theta beta b_const omega; % global variables used by 'tcsolver'				

%-----------------------------------------------------------------------------------------
% paremeters

beta    = sqrt(6.0);	
gamma   = 0.15;			
theta   = 60*pi/180;
b_const = 0.032;		

gamma_norm = gamma*sin(theta)^2;

z_min = 0;
z_max = pi/b_const;
z_interval = [z_min z_max];

f_min1 = 0.001;   % 1.st frequency band
f_max1 = 0.1;
df1    = 0.00025;

f_min2 = 0.1;     % 2.nd frequency band
f_max2 = 1.0;
df2    = 0.0025;

f_min3 = 1.0;     % 3.rd frequency band
f_max3 = 10.0;
df3    = 0.025;

%-----------------------------------------------------------------------------------------
% output file
out=fopen('tc-fdy.dat','w');

%-----------------------------------------------------------------------------------------
% 1.st loop over frequency
%-----------------------------------------------------------------------------------------
j=0;
for frequency = f_min1:df1:f_max1
    frequency	
    j=j+1;
    omega = 2.*pi*frequency;
%-----------------------------------------------------------------------------------------
% initial conditions

    uy0  = 1.;
    vx0  = 0;
    vy0  = -i*omega*uy0; % incident wave is polarized along the y-axis, vy=duy/dt
    tx0  = 0;
    ty0  = -i*omega*beta*uy0;

%----------------------------------------------------------------------------------------
% initial vector
%
    y0(1) = vx0;
    y0(2) = vy0;
    y0(3) = tx0;
    y0(4) = ty0;

%-----------------------------------------------------------------------------------------
% ODE solver

    options = odeset('AbsTol',1e-14,'RelTol',1e-10, 'Refine', 20);
    [z,y] = ode45('tcsolver',z_interval,y0,options);

    f(j) = frequency;

    ux(j) = y(length(z),1)/(i*omega); 	% calcultaion of the displacement
    uy(j) = y(length(z),2)/(i*omega);

% saving to file
    fprintf(out,'%f %f %f %f %f %f %f\n',f(j),+real(ux(j)),-imag(ux(j)),-real(uy(j)),+imag(uy(j)),abs(ux(j)),abs(uy(j)));
end

%-----------------------------------------------------------------------------------------
% 2.nd loop over frequency
%-----------------------------------------------------------------------------------------
j=0;
for frequency = f_min2:df2:f_max2
    frequency	
    j=j+1;
    omega = 2.*pi*frequency;
%-----------------------------------------------------------------------------------------
% initial conditions

    uy0  = 1.;
    vx0  = 0;
    vy0  = -i*omega*uy0; % incident wave is polarized along the y-axis, vy=duy/dt
    tx0  = 0;
    ty0  = -i*omega*beta*uy0;

%----------------------------------------------------------------------------------------
% initial vector
%
    y0(1) = vx0;
    y0(2) = vy0;
    y0(3) = tx0;
    y0(4) = ty0;

%-----------------------------------------------------------------------------------------
% ODE solver

    options = odeset('AbsTol',1e-10,'RelTol',1e-6, 'Refine', 10);
    [z,y] = ode45('tcsolver',z_interval,y0,options);

    f(j) = frequency;

    ux(j) = y(length(z),1)/(i*omega); 	% calcultaion of the displacement
    uy(j) = y(length(z),2)/(i*omega);

% saving to file
    fprintf(out,'%f %f %f %f %f %f %f\n',f(j),+real(ux(j)),-imag(ux(j)),-real(uy(j)),+imag(uy(j)),abs(ux(j)),abs(uy(j)));
end

%-----------------------------------------------------------------------------------------
% 3.rd loop over frequency
%-----------------------------------------------------------------------------------------
j=0;
for frequency = f_min3:df3:f_max3
    frequency	
    j=j+1;
    omega = 2.*pi*frequency;
%-----------------------------------------------------------------------------------------
% initial conditions

    uy0  = 1.;
    vx0  = 0;
    vy0  = -i*omega*uy0; % incident wave is polarized along the y-axis, vy=duy/dt
    tx0  = 0;
    ty0  = -i*omega*beta*uy0;

%----------------------------------------------------------------------------------------
% initial vector
%
    y0(1) = vx0;
    y0(2) = vy0;
    y0(3) = tx0;
    y0(4) = ty0;

%-----------------------------------------------------------------------------------------
% ODE solver

    options = odeset('AbsTol',1e-8,'RelTol',1e-5, 'Refine', 5);
%    [z,y] = ode45('tcsolver',z_interval,y0,options);

    f(j) = frequency;

    ux(j) = y(length(z),1)/(i*omega); 	% calcultaion of the displacement
    uy(j) = y(length(z),2)/(i*omega);

% saving to file
%    fprintf(out,'%f %f %f %f %f %f %f\n',f(j),+real(ux(j)),-imag(ux(j)),-real(uy(j)),+imag(uy(j)),abs(ux(j)),abs(uy(j)));
end

%-----------------------------------------------------------------------------------------
fclose(out);