function [ ] = Func_Print_Seismic_Cij_MTEX4(C,rho,elastic_tensor,print_file)
% Function outputs text file with seismic results
%
% Tested with MTEX 4.0.6
%
% fid3 = output text file
% C = MTEX elastic stiffness tensor (GPa)
% rho = density (gm/cm3)
% elastic_tensor = text '100% Olivine : Voigt average'
% print_file =  example '/MatLab_Programs/Seismic_Analysis.txt'
%
% David Mainprice 23/11/2014
%
%%
%***********************************************************************
% open file for output results
%***********************************************************************
%fid3=3;
[fid3] = fopen(print_file,'w');
%***********************************************************************
% Cij Elastic tensor (GPa)
%***********************************************************************
%
% print tensor and density to file
%
fprintf(fid3,'%s \n','Seismic Anisotropy Analysis');
fprintf(fid3,' \n');
fprintf(fid3,'%s %s \n','Print file:',print_file);
fprintf(fid3,' \n');
fprintf(fid3,'%s %s \n','Elastic tensor (GPa):',elastic_tensor);
fprintf(fid3,' \n');
fprintf(fid3,' \n');
fprintf(fid3,' Cij Elastic tensor (GPa) \n');
fprintf(fid3,' \n');
fprintf(' \n')
fprintf('%s %s \n','Elastic tensor (GPa):',elastic_tensor);
fprintf(' \n')
%
% convert tensor to matrix form
MM = matrix(C,'Voigt');
for i=1:6
  fprintf(fid3,'%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',MM(i,:));
  fprintf('%8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n',MM(i,:))
