% FabDiagram: contoured fabric diagram

clear all
rad = pi/180;
deg = 180/pi;

%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')

%load file with p/t data

filename = input(' Enter name of data file: ','s');
data = load(filename);
p = data(:,1);
t = data(:,2);
Ndata = length(p);  %Ndata = number of data points

%copy data point near primitive into upper hemisphere

limit = 30;
n = Ndata;
for i = 1:n
    if p(i) <= limit
        n = n + i;
        p(n) = -p(i);
        t(n) =  t(i) + 180;
        if t(n) > 360
           t(n) = t(n) - 360;
        end
    end
end 

%clean up p/t list

newP = [];
locP = find(p);   %locates non-zero elements in p
for j = 1:length(locP);
    newP = [newP  p(locP(j))];
end; 
p = newP;

newT = [];
locT = find(t);   %locates non-zero elements in t
for j = 1:length(locT);
    newT = [newT  t(locT(j))];
end;
t = newT;

Ntotal = length(p);  %Ntotal = number of augmented data points

%direction cosines of data vectors

for i = 1:Ntotal    
    Ld(i) = cos(p(i)*rad)*cos(t(i)*rad);
    Md(i) = cos(p(i)*rad)*sin(t(i)*rad);
    Nd(i) = sin(p(i)*rad);
end

%p/t of node vectors by direct calculation

clear p
clear t

n = 1;
p(n) = 90;
t(n) = 0;

for i = 1:10
    m = 6*i;
    radius = i/10;
    DeltaPhi = 360/m;
    for j = 1:m
        phi = j*DeltaPhi;
        n = n+1;         
        t(n) = phi;         
        theta = 2 * asin(radius/sqrt(2));
        p(n) = 90 - (theta*deg);
    end      
end 

%direction cosines (l,m,n) of node vectors

for i = 1:331     
     Ln(i) = cos(p(i)*rad) * cos(t(i)*rad);
     Mn(i) = cos(p(i)*rad) * sin(t(i)*rad);
     Nn(i) = sin(p(i)*rad);
end 

%(x,y) coordinates of nodes from p/t

for n = 1:331
    theta = (90-p(n))*rad;
    d = sqrt(2) * sin(theta/2);
    x(n) = d .* sin(t(n)*rad);
    y(n) = d .* cos(t(n)*rad);
end

%define counting cone

key = 2.0;
cone = (acos(Ndata/(Ndata+key^2)))*deg;

%count point density at all 331 nodes

for i = 1:331  
    z(i) = 0;

    %dot product of node vectors and data vectors
    
    for j = 1:Ntotal  
        theta = (acos(Ln(i)*Ld(j) + Mn(i)*Md(j) + Nn(i)*Ld(j)))*deg;         
        if theta <= cone
            z(i) = z(i)+1;
        end
    end
end

Zmax = 0;    
for j = 1:331
    z(j) = (z(j)/Ndata)*100;  %convert counts of percent
    if z(j) > Zmax
        Zmax = z(j);
    end    
end 

%for i = 1:331
  %  fprintf('%3i',i);
  %  fprintf('  %.5f',x(i));
 %   fprintf('  %.5f',y(i));
  %  fprintf('  %3i\n\n',z(i));
%end   
%fprintf(' Zmax = %.5f',Zmax);

plot(x,y,'k.')






%COntouring Kalsbeek
%set parameters:
%number of x and y nodes for grid
%make these smaller if you get an error about excceding the max array size
nx = 50;
ny = nx;
%contour interval:
ci = 2; 

%Check this assignment for the columns
n=x;
e=y;
hh=z;



%Grid the data linearly:

estep = abs((min(e)-max(e))/nx);
nstep = abs((min(n)-max(n))/ny);
ee = [min(e):estep:max(e)];
nn = [min(n):nstep:max(n)];
[xi,yi,zi] = griddata(e,n,hh,ee,nn');

v = floor(min(hh)):ci:max(hh);


%Here is contour map:
figure(2)
clf
%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')



[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
colormap(jet)
hold on

%show data:
plot(e,n,'o')

%show grid:
plot(xi,yi,'+')
title('Grid and observation points')
axis equal

print -depsc ptsandgrid.eps


%nice contour map
figure(3)
clf
%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')

[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
%use this one if you want to label the contours manually:
%clabel(c,h, 'manual');
colormap([0 0 0])
hold on

%show data:
plot(e,n,'.')

title('Linear interpolation and observation points')
axis equal

print -depsc contourslin.eps

%triangles for the linear interpolation
figure(4)
clf
%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')

tri = delaunay(e,n);
trimesh(tri,e,n,hh)
colormap([0 0 0])
view(2)
title('Triangles for linear interpolation')
axis equal
print -depsc triangles.eps



%nice contour map CUBIC interpolation
figure(5)
clf
%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')

[xi,yi,zi] = griddata(e,n,hh,ee,nn', 'cubic');
[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
%use this one if you want to label the contours manually:
%clabel(c,h, 'manual');
colormap([0 0 0])
hold on

%show data:
plot(e,n,'.')

title('Cubic interpolation and observation points')
axis equal

print -depsc contourscubic.eps


%imagesc
figure(6)
clf
plot(0,0,'k+')
hold on
%draw primitive circle

radius = 1;

theta = linspace(0,2*pi,180);
x = radius * cos(theta);
y = radius * sin(theta);
plot(x,y,'k-')
hold on
axis equal
axis off
plot(0,0,'k+')
plot(radius,0,'k+')
plot(0,radius,'k+')
plot(-radius,0,'k+')
plot(0,-radius,'k+')

imagesc(ee',nn',zi)
%show data:
plot(e,n,'k.')

title('Cubic interpolation and observation points')
axis equal
print -depsc image.eps  


%surface
figure(7)
clf
surfl(ee,nn,zi)
shading interp
colormap(gray)
hold on
%plot3(e,n,hh,'b.')
view(-30,60)
zlabel('elevation')   
title('Cubic interpolation')
%print -depsc vmsurface.eps 
