function makepv (POS,project,varargin)

% MAKEPARTIVIEWFILES
%
% makes OUTFILE.cf and associated files that partiview can use. 
%
% verbose = 0 gives no error messages other than Matlab's usual stuff
% verbose = 1 gives polite error messages (default)
% verbose = 2 gives less polite ones.
%
% If K > 3, redefine POS = [POS(:,1:3) POS] to make using warp command easier.
%
% Dinoj Surendran (dinoj@cs.uchicago.edu)
% http://people.cs.uchicago.edu/~dinoj/vis/partiview
%
% Modified Feb 22 to put arrows on directed edges
  
[N,K] = size (POS);            % represents N points in R^K (K is 3 for now)

OUTFILE = project;
OUTDIR = '.';              % place files in current directory by default

L=length(OUTFILE);
if (L>3) && (strcmpi(OUTFILE(L-2:L),'.cf'))   % remove .cf extension if present
  OUTFILE = OUTFILE(1:L-3); 
end

% set default values

LAB = {};          % assume no labels
INFO = {};         % assume no extra information at any point
A = 0;             % assume no attributes for any point 
ATTRIB = [];       % ditto
ATTRIBNAMES = {};  % assume no attributes have names
GROUPS = ones(N,1);% assume all particles are in 1 group
GROUPNAMES = {'class1'}; % assume all groups are g1,g2..,gN (no names), N=1 by default
G=1;               % G is number of groups... I suppose this really should be set to 1
PICS = {};         % Assume no pictures at points
COLORS = [];       % Assume default color map 
EDGES = [];        % Assume no points are joined by edges
DIREDGES = 0;      %  if 1, then edges are actually directed
SCALE = 1;         % Assume the input points have a good scale.

opts.psize = 10;
opts.ptsizemin = 10;              % what does this do? 
opts.ptsizemax = 10;
opts.slum=1;
opts.polyminmin = 0.1;
opts.polyminmax = 100;
opts.polysize = 0.03;
opts.lsize = 0.05;             % normalized later on...
opts.labelmin = 2;
opts.labelmax = 2;
opts.clipnear = 1;      
opts.clipfar = 1e8;
opts.fullscreen = 1;
opts.showpoints = 1;
opts.alpha = 0.990;
opts.fov = 45;
opts.screen = 'desktop';      % could be 'geowall' 
opts.eyeseparation = 0.005;   % for geowall, isnt used otherwise

logfile = '';
logfid = 1;
verbose = 1;
for i=1:length(varargin)
  if strcmpi('logfile',varargin{i})
     logfile = varargin{i+1}; 
  end
  if strcmpi('verbose',varargin{i})
     verbose = varargin{i+1}; 
  end
end
if length(logfile) > 0
  logfid = fopen(logfile,'wt');
  if logfid <= 0
    throwerr (verbose, sprintf('Cant open the error log file %s',logfile),1);
    throwerr (verbose, sprintf('Using screen as error log file'),1);
    logfid = 1;
  end
end