end
%
%***********************************************************************
% Density (g/cm3)
%***********************************************************************
fprintf(' \n')
fprintf(' Density (g/cm3) =  %8.4f \n', rho)
fprintf(' \n')
fprintf(fid3,' \n');
fprintf(fid3,' Density (g/cm3) =  %8.4f \n', rho);
fprintf(fid3,' \n');
%% Seismic anisotropy of an Elastic tensor (C) with density (rho)
%**************************************************************************
% Generate 1 degree pole figure grid on S2
%**************************************************************************
v_grid = plotS2Grid(C,'complete','upper');
%**************************************************************************
% Generate velocities and polarizations on a 1 degree pole figure S2Grid
[vp,vs1,vs2] = velocity(C,v_grid,rho);
%
fprintf(' \n')
fprintf(' Anisotropy Parameters \n')
fprintf(' \n')
fprintf(fid3,' \n');
fprintf(fid3,' Anisotropy Parameters \n');
fprintf(fid3,' \n');
%**************************************************************************
% P-wave Anisotropy (%)
% AVp = 200*(Vp_max_value-Vp_min_value)/(Vp_max_value+Vp_min_value)
%**************************************************************************
Vp_max = max(vp(:));
Vp_min = min(vp(:));
Aggregate_AVp = 200*(Vp_max-Vp_min)/(Vp_max+Vp_min);
% print to screen and file
fprintf(' P-wave anisotropy (%%) =  %6.2f \n', Aggregate_AVp)
fprintf(fid3,' P-wave anisotropy (%%) =  %6.2f \n', Aggregate_AVp);
%
%**************************************************************************
% S-wave anisotropy percentage for each propagation direction
% defined as AVs = 200*(Vs1-Vs2)/(Vs1+Vs2)
%**************************************************************************
% new array AVs %
a_vs = 200*(vs1-vs2)./(vs1+vs2);
Aggregate_AVs_max = max(max(a_vs(:)));
Aggregate_AVs_min = min(min(a_vs(:)));
% print to screen and file
fprintf(' S-wave maximum anisotropy (%%) = %6.2f \n',Aggregate_AVs_max)
fprintf(' S-wave minimum anisotropy (%%) = %6.2f \n',Aggregate_AVs_min)
fprintf(fid3,' S-wave maximum anisotropy (%%) = %6.2f \n',Aggregate_AVs_max);
fprintf(fid3,' S-wave minimum anisotropy (%%) = %6.2f \n',Aggregate_AVs_min);
%**************************************************************************
% dVs max and min for all propagation directions
% defined as dVs = (Vs1-Vs2) (N.B. Anisotropy is same as AVs)
%**************************************************************************
% new array dVs (km/s)
d_vs = vs1-vs2;
Aggregate_dVs_max = max(max(d_vs(:)));
Aggregate_dVs_min = min(min(d_vs(:)));
% print to screen and file
fprintf(' dVs=Vs1-Vs2 maximum (km/s) = %6.2f \n',Aggregate_dVs_max)
fprintf(' dVs=Vs1-Vs2 minimum (km/s) = %6.2f \n',Aggregate_dVs_min)
fprintf(fid3,' dVs=Vs1-Vs2 maximum (km/s) = %6.2f \n',Aggregate_dVs_max);
fprintf(fid3,' dVs=Vs1-Vs2 minimum (km/s) = %6.2f \n',Aggregate_dVs_min);
%**************************************************************************
% S1-wave anisotropy (%)
% AVs1 = 200*(Vs1_max_value-Vs1_min_value)/(Vs1_max_value+Vs1_min_value)
%**************************************************************************
Vs1_max = max(vs1(:));
Vs1_min = min(vs1(:));
% print to screen and file
Aggregate_AVs1 = 200*(Vs1_max-Vs1_min)/(Vs1_max+Vs1_min);
fprintf(' S1-wave anisotropy (%%) = %6.2f \n',Aggregate_AVs1)
fprintf(fid3,' S1-wave anisotropy (%%) = %6.2f \n',Aggregate_AVs1);
%**************************************************************************
% S2-wave anisotropy (%)
% AVs2 = 200*(Vs2_max_value-Vs2_min_value)/(Vs2_max_value+Vs2_min_value)
%**************************************************************************
Vs2_max = max(vs2(:));
Vs2_min = min(vs2(:));
% print to screen and file
Aggregate_AVs2 = 200*(Vs2_max-Vs2_min)/(Vs2_max+Vs2_min);
fprintf(' S2-wave anisotropy (%%) = %6.2f \n',Aggregate_AVs2)
fprintf(fid3,' S2-wave anisotropy (%%) = %6.2f \n',Aggregate_AVs2);
%***********************************************************************
% Vp/Vs1 anisotropy (%) ('vp./vs1')
%***********************************************************************
% new array Vp/Vs1 (km/s)
Vp_Vs1 = vp./vs1;
Vp_Vs1_max = max(Vp_Vs1(:));
Vp_Vs1_min = min(Vp_Vs1(:));
% print to screen and file
Aggregate_AVp_Vs1 = 200*(Vp_Vs1_max-Vp_Vs1_min)/(Vp_Vs1_max+Vp_Vs1_min);
fprintf(' Vp/Vs1 anisotropy (%%) = %6.2f \n',Aggregate_AVp_Vs1)
fprintf(fid3,' Vp/Vs1 anisotropy (%%) = %6.2f \n',Aggregate_AVp_Vs1);
%***********************************************************************
% Vp/Vs1 anisotropy (%) ('vp./vs2')
%***********************************************************************
% new array Vp/Vs2 (km/s)
Vp_Vs2 = vp./vs2;
Vp_Vs2_max = max(Vp_Vs2(:));
Vp_Vs2_min = min(Vp_Vs2(:));
% print to screen and file
Aggregate_AVp_Vs2 = 200*(Vp_Vs2_max-Vp_Vs2_min)/(Vp_Vs2_max+Vp_Vs2_min);
fprintf(' Vp/Vs2 anisotropy (%%) = %6.2f \n',Aggregate_AVp_Vs2)
fprintf(fid3,' Vp/Vs2 anisotropy (%%) = %6.2f \n',Aggregate_AVp_Vs2);
%***********************************************************************
% close export text file
%***********************************************************************
fclose(fid3);
end

