%
% GEH calculation code V.1.1
% Erick Andres Perez Alday, PhD, <perezald@ohsu.edu>
% Annabel Li-Pershing, BS, < lipershi@ohsu.edu>
% Muammar Kabir, PhD, < muammar.kabir@gmail.com >
% Larisa Tereshchenko, MD, PhD, < tereshch@ohsu.edu >
% Last update: February 20th, 2018

clear
clc
close all

warning('OFF');

% Input:
%      1. sampling frequency in HZ
%      2. matlab file containing the following:
%        - XYZ median beat
%        - index of QRS onset Vector Magnitudes
%        - index of R peak on Vector Magnitudes
%        - index of QRS offset Vector Magnitudes
%        - index of T peak on Vector Magnitudes
%        - index of T offset on Vector Magnitudes
%        - index of origin point

% Output:
%      1. Excel inserted with calculated GEH parameters
%      2. Mat file containing the calculated GEH parameters
%      3. 3D fig file of GEH plot
%      4. 2D jpeg of AUC VM plot




%% =========================== File and user input =============================
% sampling frequency

% sampling frequency
prompt='What is the sampling frequency in Hz';
dlg_title='Sample Frequency';
fs_d = inputdlg(prompt,dlg_title,1);


% Amplitude resolution
prompt='What is the amplitude resolution in MicroVolts';
dlg_title='Amplitude';
amp_r = inputdlg(prompt,dlg_title,1);


fs=str2num(fs_d{1});
amp_r=str2num(amp_r{1});

% load mat files
[file_name,path_name] = uigetfile('*','Select file for GEH Analysis');
path_name_saving=path_name;

%% file import
file_ID          = strsplit(file_name,'.');
file_path        = fullfile(path_name,file_name);


matfile          = matfile(file_path)

%% ===================== load variables from .mat file ========================
XYZ_median      = matfile.XYZ_M*amp_r;
R_VM            = matfile.R_VM;
q_points_VM     = matfile.q_points_VM;
s_points_VM     = matfile.s_points_VM;
tp_points_VM    = matfile.tp_points_VM;
te_points_VM    = matfile.te_points_VM;
OriginPoint_idx = matfile.Opnt_Me(1,1);

%% ============================  initiate excel file ===========================
excel_file_name = 'All_Results.xls';
Results_folder = [path_name_saving, '/', 'Results', '/'];
if (exist(Results_folder,'file') ~= 7)
    mkdir (Results_folder);
end
% write title of the variables
if ~exist(strcat(Results_folder,excel_file_name))
    fid1 = fopen(strcat(Results_folder,excel_file_name),'a');
    fprintf(fid1,'ID \t');
    fprintf(fid1,'peak QRST Angle_deg \t');
    fprintf(fid1,'area QRST Angle_deg \t');
    fprintf(fid1,'peak QRS Azimuth_deg \t');
    fprintf(fid1,'area QRS Azimuth_deg \t');
    fprintf(fid1,'peak T Azimuth_deg \t');
    fprintf(fid1,'area T Azimuth_deg \t');
    fprintf(fid1,'peak SVG Azimuth_deg \t');
    fprintf(fid1,'area SVG Azimuth_deg \t');
    fprintf(fid1,'peak QRS Elevation_deg \t');
    fprintf(fid1,'area QRS Elevation_deg \t');
    fprintf(fid1,'peak T Elevation_deg \t');
    fprintf(fid1,'area T Elevation_deg \t');
    fprintf(fid1,'peak SVG Elevation_deg \t');
    fprintf(fid1,'area SVG Elevation_deg \t');
    fprintf(fid1,'peak QRS Magnitude_uV \t');
    fprintf(fid1,'area QRS_uVms \t');
    fprintf(fid1,'peak T Magnitude_uV \t');
    fprintf(fid1,'area T_uVms \t');
    fprintf(fid1,'peak SVG Magnitude_uV \t');
    fprintf(fid1,'QT Interval_ms \t');
    fprintf(fid1,'AUC of QTVectorMagnitude_uVms \t');
    fprintf(fid1,'Wilson SVG_uVms \t');
    fprintf(fid1,'\n \n');
    fclose(fid1);
end


%% =============================================================================

% =========================== Variable Calculation ============-================
% calculation of Vector Magnitudes (Euclidian norm) using the origin point