for i=1:length(varargin)
  if strcmpi(varargin{i},'OUTDIR')
     OUTDIR = varargin{i+1};
  end

  if strcmpi(varargin{i},'nofullscreen');
    opts.fullscreen = 0;
  end
  
  if strcmpi(varargin{i},'LAB') | strcmpi(varargin{i},'LABEL') | strcmpi(varargin{i},'LABELS')
     LAB = varargin{i+1};
     if isstr (LAB)
       if exist (LAB,'file')
         LAB = loadcell(LAB);
       else
         throwerr (verbose,sprintf ('The file %s, for storing text labels for each data point, does not exist',LAB),logfid);
       end
     end
     if N ~= length(LAB)
       throwerr(verbose,sprintf('Warning: The number of labels (%d) doesnt match the number of data points (%d)',length(LAB),N),logfid);
       throwerr(verbose,sprintf('Your labels have been ignored'),logfid);
       LAB = {};
     end
   end

  if strcmpi(varargin{i},'INFO')
     INFO = varargin{i+1};
     if isstr (INFO)
       if exist (INFO,'file')
         INFO = loadcell(INFO);
       else
         throwerr (verbose,sprintf ('The file %s, for storing extra information for each data point, does not exist',INFO),logfid);
       end
     end
     if N ~= length(INFO)
       throwerr(verbose,sprintf('Warning: The number of entries in INFO (%d) doesnt match the number of data points (%d)',length(INFO),N),logfid);
       throwerr(verbose,sprintf('Your information entries have been ignored'),logfid);
       INFO = {};
     end
   end

  if strcmpi(varargin{i},'ATTRIB')
     ATTRIB = varargin{i+1};
     if isstr (ATTRIB)
       if exist (ATTRIB,'file')
         ATTRIB = load(ATTRIB);
       else
         throwerr (verbose,sprintf ('The file %s, for storing attributes for each data point, does not exist',ATTRIB),logfid);
       end
     end
         
     if N ~= size(ATTRIB,1)
       throwerr(verbose,sprintf('Warning: The number of rows in ATTRIB (%d) doesnt match the number of data points (%d)',size(ATTRIB,1),N),logfid);
       throwerr(verbose,sprintf('The attributes you provided for your data points have been ignored'),logfid);
       ATTRIB = [];
     end
 
     A = size(ATTRIB,2);
   end

   if strcmpi(varargin{i},'ATTRIBNAMES')
     ATTRIBNAMES = varargin{i+1};
     if isstr (ATTRIBNAMES)
       if exist (ATTRIBNAMES,'file')
         ATTRIBNAMES = loadcell(ATTRIBNAMES);
       else
         throwerr (verbose,sprintf ('The file %s, for storing attribute names, does not exist',ATTRIBNAMES),logfid);
         throwerr (verbose,'Ignoring it',logfid);
       end
     end
   end     

   if strcmpi(varargin{i},'GROUPS')
     GROUPS = varargin{i+1};
     if isstr (GROUPS)
       if exist (GROUPS,'file')
         GROUPS = load(GROUPS);
       else
         throwerr (verbose,sprintf ('The file %s, that says which group each data point belongs to, does not exist',GROUPS),logfid);
         throwerr (verbose,'Assuming all points are in the same group',logfid);
	 GROUPS = ones(N,1);
       end 
     end
     if (max(abs(mod(GROUPS,1))) > 0)
         throwerr (verbose,'Some values of your GROUPS vector arent integers. I am rounding them off for you',logfid);
	 GROUPS = round(GROUPS);
     end
     
     if min(GROUPS) < 1
         throwerr (verbose,sprintf ('Some values of your GROUPS vector are less than 1. Warning: the corresponding points will not show up'),logfid);
     end	 
   end     
   

   if strcmpi(varargin{i},'GROUPNAMES')
     GROUPNAMES = varargin{i+1};
     if isstr (GROUPNAMES)
       if exist (GROUPNAMES,'file')
         GROUPNAMES = loadcell(GROUPNAMES);
       else
         throwerr (verbose,sprintf ('The file %s, that gives names for each group, does not exist',GROUPNAMES),logfid);
         throwerr (verbose,'Ignoring it, calling groups g1,g2,etc',logfid);
       end
     end
   end     

   if strcmpi(varargin{i},'PICS') | strcmpi(varargin{i},'PICTURES') | strcmpi(varargin{i},'PIC') | strcmpi(varargin{i},'PICTURE')
     PICS = varargin{i+1};
     if isstr (PICS)
       if exist (PICS,'file')
         PICS = loadcell(PICS);
       else
         throwerr (verbose,sprintf('The file %s, that defines pictures for each point, does not exist',PICS),logfid);
         throwerr (verbose,'Assuming no pictures exist',logfid);
       end
     end
   end     

   if strcmpi(varargin{i},'PICSIZE') | strcmpi(varargin{i},'PICTURESIZE')
     PICSIZE = varargin{i+1};
   end     

   if strcmpi(varargin{i},'COLORS') | strcmpi(varargin{i},'COLOR')
     COLORS = varargin{i+1};
     if isstr (COLORS)
       if exist (COLORS,'file')
         COLORS = loadcell(COLORS);
       else
         throwerr (verbose,sprintf ('The file %s, that defines a color mapping, does not exist',COLORS),logfid);
         throwerr (verbose,'Using default color map',logfid);
       end
     end
   end     

   if strcmpi(varargin{i},'EDGES')
     EDGES = varargin{i+1};
     if isstr (EDGES)
       if exist (EDGES,'file')
         EDGES = load (EDGES);
       else
         throwerr (verbose,sprintf ('The file %s, for storing edges, does not exist',ATTRIBNAMES),logfid);
         throwerr (verbose,'Ignoring it',logfid);
       end
     end
     if (length(EDGES)>0) & issparse(EDGES)
       [a,b]=find(EDGES);
       EDGES = [a b];
     end 
   end     

   if strcmpi(varargin{i},'DIREDGES')
     EDGES = varargin{i+1};
     DIREDGES = 1;
     if isstr (EDGES)
       if exist (EDGES,'file')
         EDGES = load (EDGES);
       else
         throwerr (verbose,sprintf ('The file %s, for storing edges, does not exist',ATTRIBNAMES),logfid);
         throwerr (verbose,'Ignoring it',logfid);
       end
     end
   end     
   
   
   if strcmpi(varargin{i},'SCALE')
     SCALE = varargin{i+1};
   end     

   if strcmpi(varargin{i},'SCREEN')
     opts.screen = varargin{i+1};
   end     

   if strcmpi(varargin{i},'geowall') |  strcmpi(varargin{i},'stereo')
     opts.screen = 'geowall';
   end
   
   if strcmpi(varargin{i},'EYESEPARATION') | strcmpi(varargin{i},'EYESEP')
     opts.eyeseparation = varargin{i+1};
   end     

