
n = 10;

h = 1;
b = 1;

[U,THETA] = meshgrid([0:h/(n-1):h],[0:(2*pi)/(n-1):2*pi]);

# kon
kx = @(u,theta) (b/2)*(h-u).*cos(theta);
ky = @(u,theta) (b/2)*(h-u).*sin(theta);
kz = @(u,theta) u;

ROTATION

alfa = pi/2;

# normal och kon-gemensam punkt till

# cirkel
a1 = [kx((8/10)*h,alfa) ky((8/10)*h,alfa) kz((8/10)*h,alfa)]';
pn1 = [0 0 1]'; 
pn1 = pn1 / norm(pn1);

# ellips
a2 = [kx((7/10)*h,alfa) ky((7/10)*h,alfa) kz((7/10)*h,alfa)]';
pn2 = [0 .3 1]'; 
pn2 = pn2 / norm(pn2);

# parabel
a3 = [kx((5/10)*h,alfa) ky((5/10)*h,alfa) kz((5/10)*h,alfa)]';
pn3 = [0 0.6 1]'; 
pn3 = pn3 / norm(pn3);

# hyperbel
a4 = [kx((3/10)*h,alfa) ky((3/10)*h,alfa) kz((3/10)*h,alfa)]';
pn4 = [0 0.7 1]'; 
pn4 = pn4 / norm(pn4);

# parametrar
[S,T] = meshgrid([-2:1/(n-1):2],[-2:1/(n-1):2]);

# vektorer parallella med planet
v1 = R(pi/2,0)*pn1;
v1 = v1 ./ norm(v1);
w1 = cross(v1,pn1);

# plan1
p1x = @(s,t) a1(1) + s*v1(1) + t*w1(1);
p1y = @(s,t) a1(2) + s*v1(2) + t*w1(2);
p1z = @(s,t) a1(3) + s*v1(3) + t*w1(3);

# vektorer i planet
v2 = R(pi/2,0)*pn2;
v2 = v2 ./ norm(v2);
w2 = cross(v2,pn2);

# plan2
p2x = @(s,t) a2(1) + s*v2(1) + t*w2(1);
p2y = @(s,t) a2(2) + s*v2(2) + t*w2(2);
p2z = @(s,t) a2(3) + s*v2(3) + t*w2(3);

# vektorer i planet
v3 = R(pi/2,0)*pn3;
v3 = v3 ./ norm(v3);
w3 = cross(v3,pn3);

# plan3
p3x = @(s,t) a3(1) + s*v3(1) + t*w3(1);
p3y = @(s,t) a3(2) + s*v3(2) + t*w3(2);
p3z = @(s,t) a3(3) + s*v3(3) + t*w3(3);

# vektorer i planet
v4 = R(pi/2,0)*pn4;
v4 = v4 ./ norm(v4);
w4 = cross(v4,pn4);

# plan3
p4x = @(s,t) a4(1) + s*v4(1) + t*w4(1);
p4y = @(s,t) a4(2) + s*v4(2) + t*w4(2);
p4z = @(s,t) a4(3) + s*v4(3) + t*w4(3);

# kon
K = zeros(n,n,3);

t = [0:h/(n-1):h];

for j=[1:n]
  for jj=[1:n]
    K(j,jj,:) =(b/2) * (t(j)) .* R(0,gamma(jj))*[1 0 0]' + [0 0 (h-t(j))]';
  endfor
endfor

# Konisk sektion 
omega = [0:2*pi/(n-1):2*pi];
gamma = [0:2*pi/(n-1):2*pi];

u1 = @(theta) (pn1(1)*(h.*cos(theta)-a1(1)) + pn1(2)*(h.*sin(theta)-a1(2)) - pn1(3).*a1(3)) ./ (pn1(1).*cos(theta) + pn1(2).*sin(theta) - pn1(3));
u2 = @(theta) (pn2(1)*(h.*cos(theta)-a2(1)) + pn2(2)*(h.*sin(theta)-a2(2)) - pn2(3).*a2(3)) ./ (pn2(1).*cos(theta) + pn2(2).*sin(theta) - pn2(3));
u3 = @(theta) (pn3(1)*(h.*cos(theta)-a3(1)) + pn3(2)*(h.*sin(theta)-a3(2)) - pn3(3).*a3(3)) ./ (pn3(1).*cos(theta) + pn3(2).*sin(theta) - pn3(3));
u4 = @(theta) (pn4(1)*(h.*cos(theta)-a4(1)) + pn4(2)*(h.*sin(theta)-a4(2)) - pn4(3).*a4(3)) ./ (pn4(1).*cos(theta) + pn4(2).*sin(theta) - pn4(3));

figure(1);
clf;
hold on;
surf(p1x(S,T),p1y(S,T),p1z(S,T));
plot3(kx(u1(omega),omega),ky(u1(omega),omega),kz(u1(omega),omega)+.1);


# surf(kx(u1(THETA),THETA),ky(u1(THETA),THETA),kz(u1(THETA),THETA));
# surf(kx(u2(THETA),THETA),ky(u2(THETA),THETA),kz(u2(THETA),THETA));
# surf(kx(u3(THETA),THETA),ky(u3(THETA),THETA),kz(u3(THETA),THETA));
# surf(kx(u4(THETA),THETA),ky(u4(THETA),THETA),kz(u4(THETA),THETA));

# plot3(a1(1),a1(2),a1(3),'*r');
# plot3(a2(1),a2(2),a2(3),'*r');
# plot3(a3(1),a3(2),a3(3),'*r');
# plot3(a4(1),a4(2),a4(3),'*r');

# axis([-1 1 -1 1 0 1]);
# view([1 1 1]);


