function export2stl(D,offset)
% exports a data matrix to an *.StL file that can be used for a 3D printer
% 
% input:
% D ....... (k x n)-matrix, with
%      k is the number of spectra, 
%      n is the number of channels
% 
% optional input:
% offset .... lift up D by the offset value (in % with respect to the maximum
%         height of D)
% 
% output  print.stl in a 10cm x 10cm x 5*(1+offset)cm scale

if ~exist('offset')
    offset = 1;
end

fileID = fopen('print.stl','w');
fprintf(fileID,'solid print \n');

maxD = max(max(D));
D = 1./maxD*D;
D = 50*(D+offset); % scale to 5cm

t = 1:size(D,1);
x = 1:size(D,2);

lenx = length(x);
lent = length(t);

x = x-min(x);
x = 100*x./max(max(x)); % scale to 10cm
t = t-min(t);
t = 100*t./max(max(t)); % scale to 10cm

%surface
for j=1:lent-1
    disp(j)
    for i=1:lenx-1
        %tri1
        p = [x(i)   t(j)   D(j,i)];
        q = [x(i+1) t(j)   D(j,i+1)];
        r = [x(i+1) t(j+1) D(j+1,i+1)];
        n = cross(q-p,r-p);
        maketri(p,q,r,n,fileID);
        
        %tri2
        p = [x(i)   t(j)   D(j,i)];
        q = [x(i+1) t(j+1) D(j+1,i+1)];
        r = [x(i)   t(j+1) D(j+1,i)];
        n = cross(q-p,r-p);
        maketri(p,q,r,n,fileID);
    end
end

%sides
for i=1:lenx-1
    %tri1
    p = [x(i)   t(1)   0];
    q = [x(i+1) t(1)   0];
    r = [x(i+1) t(1)   D(1,i+1)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri2
    p = [x(i)   t(1)   0];
    q = [x(i+1) t(1)   D(1,i+1)];
    r = [x(i)   t(1)   D(1,i)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri3
    p = [x(i)   t(end)   0];
    r = [x(i+1) t(end)   0];
    q = [x(i+1) t(end)   D(end,i+1)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri4
    p = [x(i)   t(end)   0];
    r = [x(i+1) t(end)   D(end,i+1)];
    q = [x(i)   t(end)   D(end,i)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
end

for i=1:lent-1
    %tri1
    p = [x(1) t(i)    0];
    q = [x(1) t(i+1)  D(i+1,1)];
    r = [x(1) t(i+1)  0];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri2
    p = [x(1) t(i)   0];
    q = [x(1) t(i)   D(i,1)];
    r = [x(1) t(i+1) D(i+1,1)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri3
    p = [x(end) t(i)    0];
    r = [x(end) t(i+1)  D(i+1,end)];
    q = [x(end) t(i+1)  0];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    %tri4
    p = [x(end) t(i)   0];
    r = [x(end) t(i)   D(i,end)];
    q = [x(end) t(i+1) D(i+1,end)];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
end

%bottom
m = [x(floor(lenx/2)) t(floor(lent/2)) 0];
for i=1:lent-1
    p = [x(1) t(i)    0];
    q = [x(1) t(i+1)  0];
    r = m;
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    p = [x(end) t(i)    0];
    q = m;
    r = [x(end) t(i+1)  0];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
end

for i=1:lenx-1
    p = [x(i)   t(1)  0];
    q = m;
    r = [x(i+1) t(1)  0];
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
    
    p = [x(i)   t(end)  0];
    q = [x(i+1) t(end)  0];
    r = m;
    n = cross(q-p,r-p);
    maketri(p,q,r,n,fileID);
end

fprintf(fileID,'endsolid print \n');
fclose(fileID);

function [  ] = maketri( p,q,r,n,fileID )
fprintf(fileID,'facet normal %8.5f %8.5f %8.5f\n',n(1),n(2),n(3));
fprintf(fileID,'outer loop\n');
fprintf(fileID,'vertex %8.5f %8.5f %8.5f\n',p(1),p(2),p(3)');
fprintf(fileID,'vertex %8.5f %8.5f %8.5f\n',q(1),q(2),q(3)');
fprintf(fileID,'vertex %8.5f %8.5f %8.5f\n',r(1),r(2),r(3)');
fprintf(fileID,'endloop\n');
fprintf(fileID,'endfacet\n');