end

if ~isdir(OUTDIR)
  success = mkdir (OUTDIR);
  if (0 == success)
    throwerr (verbose,sprintf('Failed to make new directory %s',OUTDIR),logfid);
    throwerr (verbose,sprintf('Placing output files in current directory'),logfid);
    OUTDIR = '.'  
  end
end

% Fill up attribute names if required

for a = 1 : length(ATTRIBNAMES)
  if ~length(ATTRIBNAMES{a})
     ATTRIBNAMES{a} = ['a' num2str(a)];
  end
end
for a = 1+length(ATTRIBNAMES) : A
  ATTRIBNAMES{a} = ['a' num2str(a)];
end

G = max(max(GROUPS),length(GROUPNAMES));
for a = 1 : G
  if (length(GROUPNAMES) < a) | (length(GROUPNAMES{a}) == 0)
    if length(EDGES)>0
      GROUPNAMES{a} = sprintf ('group%d',a);
    else
      GROUPNAMES{a} = sprintf ('g%d',a);
    end
  end
end

if length(COLORS) & (3 ~= size(COLORS,2))
  throwerr (verbose, 'The color mapping matrix (COLORS) you provided does not have three columns. Using default color mapping',logfid);
end

if ~length(COLORS)
  if G > 1
    COLORS = rand(G,3);
    COLORS = COLORS./repmat(sum(COLORS,2),1,size(COLORS,2));
  else
    % if attribute values are ordinal
    COLORS = horzcat ( [0:0.01:1]', [1:-0.01:0]', zeros(101,1) );
  end
end

if ~length(OUTDIR)
  OUTDIR = '.';
end
%%%%%%%%%% inputting done %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

POS = POS * SCALE;

CMAP0 = fopen ([OUTDIR '/' 'white.cmap'],'wt');
fprintf (CMAP0, '1\n1 1 1\n');
fclose(CMAP0);
CMAP0 = fopen ([OUTDIR '/' 'cyan.cmap'],'wt');
fprintf (CMAP0, '1\n0 1 1\n');
fclose(CMAP0);

cmapfile = [OUTDIR '/' OUTFILE '.cmap'];
CMAP = fopen (cmapfile,'wt');
fprintf (CMAP,'%d\n', size(COLORS,1));
for i=1:size(COLORS,1)
  fprintf (CMAP,'%0.4f %0.4f %0.4f\n',COLORS(i,1), COLORS(i,2), COLORS(i,3));
end
fclose(CMAP);

