%% Import Script for EBSD Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = crystalSymmetry('432', [3.52 3.52 3.52], 'mineral', 'Nickel', 'color', 'light blue');

% plotting convention
setMTEXpref('xAxisDirection','north');
setMTEXpref('zAxisDirection','outOfPlane');

%% Define elastic property

pname = 'data_tensor\';
fname = [pname 'IN739LC.GPa'];

C = loadTensor(fname,CS,'propertyname','elastic stiffness','unit','GPa','interface','generic')

figure('position',[0 0 500 400])
plot(C,'plotType','YoungsModulus','3d','complete');
colorbar;
%saveFigure('Nickel_SC_3d.pdf')

figure('position',[0 0 500 400])
plot(C,'plotType','YoungsModulus','contourf',120:20:300,'projection','stereo','complete')
annotate([xvector,-yvector,zvector],'label',{'010','100','001'},'color','k');
colormap('jet');
colorbar;
%saveFigure('Nickel_SC_2d.pdf')

%% Voigt / Reuss / Hill bounds
[C_v,C_r,C_h] = calcTensor(odf,C)
figure('position',[0 0 500 400])
plot(C_v,'plotType','YoungsModulus','contourf',120:20:300,'projection','stereo');
annotate([xvector,-yvector,zvector],'label',{'P','T','B'});
colormap('jet')
%saveFigure('Voigt.pdf')
colorbar;

% isotropic Voigt bounds
cijij = double(EinsteinSum(C_v,[-1 -2 -1 -2]));
ciijj = double(EinsteinSum(C_v,[-1 -1 -2 -2]));
lambda_v = (2*ciijj - cijij)/15;
mue_v = (-ciijj + 3*cijij)/30;
E_v = mue_v *(2*mue_v + 3*lambda_v) / (mue_v + lambda_v)

% Reuss tensor
S_r = inv(C_r)
figure('position',[0 0 500 400])
plot(C_r,'plotType','YoungsModulus','contourf',120:20:300,'projection','stereo');
annotate([xvector,-yvector,zvector],'label',{'P','T','B'});
colormap('jet')
%saveFigure('Reuss.pdf')
colorbar;

% isotropic Reuss bounds
sijij = double(EinsteinSum(S_r,[-1 -2 -1 -2]));
siijj = double(EinsteinSum(S_r,[-1 -1 -2 -2]));
E_r = 15 / (siijj + 2*sijij)
mue_r = 7.5 / (-siijj + 3*sijij);

%% profile BT
ang = 0.0:pi/90.:0.5*pi;
ang0_90 = 0.0:pi/2.:0.5*pi;
vec = vector3d(0.0,-sin(ang),cos(ang));
diag = vector3d(0.7071*sin(ang),0.7071*sin(ang),cos(ang));

figure('position',[0 0 500 400])
hold on
plot(ang,YoungsModulus(C_v,vec),'-b','LineWidth',2)
plot(ang,YoungsModulus(C_r,vec),'-b','LineWidth',2)
plot(ang,YoungsModulus(C,diag),'--r','LineWidth',2)
plot(ang,E_v*abs(vec),'-.k','LineWidth',2)
plot(ang,E_r*abs(vec),'-.k','LineWidth',2)

plot(ang,200,'.k','LineWidth',2)
plot(0,158,'s','MarkerSize',10,'MarkerFaceColor','b','MarkerEdgeColor','k')
plot(pi/2.,237,'s','MarkerSize',10,'MarkerFaceColor','b','MarkerEdgeColor','k')
set(gca,'Ylim',[120 300])
set(gca,'XTick',0.0:pi/12.:pi/2.)
set(gca,'XTickLabel',{'B','15°','30°','45°','60°','75°','T'})
hold off
%saveFigure('BT.pdf')

%% profile TP
figure('position',[0 0 500 400])
hold on
vec = vector3d(-sin(ang),cos(ang),0.);
diag = vector3d(sin(ang+pi/4.),cos(ang+pi/4.),0.0);

plot(ang,YoungsModulus(C_v,vec),'-b','LineWidth',2)
plot(ang,YoungsModulus(C_r,vec),'-b','LineWidth',2)
plot(ang,YoungsModulus(C,diag),'--r','LineWidth',2)
plot(ang,E_v*abs(vec),'-.k','LineWidth',2)
plot(ang,E_r*abs(vec),'-.k','LineWidth',2)

plot(ang,200,'.k','LineWidth',2)
plot(0,237,'s','MarkerSize',10,'MarkerFaceColor','b','MarkerEdgeColor','k')
plot(pi/2.,237,'s','MarkerSize',10,'MarkerFaceColor','b','MarkerEdgeColor','k')

set(gca,'Ylim',[120 300])
set(gca,'XTick',0.0:pi/12.:pi/2.)
set(gca,'XTickLabel',{'T','15°','30°','T+P','60°','75°','P'})
hold off
%saveFigure('TP.pdf')

    