for ii=1:length(XYZ_median (:,1))
	VecMag(ii)=norm([XYZ_median(OriginPoint_idx,1)-XYZ_median(ii,1),XYZ_median(OriginPoint_idx,2)-XYZ_median(ii,2),XYZ_median(OriginPoint_idx,3)-XYZ_median(ii,3)]);
end



% find R peak in median XYZ beat
[Rx_val, Rx]  = max(XYZ_median(1:500,1));
[Ry_val, Ry]  = max(XYZ_median(1:500,2));
[Rz_val, Rz]  = max(XYZ_median(1:500,3));

% define R peak and T peak as R axis and T axis
Raxis = XYZ_median(R_VM,:);
Taxis = XYZ_median(tp_points_VM(1),:);

% ========================== Calculate AUC on Vector Magnitude ===========================
spac_incr = 1000/fs; % spacing increment for trapz calculation
AUC_VM_QT=0;
AUC_VM_QT = trapz(abs(VecMag(q_points_VM(1,1):te_points_VM(1,1)))) * spac_incr;



% ======================= GEH Variable Calculation =============================
%  origin point
CP=XYZ_median(OriginPoint_idx,:);
% Y axis vector
Ynew=[0,1,0];



% QRS and T integration: for Wilson SVG calculation
SumVGx=trapz(XYZ_median(q_points_VM(1,1):te_points_VM(1,1),1))*spac_incr;
SumVGy=trapz(XYZ_median(q_points_VM(1,1):te_points_VM(1,1),2))*spac_incr;
SumVGz=trapz(XYZ_median(q_points_VM(1,1):te_points_VM(1,1),3))*spac_incr;


% QRS and T integration for area vectors
meanVxQ=trapz(XYZ_median(q_points_VM(1,1):s_points_VM(1,1),1))*spac_incr;
meanVxT=trapz(XYZ_median(s_points_VM(1,1):te_points_VM(1,1),1))*spac_incr;
meanVyQ=trapz(XYZ_median(q_points_VM(1,1):s_points_VM(1,1),2))*spac_incr;
meanVyT=trapz(XYZ_median(s_points_VM(1,1):te_points_VM(1,1),2))*spac_incr;
meanVzQ=trapz(XYZ_median(q_points_VM(1,1):s_points_VM(1,1),3))*spac_incr;
meanVzT=trapz(XYZ_median(s_points_VM(1,1):te_points_VM(1,1),3))*spac_incr;



% QRS area and T area vectors based on integrals
MEAN_QRSO=[meanVxQ meanVyQ meanVzQ];
MEAN_TO=[meanVxT meanVyT meanVzT];




%% QT interval
timeM=((1:length(VecMag))/fs)*1000;
QT_interval=timeM(te_points_VM(1,1))-timeM(q_points_VM(1,1));


% peak vectors QRS and T amplitude
QRS_amp=sqrt(Raxis(1)^2+Raxis(2)^2+Raxis(3)^2);
T_amp=sqrt(Taxis(1)^2+Taxis(2)^2+Taxis(3)^2);

% peak SVG vector and mean SVG vector calculation as vector sum of QRS and T vectors 

SVG_axis=sum([Taxis ; Raxis]);
SVG_MO=sum([MEAN_TO; MEAN_QRSO]);




% Origin Point-P, Q-S and S-T vector calculation
qs3 = XYZ_median(q_points_VM(1,1):s_points_VM(1,1),3);
qs1 = XYZ_median(q_points_VM(1,1):s_points_VM(1,1),1);
qs2 = XYZ_median(q_points_VM(1,1):s_points_VM(1,1),2);
st3 = XYZ_median(s_points_VM(1,1):te_points_VM(1,1),3);
st1 = XYZ_median(s_points_VM(1,1):te_points_VM(1,1),1);
st2 = XYZ_median(s_points_VM(1,1):te_points_VM(1,1),2);



% ============================= Angle Claculation ==============================

% peak QRS-T angle
QRSTang=rad2deg(acos(dot(Raxis,Taxis)/(sqrt(Raxis(1)^2+Raxis(2)^2+Raxis(3)^2)*sqrt(Taxis(1)^2+Taxis(2)^2+Taxis(3)^2))));