cffile = [OUTDIR '/' OUTFILE '.cf'];
CFID = fopen (cffile,'wt');

fprintf (CFID,'filepath +:.:./images:./data\n');            % can place .sgi picture files in images directory if wished

if 1==opts.fullscreen 
  fprintf (CFID, 'eval detach \n');
end

if length(EDGES)>0                    % should be a numedges x 2 array by now
 fprintf (CFID,'\nobject g1=edges\n');
 fprintf (CFID,'include %s_edges.speck\n', OUTFILE);
 fprintf (CFID, 'eval cmap cyan.cmap\n');
 fprintf (CFID, 'eval alpha 1\n');
 fprintf (CFID, '\n');

 EDGEFID = fopen ([OUTDIR '/' OUTFILE '_edges.speck'], 'wt');
 for n=1:size(EDGES,1)

   a = EDGES(n,1);
   b = EDGES(n,2);
   from = POS(a,1:3);
   to = POS(b,1:3);

   if DIREDGES
     fprintf (EDGEFID, '%s', diredge(from,to));
   else
     fprintf (EDGEFID, '%s', undiredge(from,to));
   end
   
 end
end

[R,r,s] = analyzepoints (POS(:,1:3), [0 0 0], 1);              % radius of sphere centered at (0,0,0) containing all points
%[R r s s/R]
opts.clipnear = min (opts.clipnear,r/5);
opts.clipfar = max (opts.clipfar,R*5);

