fc_plotPajek(net, netLabels, pajekDir, pajekGraphType)

function [pajekDir] = fc_plotPajek(net, netLabels, pajekDir, pajekGraphType)

% Function to call on pajek.exe (an external program, see
% http://vlado.fmf.uni-lj.si/pub/networks/pajek/) to create a 3-D plot
% pajekDir must be set to the directory that pajek.exe resides - this can
% be done using the menu. The function first writes a command script in the pajek
% directory (pajek.log) and a .bat file to run pajek (runPajek.bat) in the
% pajek directory, to instruct it to read this file and save a .mol file containing 3-D
% coordinates (pajek.mol). Then saves the connectivity data
% as a .net file, and is then read and plotted as a 3-D graph in
% a new figure window. There are 2 pajek graph types: Kamada and Fruchterman
% Reingold 3D, which can be selected using the pajek menu.

if isempty(net)
    errordlg('Cant plot network until network data has been loaded!')
    pajekDir = []
    return;
end

pajeksave(net, netLabels, [], 'pajek.net');

% Check the given pajekDir path for existence of a directory (7)
if exist(pajekDir) ~= 7
    errordlg('Pajek directory not found. Please select the directory the main Pajek executable resides in using the Pajek menu option.')
    pajekDir = []
    return;
end

%Make sure there's a runPajek.bat in the pajek directory containing:
% cd [path to pajek dir]
% pajek

if exist([pajekDir '\runPajek.log']) == 0
    runPajekFile = fopen([pajekDir '\runPajek.bat'], 'w');
    fprintf(runPajekFile, 'cd %s \n pajek', pajekDir);
    fclose(runPajekFile);
end

noOfNodes = num2str(length(net));

if strcmp(pajekGraphType, 'Kamada-Kawai') == 1
    graphType = 'KAMADA';
elseif strcmp(pajekGraphType, 'Fruchterman Reingold 3D') == 1
    graphType = 'FRUCHT3D';
else
    errordlg('Error - Couldnt ascertain pajek graph type');
    pajekDir = [];
    return;
end

currentDirectory = pwd;
pajekLogFile = fopen([pajekDir '\pajek.log'], 'w'); % w = create new file for writing, discard any existing file
pajekInstructions = {'NETBEGIN 1'; 'CLUBEGIN 1'; 'PERBEGIN 1'; 'CLSBEGIN 1'; 'HIEBEGIN 1';...
    'VECBEGIN 1'; ''; ['Msg Reading Network   ---    ' currentDirectory '\pajek.net'];...
    ['N 1 RDN ' currentDirectory '\pajek.net (' noOfNodes ')'];...
    'E 1 CIRCULAR'; 'E 1 DRAW 0 0 0';...
    ['E 1 ' graphType];...
    ['E 1 MDLMOL 0 0 0 ' currentDirectory '\pajek.mol 1.000']; 'EXIT'};

for printCount = 1:14
    fprintf(pajekLogFile, '%s \r', pajekInstructions{printCount});
end

fclose(pajekLogFile);

try
    dos([pajekDir '\runPajek.bat']);
catch
    errordlg('Error running Pajek!')
    set(gca, 'visible', 'off');
    pajekDir = []
    return;
end

%---------------------------------
% Extract data from .mol

fid = fopen('pajek.mol');
mol_data = {};

while 1
    line = fgetl(fid);
    if line == -1
        break;
    end
    mol_data{end+1, 1} = line;
end
fclose(fid);

%-------------------------------------------------------------
% A .mol file has, on its 4th line, a series of numbers showing how many
% nodes in the file, and how many edges connecting those nodes.
% If there are 6 characters, the first 3 are the no.of nodes, last 3 are
% no. of edges, e.g.:
% _64_120
% If there are 7 characters, first 3 = no. of nodes, last 4 = no. of edges.
% Can only assume .mol file not designed for such large numbers. e.g.:
% 7023073 means 702 nodes, 3073 edges

% Work this out:

number_of_chars = length(mol_data{4});
if number_of_chars == 6
    noOfNodes = str2num(mol_data{4}(1:3));
    noOfEdges = str2num(mol_data{4}(4:6));
elseif number_of_chars == 7
    noOfNodes = str2num(mol_data{4}(1:3));
    noOfEdges = str2num(mol_data{4}(4:7));
else
    errordlg(['Error processing .mol file - unexpected no. of characters on 4th line: ' number_of_chars])
end

xcoords = [];
ycoords = [];
zcoords = [];

for nodeCounter = 1:noOfNodes

    temp = (sscanf (mol_data{4 + nodeCounter}, '%f %f %f %*c'));
    xcoords(nodeCounter,1) = temp(1,1);
    ycoords(nodeCounter,1) = temp(2,1);
    zcoords(nodeCounter,1) = temp(3,1);

end

edgeBegin = [];
edgeEnd = [];

for edgeCounter = 1:noOfEdges
    edgeBegin(edgeCounter,1) = str2num(mol_data{4 + noOfNodes + edgeCounter}(1:3));
    edgeEnd(edgeCounter,1) = str2num(mol_data{4 + noOfNodes + edgeCounter}(4:6));
end

pajek_fig = figure;
plot3(xcoords, ycoords, zcoords, 'r.');
hold on;
xEdgeCoords = [];
yEdgeCoords = [];
zEdgeCoords = [];

for edgePlotCounter = 1:noOfEdges
    xEdgeCoords(1, 1) = xcoords(edgeBegin (edgePlotCounter));
    yEdgeCoords(1, 1) = ycoords(edgeBegin (edgePlotCounter));
    zEdgeCoords(1, 1) = zcoords(edgeBegin (edgePlotCounter));
    xEdgeCoords(2, 1) = xcoords(edgeEnd (edgePlotCounter));
    yEdgeCoords(2, 1) = ycoords(edgeEnd (edgePlotCounter));
    zEdgeCoords(2, 1) = zcoords(edgeEnd (edgePlotCounter));
    plot3(xEdgeCoords, yEdgeCoords, zEdgeCoords, 'b-');
end

hold off;
rotate3d(gca);