% mean QRS-T angle
QRSTang_M=rad2deg(acos(dot(MEAN_QRSO,MEAN_TO)/(sqrt(MEAN_QRSO(1)^2+MEAN_QRSO(2)^2+MEAN_QRSO(3)^2)*sqrt(MEAN_TO(1)^2+MEAN_TO(2)^2+MEAN_TO(3)^2))));

% Azimuth of QRS: peak, area 
AZ_OQ=(rad2deg(acos(Raxis(1)/sqrt(Raxis(1)^2+Raxis(3)^2))))*(((Raxis(3)<0)*-1)+((Raxis(3)>0)*1));
AZ_OQM=(rad2deg(acos( MEAN_QRSO(1)/sqrt(MEAN_QRSO(1)^2+MEAN_QRSO(3)^2))))*(((MEAN_QRSO(3)<0)*-1)+((MEAN_QRSO(3)>0)*1));

% Azimuth of T: peak, area 
AZ_OT=(rad2deg(acos(Taxis(1)/sqrt(Taxis(1)^2+Taxis(3)^2))))*(((Taxis(3)<0)*-1)+((Taxis(3)>0)*1));
AZ_OTM=(rad2deg(acos(MEAN_TO(1)/sqrt(MEAN_TO(1)^2+MEAN_TO(3)^2))))*(((MEAN_TO(3)<0)*-1)+((MEAN_TO(3)>0)*1));

% Azimuth of SVG: peak, area 
AZ_SVG=(rad2deg(acos(SVG_axis(1)/sqrt(SVG_axis(1)^2+SVG_axis(3)^2))))*(((SVG_axis(3)<0)*-1)+((SVG_axis(3)>0)*1));
AZ_SVG_M=(rad2deg(acos(SVG_MO(1)/sqrt(SVG_MO(1)^2+SVG_MO(3)^2))))*(((SVG_MO(3)<0)*-1)+((SVG_MO(3)>0)*1));