for g=1:G
  gpfilename = sprintf ('%s_%d.speck', OUTFILE,g);
  GPFID = fopen ([OUTDIR '/' gpfilename],'wt');
  if length(LAB)
    labfilename = sprintf ('%s_%d.label', OUTFILE, g);
    LABFID = fopen ([OUTDIR '/' labfilename],'wt');
  end
  sizevar = '';
  groupvar = 0;
  
  if G > 1 
   indices = find (GROUPS == g);
  else
    indices = 1:N;          % required? Or equivalent?
  end

  if K > 3
    offset = K;
  else
    offset = 0;
  end

  for i=1:size(ATTRIB,2)
    if (length(ATTRIBNAMES) >= i) & (length(ATTRIBNAMES{i})>0)
      fprintf (GPFID,'datavar %d %s\n', i-1+offset,ATTRIBNAMES{i});
      if strcmpi (ATTRIBNAMES{i},'size')
         sizevar = ATTRIBNAMES{i};
      end
      if strcmpi (ATTRIBNAMES{i},'gp') | strcmpi (ATTRIBNAMES{i},'group') 
         groupvar = i;
         indices = find (ATTRIB(:,i)==i);
      end
    end
  end

  if length(PICS)>0
    fprintf (GPFID, 'texturevar %d\n', offset+size(ATTRIB,2));
  end

  if length(sizevar)>0
    fprintf (GPFID, 'lumvar %s\n', sizevar);            % lumvar or polylumvar ?
  end

  if length(PICS)
    for ii=1:length(PICS)
      fprintf (GPFID, 'texture -M %d %s\n', ii, PICS{ii});
    end
    ATTRIB = horzcat(ATTRIB, [1:length(PICS)]');
  end

  for ii=1:length(indices)         % output values for each point in this group 
    n=indices(ii);
    for k=1:3
      fprintf (GPFID,'%0.6f ',POS(n,k));
    end
    if K>3
      for k=1:K
        fprintf (GPFID,'%0.6f ',POS(n,k));
      end
    end
    for k = 1 : size (ATTRIB,2)
      fprintf (GPFID,'%0.5f ',ATTRIB(n,k));
    end
    if (length(INFO) >= n) 
     fprintf (GPFID,' # %s', INFO{n});
    end
    fprintf (GPFID, '\n');

    if length(LAB)
      for k=1:3
        fprintf (LABFID,'%0.6f ',POS(n,k));
      end
      fprintf (LABFID,' text %s \n',LAB{n});
    end
  end

  fclose (GPFID);
  if (length(LAB)>0)
    fclose (LABFID);
  end

  if length(EDGES)
    specknum = g+1;
  else
    specknum = g;
  end

  fprintf (CFID, '\nobject g%d', g+double(length(EDGES)>0));
  if (length (GROUPNAMES) >= g) & (length (GROUPNAMES{g}) > 0)
    fprintf (CFID, '=%s', GROUPNAMES{g});    
  end

  fprintf (CFID,'\n');
  fprintf (CFID, '  include %s\n', gpfilename);
  if length(LAB)>0
    fprintf (CFID, '  include %s\n', labfilename);
  end  

  fprintf (CFID,'eval psize %0.10f\n', opts.psize * R * R);   
  fprintf (CFID,'eval ptsize %0.10f %0.10f\n', opts.ptsizemin, opts.ptsizemax);

  if length(LAB)>0
    opts.lsize
    fprintf (CFID, 'eval lsize %0.10f\n',opts.lsize * R);
    fprintf (CFID, 'eval labelmin %0.10f %0.10f\n',opts.labelmin, opts.labelmax);
    fprintf (CFID, 'eval laxes off\n');
    fprintf (CFID, 'eval labels on\n');
  end    

  if opts.showpoints % & ~length(PICS)
    fprintf (CFID,'eval points on\n');
  else
    fprintf (CFID,'eval points off\n');
  end

  fprintf (CFID, 'eval slum %0.10f\n', opts.slum);

  fprintf (CFID, 'eval lum const 1\n');              % change to take into account const points

  if length(PICS)
     fprintf (CFID, 'eval poly on\n');
     fprintf (CFID, 'eval points off\n');
     fprintf (CFID, 'eval polysides 4\n');
     fprintf (CFID, 'eval polysize %0.10f\n', opts.polysize * s);
%     fprintf (CFID, 'eval polymin %0.10f %0.5f\n', opts.polyminmin, opts.polyminmax);    % no need to use this line
     fprintf (CFID, 'eval textures on\n');
     fprintf (CFID, 'eval alpha %0.3f\n', opts.alpha);
   end

   if G > 1                                % color each group separately
     fprintf (CFID, 'eval color rgb %0.3f %0.3f %0.3f\n', COLORS(g,1),COLORS(g,2), COLORS(g,3));
   else
     if A
       fprintf (CFID, 'eval cmap %s.cmap\n', OUTFILE);
       fprintf (CFID, 'eval color %d\n',offset);
     else
       fprintf (CFID, 'eval color rgb 1 1 1 \n');
     end
   end
end

if length(EDGES)
  fclose(EDGEFID);
end

fprintf (CFID, '\n');
fprintf (CFID, 'eval focalpoint on\n');
fprintf (CFID, 'eval fov %0.3f\n', opts.fov);
if strcmpi(opts.screen,'geowall')
  fprintf (CFID, 'eval stereo %0.10f crosseyed\n',-1*opts.eyeseparation);
else
  fprintf (CFID, 'eval stereo off\n');
end
fprintf (CFID, 'eval clip %0.10f %0.10f\n', opts.clipnear, opts.clipfar); 

if isfield(opts,'axes')
  fprintf (CFID, 'eval censize %0.10f\n', opts.axes);
else
  fprintf (CFID, 'eval censize %0.10f\n', 10^floor(log10(R)));
end

fprintf (CFID, 'warp -wx x:1 -wy y:1 -wz z:1\n'); 

if 1==opts.fullscreen 
  fprintf (CFID, 'eval detach full\n');
end


fprintf (CFID, 'eval jump 0 0 %0.10f \n', R*5);
fclose (CFID);

pvexec = 'partiview';
if isunix
  pvexec = ['./' pvexec];
else
  pvexec = [pvexec '.exe']; 
end

system (sprintf('%s %s &',pvexec,cffile));
if length(logfile) > 0
  fclose(logfid);
end

%%%%%%%%%% things to do when you have too much procrastination on your hands %%%%%%%%%

function throwerr (verbose,errmsg,logfid)

badwords = {'annoying person', 'twit', 'half-baked zinjanthropus', 'poor excuse for a mutation', 'poorly cooked cauliflower', 'rednecked rhinoceros', 'Oedipus', 'telligent beast', 'goddamn mudhole','salt-crested skunk', 'rat', 'louse', 'moron', 'slimy squirt', 'poor excuse for a mistake', 'outlier', 'monstrosity', 'rabbit pellet', 'piece of shit', 'blistering barracuda'};

%  documentation on some of the lesser known insults above:
%    telligent: opposite of intelligent, i.e. stupid
%    Oedipus: motherfucker
%    rabbit pellet: a piece of shit

rand ('state',sum(100*clock));
B = length(badwords);
b = min(B,max(1,ceil(B*rand)));
s = badwords{b};

if (verbose == 1)
  fprintf (logfid,'%s.\n',errmsg);
elseif (verbose == 2)
  fprintf (logfid,'%s, you %s.\n',errmsg,s);
end

%%%%%%%%%%%%%%%% compute initial jump parameters

function [R,r,s] = analyzepoints (positions, center,fracpoints)

% computes inner and outer radius of the sphere (centered at center) 
% R is radius of smallest sphere containing at least a fraction fracpoints of the points 
% r is radius of largest sphere containing at most a fraction 1-fracpoints of the points 
% fracpoints is a real number between 0 and 1 inclusive
% positions is a N x K matrix
% center is a 1 x K matrix
% Also returns the radius of the sphere 
% Also returns the a representative value for the separation of a point from its nearest neighbor
%
% Computations are based only on the first 3 dimensions of positions

R=0;
r = 0;
s = 0;
[N,K] = size(positions);
if K > 3
  positions = positions(:,1:3);
elseif K < 3
  positions = horzcat(positions,zeros(size(positions,1),3-K));
end

if ~N
  return;
end  

AA = positions - repmat(center,N,1);
AA = (sum(AA.^2,2)).^.5;
AA = sort (AA);
R = AA(ceil(fracpoints*N));
r = AA(1+ceil((1-fracpoints)*N));

blocksize = 100;

if N <= blocksize
  D2 = L2_distance (positions',positions');
  D2 = sort(D2(D2>0));
  s = D2 (ceil(0.01 * length(D2)));
else
  rand('state',sum(100*clock));

  % finds mean distance to their nearest neighbor of some randomly selected particles

  some = round(1+(N-1)*rand(blocksize,1));       % indices of some particles
  distnn = repmat(Inf,length(some),1);  

  numblocks = ceil(N/blocksize);
  for i = 1 : numblocks
    b1 = (i-1)*blocksize + 1;
    b2 = min(N,i*blocksize);
    E = L2_distance (positions(some,:)',positions(b1:b2,:)');     % E is length(some) x length(b1:b2)
    for j = 1 : length(some)
       x=E(j,:);
       distnn(j) = min(distnn(j),min(x(x>0)));
    end
  end
  s = mean(distnn);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function d = L2_distance(a,b,df)
% L2_DISTANCE - computes Euclidean distance matrix
%
% E = L2_distance(A,B)
%
%    A - (DxM) matrix 
%    B - (DxN) matrix
%    df = 1, force diagonals to be zero; 0 (default), do not force
% 
% Returns:
%    E - (MxN) Euclidean distances between vectors in A and B
%
%
% Description : 
%    This fully vectorized (VERY FAST!) m-file computes the 
%    Euclidean distance between two vectors by:
%
%                 ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
%
% Example : 
%    A = rand(400,100); B = rand(400,200);
%    d = distance(A,B);

% Author   : Roland Bunschoten
%            University of Amsterdam
%            Intelligent Autonomous Systems (IAS) group
%            Kruislaan 403  1098 SJ Amsterdam
%            tel.(+31)20-5257524
%            bunschot@wins.uva.nl
% Last Rev : Wed Oct 20 08:58:08 MET DST 1999
% Tested   : PC Matlab v5.2 and Solaris Matlab v5.3

% Copyright notice: You are free to modify, extend and distribute 
%    this code granted that the author of the original code is 
%    mentioned as the original author of the code.

% Fixed by JBT (3/18/00) to work for 1-dimensional vectors
% and to warn for imaginary numbers.  Also ensures that 
% output is all real, and allows the option of forcing diagonals to
% be zero.  

if (nargin < 2)
   error('Not enough input arguments');
end

if (nargin < 3)
   df = 0;    % by default, do not force 0 on the diagonal
end

if (size(a,1) ~= size(b,1))
   error('A and B should be of same dimensionality');
end

if ~(isreal(a)*isreal(b))
   disp('Warning: running distance.m with imaginary numbers.  Results may be off'); 
end

if (size(a,1) == 1)
  a = [a; zeros(1,size(a,2))]; 
  b = [b; zeros(1,size(b,2))]; 
end

aa=sum(a.*a); bb=sum(b.*b); ab=a'*b; 
d = sqrt(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab);

% make sure result is all real
d = real(d); 

% force 0 on the diagonal? 
if (df==1)
  d = d.*(1-eye(size(d)));
end

%%%%%%%%%%%%%%%%%

function z = loadcell (filename)

% LOADCELL reads an ascii file into a cell array, one line per element.
%          
% z = loadcell (filename);
%
% Note: anything after % in the file is not read

if ~ischar(filename), 
  error (sprintf ('The first argument to READLINE must be a string (representing a file name)',filename));
elseif ~exist (filename)
  error (sprintf ('File %s not found',filename));
end;

z = textread (filename, '%s', 'commentstyle', 'matlab', 'headerlines',0, 'delimiter', '\n' );

%%%%%%%%%%%% EDGE DRAWING FUNCTIONS %%%%%%%%%%%

function s = undiredge(from, to);
  
  s = sprintf ('mesh -c 0 -s wire {\n1 2\n%f %f %f\n%f %f %f\n}\n', from(1), from(2), from(3), to(1), to(2), to(3));

function s = diredge(from,to)

  if ~min(abs(from-to))       % if from equals to, then dont draw an edge! (add loops later)
    s = '';
    return;
  end
  
  theta = 5; % angle of the arrowcone in degrees
  theta = (pi/180) * theta; % conv to radians
  
  % from,to are 1 x 3 matrices

  h = from + (to-from)*0.8;   % end of arrowhead, i.e. draw headend->to
  
  % c is some point such that c-h is orthogonal to the edge/arrowaxis
  
  
  if min(abs(to(1:2)-h(1:2)))   % not parallel to z-axis
    c = [0 0 dot(h,to-h)/(to(3)-h(3))];
  elseif (to(1)~=h(1))
    c = h + [0 1 0];
  else
    c = h + [1 0 0];
  end
  
  % d is a vector from the arrowconebase's center to some point on the circumference of the arrowcone
  
  d = sqrt(dot(to-h,to-h)) * tan(theta) *(c-h)/sqrt(dot(c-h,c-h));
  
  a = h + d;
  b = h - d;

  c2 = h+cross(c-h,to-h);

  if (0 == dot(c2-h,c2-h))
    from
    to
    c-h
    to-h
    c
    h
    c2
  end
  
  d2 = sqrt(dot(to-h,to-h)) * tan(theta) *(c2-h)/sqrt(dot(c2-h,c2-h));

  a2 = h+d2;
  b2 = h-d2;
  
  % draw first arrow
  
  s = sprintf ('mesh -c 0 -s solid {\n2 2\n%f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n}\n', to(1), to(2), to(3),  to(1), to(2), to(3), b(1), b(2), b(3), a(1), a(2), a(3));

  % draw second arrow
  
  s= sprintf ('mesh -c 0 -s solid {\n2 2\n%f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n}\n%s\n', to(1), to(2), to(3),  to(1), to(2), to(3), b2(1), b2(2), b2(3), a2(1), a2(2), a2(3),s);
  
  % draw arrow base
  s= sprintf ('mesh -c 0 -s solid {\n2 2\n%f %f %f\n%f %f %f\n%f %f %f\n%f %f %f\n}\n%s\n', a(1),a(2),a(3), a2(1),a2(2),a2(3), b2(1),b2(2),b2(3), b(1),b(2),b(3), s);
  
  % draw shank i.e. edge
  
  s = sprintf ('mesh -c 0 -s wire {\n1 2\n%f %f %f\n%f %f %f\n}\n%s', from(1), from(2), from(3), to(1), to(2), to(3), s);