clear all
close all

% disegna toro
ls=30;
lphi=50;
[s,phi]=meshgrid(linspace(0,2*pi,ls),linspace(0,2*pi,lphi));
gx=cos(s);
gy=sin(s);
gz(1:lphi,1:ls)=0;
nx=-cos(s);
ny=-sin(s);
nz=gz;
bx=gz;
by=gz;
bz(1:lphi,1:ls)=1;
rho=0.4;
x=gx+rho*cos(phi).*nx+rho*sin(phi).*bx;
y=gy+rho*cos(phi).*ny+rho*sin(phi).*by;
z=gz+rho*cos(phi).*nz+rho*sin(phi).*bz;
surf(x,y,z)
axis equal
axis tight

figure
clear all

% disegna tubo con sezione circolare attorno ad elica
[s,phi]=meshgrid(linspace(0,6*pi,40),linspace(0,2*pi,10));
[ls,lphi]=size(s);
s2=s/sqrt(2);
gx=cos(s2);
gy=sin(s2);
gz=s2;
nx=-cos(s2);
ny=-sin(s2);
nz(1:ls,1:lphi)=0;
bx=sin(s2)/sqrt(2);
by=-cos(s2)/sqrt(2);
bz(1:ls,1:lphi)=1/sqrt(2);
rho=0.8;
x=gx+rho*cos(phi).*nx+rho*sin(phi).*bx;
y=gy+rho*cos(phi).*ny+rho*sin(phi).*by;
z=gz+rho*cos(phi).*nz+rho*sin(phi).*bz;
surf(x,y,z)
axis tight
axis equal

figure
clear all

% disegna tubo con sezione ellittica attorno ad elica
[s,phi]=meshgrid(linspace(0,6*pi,40),linspace(0,2*pi,20));
[lphi,ls]=size(s);
% matrice delle giuste dimensioni composta da zeri
z0(1:lphi,1:ls)=0;
% matrice delle giuste dimensioni composta da 1/sqrt(2)
isq2(1:lphi,1:ls)=1/sqrt(2);
s2=isq2.*s;
gx=cos(s2);
gy=sin(s2);
gz=s2;
nx=-cos(s2);
ny=-sin(s2);
nz=z0;
bx=isq2.*sin(s2);
by=-isq2.*cos(s2);
bz=isq2;
rho1=0.8;
rho2=0.3;
x=gx+rho1*cos(phi).*nx+rho2*sin(phi).*bx;
y=gy+rho1*cos(phi).*ny+rho2*sin(phi).*by;
z=gz+rho1*cos(phi).*nz+rho2*sin(phi).*bz;
surf(x,y,z)
axis tight
axis equal

clear all
figure

% disegna toro a sezione "a fiore"
[s,phi]=meshgrid(linspace(0,2*pi,30),linspace(0,2*pi,80));
[lphi,ls]=size(s);
z0(1:lphi,1:ls)=0;
o1(1:lphi,1:ls)=1;
gx=cos(s);
gy=sin(s);
gz=z0;
nx=-cos(s);
ny=-sin(s);
nz=z0;
bx=z0;
by=z0;
bz=o1;
rho=0.7*o1+0.1*sin(8*phi);
x=gx+rho.*cos(phi).*nx+rho.*sin(phi).*bx;
y=gy+rho.*cos(phi).*ny+rho.*sin(phi).*by;
z=gz+rho.*cos(phi).*nz+rho.*sin(phi).*bz;
surf(x,y,z)
axis equal

figure
clear all

% disegna nastro di Moebius
[t,s]=meshgrid(linspace(0,2*pi,30),linspace(-0.3,0.3,8));
gx=cos(t);
gy=sin(t);
gz(1:8,1:30)=0;
nx=-cos(t);
ny=-sin(t);
nz=gz;
bx=gz;
by=gz;
bz(1:8,1:30)=1;
x=gx+s.*cos(t/2).*bx+s.*sin(t/2).*nx;
y=gy+s.*cos(t/2).*by+s.*sin(t/2).*ny;
z=gz+s.*cos(t/2).*bz+s.*sin(t/2).*nz;
surf(x,y,z)
axis equal