% Elevation of QRS: peak, area 
EL_OQ=(rad2deg(acos(dot(Raxis,Ynew)/(sqrt(Raxis(1)^2+Raxis(2)^2+Raxis(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));
EL_OQM=(rad2deg(acos(dot(MEAN_QRSO,Ynew)/(sqrt(MEAN_QRSO(1)^2+MEAN_QRSO(2)^2+MEAN_QRSO(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));

% Elevation of T: peak, area 
EL_OT=(rad2deg(acos(dot(Taxis,Ynew)/(sqrt(Taxis(1)^2+Taxis(2)^2+Taxis(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));
EL_OTM=(rad2deg(acos(dot(MEAN_TO,Ynew)/(sqrt(MEAN_TO(1)^2+MEAN_TO(2)^2+MEAN_TO(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));

% Elevation of SVG: peak, area 
EL_SVG=(rad2deg(acos(dot(SVG_axis,Ynew)/(sqrt(SVG_axis(1)^2+SVG_axis(2)^2+SVG_axis(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));
EL_SVG_M=(rad2deg(acos(dot(SVG_MO,Ynew)/(sqrt(SVG_MO(1)^2+SVG_MO(2)^2+SVG_MO(3)^2)*sqrt(Ynew(1)^2+Ynew(2)^2+Ynew(3)^2)))));


% ========================== Magnitudes Calculation ============================
% Magnitude of QRS: peak, area
QRS_Mag=sqrt(Raxis(1)^2+Raxis(2)^2+Raxis(3)^2);
QRS_Mag_M=sqrt(MEAN_QRSO(1)^2+MEAN_QRSO(2)^2+MEAN_QRSO(3)^2);

% Magnitude of T: peak, area
T_Mag=sqrt(Taxis(1)^2+Taxis(2)^2+Taxis(3)^2);
T_Mag_M=sqrt(MEAN_TO(1)^2+MEAN_TO(2)^2+MEAN_TO(3)^2);

% Magnitude of SVG: peak
SVG_Mag=sqrt(SVG_axis(1)^2+SVG_axis(2)^2+SVG_axis(3)^2);

% Magnitude of WVG: Wilson's Ventricular Gradient

WVG=sqrt((SumVGx^2) + (SumVGy^2) + (SumVGz^2));



%% =============================================================================
% ================================== Plotting ==================================





%% ============================ plot AUC on Vector Magnitude =============================



bLine_VM=ones(length(VecMag),1)*VecMag(OriginPoint_idx);
FT=[timeM(q_points_VM(1,1):te_points_VM(1,1)),fliplr(timeM(q_points_VM(1,1):te_points_VM(1,1)))];
F_VM=[VecMag(q_points_VM(1,1):te_points_VM(1,1))';fliplr(bLine_VM(q_points_VM(1,1):te_points_VM(1,1)))];



VM_AUC=figure('Name','AUC on VecMag','visible','on');

subplot(211)
v1=plot(timeM,VecMag,'k','DisplayName','Vector Magnitude');
hold on
fill(FT,F_VM,[0 0 0]+0.75)
hold off
xlabel('Time (ms)','FontSize',10);
ylabel('Amplitude (mV)','FontSize',10);
title('AUC of VectorMagnitude')
subplot(212)
x1=plot(timeM,XYZ_median(:,1),'DisplayName','X Lead');
hold on
x2=plot(timeM,XYZ_median(:,2),'DisplayName','Y Lead');
x3=plot(timeM,XYZ_median(:,3),'DisplayName','Z Lead');
xlabel('Time (ms)','FontSize',10);
ylabel('Amplitude (mV)','FontSize',10);

lgd = legend([v1 x1 x2 x3]);
hold off

%% =============================== GEH Plot ====================================

map_c=hsv;
map_c1=map_c(1:end-8,:);


Fig3D=figure('visible','on','outerposition',[0 0 1400 1000]);
ax3D = axes('Parent',Fig3D);
%% plotting GEH vectors and loops
hold on
p1 = plot3([CP(:,3) Raxis(3)],[CP(:,1) Raxis(1)],[CP(:,2) Raxis(2)],'r','LineWidth',2, 'DisplayName', 'Peak QRS');
p2 = plot3([CP(:,3) Taxis(3)],[CP(:,1) Taxis(1)],[CP(:,2) Taxis(2)],'g','LineWidth',2, 'DisplayName', 'Peak T');
p7 = plot3(CP(:,3),CP(:,1),CP(:,2),'m+','LineWidth',3,'DisplayName', 'Ori. point');
p5 = plot3([CP(:,3) SVG_axis(3)],[CP(:,1) SVG_axis(1)],[CP(:,2) SVG_axis(2)],'b','LineWidth',2, 'DisplayName', 'Peak SVG');


% for plotting in patch function, need to decimate the last data point in the following 2 leads
qs1(end) = nan; qs2(end) = nan;
qrs_loop=patch(qs3,qs1,qs2,(1:spac_incr:spac_incr*length(qs1)),'EdgeColor','interp','DisplayName','QRS Loop');

% plot QRS loop point with the following
for c_ind = 1:5:length(qs1)
scatter3(qs3(c_ind),qs1(c_ind),qs2(c_ind),20,spac_incr*c_ind,'filled' );
end

% for plotting in patch function, need to decimate the last data point in the following 2 leads
st1(end) = nan; st2(end) = nan;
st_loop=patch(st3,st1,st2,(spac_incr*(length(qs1)+1):spac_incr:spac_incr*(length(qs1)+length(st1))),'EdgeColor','interp','DisplayName','T Loop');


% plot T loop with the following
%st_color = [(0.1:0.9/(length(st3)-1):1)',zeros(length(st3),1),flipud((0.1:0.9/(length(st3)-1):1)')];


for c_ind = 1:5:length(st3)
scatter3(st3(c_ind),st1(c_ind),st2(c_ind),20,spac_incr*(length(qs1)+c_ind),'filled');
end

% plot the projected lines to show SVG is the sum of QRS and T vectors
plot3([Raxis(3) SVG_axis(3)], [Raxis(1) SVG_axis(1)],[Raxis(2) SVG_axis(2)],'--','Color', [0.8 0.8 0.8],'LineWidth',2);
plot3([Taxis(3) SVG_axis(3)], [Taxis(1) SVG_axis(1)],[Taxis(2) SVG_axis(2)],'--','Color', [0.8 0.8 0.8],'LineWidth',2);

map_c=hsv;
map_c1=map_c(1:end-8,:);
caxis([0 2*(length(qs1)+length(st1))]);
colormap(map_c1)
colorbar



set(gca,'YDir','reverse')
set(gca,'ZDir','reverse')


%% =========================== Additional Graphic ==============================

% ================= smiley face ==================
% location reference at the center of smiley face and coordinate
x_face = 400;
y_face = -400;
z_face = -400;
r_face = 150;
r_smiley = 110;
THK = 10; % thickness of face
d_eye = 75; % distance from the center of circle
h_eye = -50; % height from the center of the circle
x_coor = 400;
y_coor = -400;
z_coor = -200;
l_arrow = 150;
r_ctr = 10;

% change x y z above to move the smiley face
[x, y, z] = cylinder(150);
surf(z*THK+x_face, x+y_face,y+z_face, 'EdgeColor','none','FaceColor','y');
theta = -pi:0.05:pi;
x_circle = r_face*cos(theta);
y_circle = r_face*sin(theta);
% ================== front face ====================
fill3(zeros(1,length(x_circle))+x_face,x_circle+y_face,y_circle+z_face,'y')
% ================== back face =====================
fill3(zeros(1,length(x_circle))+x_face+THK,x_circle+y_face,y_circle+z_face,[0.5 0.5 0.5])
% =================== make eyes ====================
[x_sphere, y_sphere, z_sphere] = sphere; % unit sphere
 surf(x_sphere*5+x_face, y_sphere*THK+y_face+d_eye,z_sphere*THK+z_face+h_eye, 'EdgeColor', 'none', 'FaceColor', 'k')
 surf(x_sphere*5+x_face, y_sphere*THK+y_face-d_eye,z_sphere*THK+z_face+h_eye, 'EdgeColor', 'none', 'FaceColor', 'k')

% ==================== make nose ====================
t_cone = 0:0.005:0.5;
cone_clr = [255 171 0]./255;
cone_len = -40;
[x_cone, y_cone,z_cone] = cylinder(t_cone);
surf(-z_cone*cone_len+x_face+cone_len,x_cone*2*cone_len+y_face,y_cone*2*cone_len+z_face,'EdgeColor', 'none', 'FaceColor', [255 171 0]./255)

% ==================== make smile ====================
theta_smile = pi*0.2:0.05:pi*0.8;
x_smiley = r_smiley*cos(theta_smile);
y_smiley = r_smiley*sin(theta_smile);
 plot3(ones(1,length(x_smiley))+x_face-THK, x_smiley+y_face, y_smiley+z_face+10, 'k', 'LineWidth', 6)

% ===================== coordinate ====================
% coordinate center
surf(x_sphere*r_ctr+x_coor, y_sphere*r_ctr+y_coor, z_sphere*r_ctr+z_coor, 'EdgeColor', 'none', 'FaceColor', 'k');
% arrows
l = mArrow3([x_coor,y_coor,z_coor], [x_coor,y_coor+l_arrow,z_coor], 'color', 'r', 'stemWidth', 3);
text(x_coor,y_coor+l_arrow,z_coor, '\bf Left', 'FontSize', 6, 'HorizontalAlignment','right');
f = mArrow3([x_coor,y_coor,z_coor], [x_coor+l_arrow,y_coor,z_coor], 'color', 'b', 'stemWidth', 3);
text(x_coor+l_arrow+10,y_coor,z_coor, '\bf Posterior', 'FontSize', 6, 'HorizontalAlignment','left');
d = mArrow3([x_coor,y_coor,z_coor], [x_coor,y_coor,z_coor+l_arrow+20], 'color', 'g', 'stemWidth', 3);
text(x_coor,y_coor,z_coor+l_arrow+15, '\bf Inferior', 'FontSize', 6, 'HorizontalAlignment','center');

% ============================== plot properties ===============================
axis equal
grid on
box on
xlabel('Z (mV)','FontSize',20);
ylabel('X (mV)','FontSize',20);
zlabel('Y (mV)','FontSize',20);
zTickLabel = ax3D.XTick;
zTick = cellfun(@str2num, ax3D.XTickLabel);
set(gca,'XTick',zTick );
set(gca,'XTickLabel',zTickLabel);
view(290,60);

% plot legend and file ID
hold off


lgd = legend([p1 p2 p5 p7 qrs_loop st_loop],'location','bestoutside');
title(lgd,file_ID{1});

%% =============================================================================

%% ==================== Write Calculated Variables to Excel ====================
fid1 = fopen(strcat(Results_folder,excel_file_name),'a');
fprintf(fid1,'%s \t',file_ID{1});
fprintf(fid1,'%f \t',QRSTang);
fprintf(fid1,'%f \t',QRSTang_M);
fprintf(fid1,'%f \t',AZ_OQ);
fprintf(fid1,'%f \t',AZ_OQM);
fprintf(fid1,'%f \t',AZ_OT);
fprintf(fid1,'%f \t',AZ_OTM);
fprintf(fid1,'%f \t',AZ_SVG);
fprintf(fid1,'%f \t',AZ_SVG_M);
fprintf(fid1,'%f \t',EL_OQ);
fprintf(fid1,'%f \t',EL_OQM);
fprintf(fid1,'%f \t',EL_OT);
fprintf(fid1,'%f \t',EL_OTM);
fprintf(fid1,'%f \t',EL_SVG);
fprintf(fid1,'%f \t',EL_SVG_M);
fprintf(fid1,'%f \t',QRS_Mag);
fprintf(fid1,'%f \t',QRS_Mag_M);
fprintf(fid1,'%f \t',T_Mag);
fprintf(fid1,'%f \t',T_Mag_M);
fprintf(fid1,'%f \t',SVG_Mag);
fprintf(fid1,'%f \t',QT_interval);
fprintf(fid1,'%f \t',AUC_VM_QT);
fprintf(fid1,'%f \t',WVG);
fprintf(fid1,'\n');
fclose(fid1);

%% =============================================================================

%% ======================== Save Plot to Image folder ==========================
Images_folder = [Results_folder 'All Figures' '/'];
if (exist(Images_folder) == 0)
  mkdir (Images_folder);
end

filepath = [Images_folder  strcat(file_ID{1},'_VM_AUC_QT.jpg')];
saveas(VM_AUC, filepath, 'jpg');

filepath = [Images_folder  strcat(file_ID{1},'_3D.fig')];
saveas(Fig3D, filepath, 'fig');
filepath = [Images_folder  strcat(file_ID{1},'_3D.jpg')];
saveas(Fig3D, filepath, 'jpg');


%% =============================================================================

%% ================== Write Calculated Variables to .mat file ==================
Mat_folder = [Results_folder 'Mat Files' '/'];
if (exist(Mat_folder) == 0)
mkdir (Mat_folder);
end

% Save to a new MAT file
save(strcat(Mat_folder,'/',file_ID{1},'.mat'), 'CP','AZ_OQ', 'AZ_OQM', 'AZ_OT', 'AZ_OTM', 'AZ_SVG', 'AZ_SVG_M', ...
    'EL_OQ', 'EL_OQM','EL_OT', 'EL_OTM', 'EL_SVG', 'EL_SVG_M', 'MEAN_QRSO', 'MEAN_TO', ...
    'QRSTang', 'QRSTang_M', 'T_Mag', 'T_Mag_M', 'QRS_Mag', 'QRS_Mag_M', 'SVG_axis',  'SVG_Mag', 'Taxis', 'AUC_VM_QT', 'WVG');

%% =============================================================================

%% ============================ Additional function ============================
    function h = mArrow3(p1,p2,varargin)
    %mArrow3 - plot a 3D arrow as patch object (cylinder+cone)
    %
    % syntax:   h = mArrow3(p1,p2)
    %           h = mArrow3(p1,p2,'propertyName',propertyValue,...)
    %
    % with:     p1:         starting point
    %           p2:         end point
    %           properties: 'color':      color according to MATLAB specification
    %                                     (see MATLAB help item 'ColorSpec')
    %                       'stemWidth':  width of the line
    %                       'tipWidth':   width of the cone
    %
    %           Additionally, you can specify any patch object properties. (For
    %           example, you can make the arrow semitransparent by using
    %           'facealpha'.)
    %
    % example1: h = mArrow3([0 0 0],[1 1 1])
    %           (Draws an arrow from [0 0 0] to [1 1 1] with default properties.)
    %
    % example2: h = mArrow3([0 0 0],[1 1 1],'color','red','stemWidth',0.02,'facealpha',0.5)
    %           (Draws a red semitransparent arrow with a stem width of 0.02 units.)
    %
    % hint:     use light to achieve 3D impression
    % Author: Georg Stillfried

    propertyNames = {'edgeColor'};
    propertyValues = {'none'};

    %% evaluate property specifications
    for argno = 1:2:nargin-2
        switch varargin{argno}
            case 'color'
                propertyNames = {propertyNames{:},'facecolor'};
                propertyValues = {propertyValues{:},varargin{argno+1}};
            case 'stemWidth'
                if isreal(varargin{argno+1})
                    stemWidth = varargin{argno+1};
                else
                    warning('mArrow3:stemWidth','stemWidth must be a real number');
                end
            case 'tipWidth'
                if isreal(varargin{argno+1})
                    tipWidth = varargin{argno+1};
                else
                    warning('mArrow3:tipWidth','tipWidth must be a real number');
                end
            otherwise
                propertyNames = {propertyNames{:},varargin{argno}};
                propertyValues = {propertyValues{:},varargin{argno+1}};
        end
    end

    %% default parameters
    if ~exist('stemWidth','var')
        ax = axis;
        if numel(ax)==4
            stemWidth = norm(ax([2 4])-ax([1 3]))/300;
        elseif numel(ax)==6
            stemWidth = norm(ax([2 4 6])-ax([1 3 5]))/300;
        end
    end
    if ~exist('tipWidth','var')
        tipWidth = 3*stemWidth;
    end
    tipAngle = 22.5/180*pi;
    tipLength = tipWidth/tan(tipAngle/2);
    ppsc = 50;  % (points per small circle)
    ppbc = 250; % (points per big circle)

    %% ensure column vectors
    p1 = p1(:);
    p2 = p2(:);

    %% basic lengths and vectors
    x = (p2-p1)/norm(p2-p1); % (unit vector in arrow direction)
    y = cross(x,[0;0;1]);    % (y and z are unit vectors orthogonal to arrow)
    if norm(y)<0.1
        y = cross(x,[0;1;0]);
    end
    y = y/norm(y);
    z = cross(x,y);
    z = z/norm(z);

    %% basic angles
    theta = 0:2*pi/ppsc:2*pi; % (list of angles from 0 to 2*pi for small circle)
    sintheta = sin(theta);
    costheta = cos(theta);
    upsilon = 0:2*pi/ppbc:2*pi; % (list of angles from 0 to 2*pi for big circle)
    sinupsilon = sin(upsilon);
    cosupsilon = cos(upsilon);

    %% initialize face matrix
    f = NaN([ppsc+ppbc+2 ppbc+1]);

    %% normal arrow
    if norm(p2-p1)>tipLength
        % vertices of the first stem circle
        for idx = 1:ppsc+1
            v(idx,:) = p1 + stemWidth*(sintheta(idx)*y + costheta(idx)*z);
        end
        % vertices of the second stem circle
        p3 = p2-tipLength*x;
        for idx = 1:ppsc+1
            v(ppsc+1+idx,:) = p3 + stemWidth*(sintheta(idx)*y + costheta(idx)*z);
        end
        % vertices of the tip circle
        for idx = 1:ppbc+1
            v(2*ppsc+2+idx,:) = p3 + tipWidth*(sinupsilon(idx)*y + cosupsilon(idx)*z);
        end
        % vertex of the tiptip
        v(2*ppsc+ppbc+4,:) = p2;

        % face of the stem circle
        f(1,1:ppsc+1) = 1:ppsc+1;
        % faces of the stem cylinder
        for idx = 1:ppsc
            f(1+idx,1:4) = [idx idx+1 ppsc+1+idx+1 ppsc+1+idx];
        end
        % face of the tip circle
        f(ppsc+2,:) = 2*ppsc+3:(2*ppsc+3)+ppbc;
        % faces of the tip cone
        for idx = 1:ppbc
            f(ppsc+2+idx,1:3) = [2*ppsc+2+idx 2*ppsc+2+idx+1 2*ppsc+ppbc+4];
        end

    %% only cone v
    else
        tipWidth = 2*sin(tipAngle/2)*norm(p2-p1);
        % vertices of the tip circle
        for idx = 1:ppbc+1
            v(idx,:) = p1 + tipWidth*(sinupsilon(idx)*y + cosupsilon(idx)*z);
        end
        % vertex of the tiptip
        v(ppbc+2,:) = p2;
        % face of the tip circle
        f(1,:) = 1:ppbc+1;
        % faces of the tip cone
        for idx = 1:ppbc
            f(1+idx,1:3) = [idx idx+1 ppbc+2];
        end
    end

    %% draw
    fv.faces = f;
    fv.vertices = v;
    h = patch(fv);
    for propno = 1:numel(propertyNames)
        try
            set(h,propertyNames{propno},propertyValues{propno});
        catch
            disp(lasterr)
        end
    end
end