My code does not work for double helix on toroid

Hello,

I would appreciate if anyone could help me on this.

My code does not work for double helix on toroid: I want to control the radius of the smaller helix and the number of turns. I would appreciate any suggestions.

%HELIX ON HELIX_CYLINDRICAL
clear
clc
%
syms t r R h
assume(t,'real')
assume(R,'real')
assume(h,'real')
% CYLINDRICAL HELIX
x=h*t;
y=R*cos(t);
z=R*sin(t);

% % HELIX ON TORUS
% %       VARIABLES
% R=3;% MAJOR RADIUS
% r=1;% MINOR RADIUS
% n=6;% No. of loops
% x = (R + r*cos(n*t)).*cos(t);
% y = r*sin(n*t);
% z = (R + r*cos(n*t)).*sin(t);

% r_vec is the equation of the curve of interest: here helix
r_vec=[x;y;z];
dx=diff(x);
dy=diff(y);
dz=diff(z);
tan_vec=[dx;dy;dz];
nor_vec=diff(tan_vec)/simplify(norm(diff(tan_vec)));
binor_vec=simplify((cross( tan_vec, nor_vec))/norm(tan_vec));
%parametrized for the next curve
syms a
assume(a,'real')
%Number of turns of the second helix
n2=12;
u=n2*t;
S= r_vec + a.*nor_vec.*cos(u)+a.*binor_vec.*sin(u);
% VARIABLES FOR THE SECOND LEVEL HELIX
% a= radius 
h_num=1; R_num=3; a_num=0.5;
S1=subs(S, [h, R, a], [h_num, R_num, a_num]);
% %%Plotting
% R=6; r=2; n=5;
t1 = 0:pi/100:40*pi;
% u1 = 0:pi/15:2*pi;
% x1=R.*cos(t)+r.*cos(t).*cos(n.*t);
% y1=R*sin(t)+r*cos(n.*t).*sin(t);
% z1=r*sin(n.*t);
% x1=S1(1)
% y1=S1(2)
% z1=S1(3)
x2=double(subs(S1(1),t ,t1)); y2=double(subs(S1(2),t ,t1)); z2=double(subs(S1(3),t ,t1));
% x3=double(subs(x2,u ,u1)); y3=double(subs(y2,u ,u1)); z3=double(subs(z2,u ,u1));

%x2=subs(x1,[t, u], [t1, u1])         
% %subs(x, t)

plot3(x2,y2,z2)
axis equal
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')


%   HELIX ON TORUS DIRECT FORUMLA
%
t = 0:pi/500:40*pi;
%       VARIABLES
R=3;% MAJOR RADIUS
r=1;% MINOR RADIUS
n=6;% No. of loops

xt = (R + r*cos(n*t)).*cos(t);
yt = r*sin(n*t);
zt = (R + r*cos(n*t)).*sin(t);

plot3(xt,yt,zt)
axis equal
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')


%
%HELIX ON HELIX_CYLINDRICAL
clear
clc
%
syms t r R h
assume(t,'real')
assume(R,'real')
assume(h,'real')
% CYLINDRICAL HELIX
% x=h*t;
% y=R*cos(t);
% z=R*sin(t);

% % HELIX ON TORUS
% %       VARIABLES
R=3;% MAJOR RADIUS
r=1;% MINOR RADIUS
n=6;% No. of loops
x = (R + r*cos(n*t)).*cos(t);
y = r*sin(n*t);
z = (R + r*cos(n*t)).*sin(t);

% r_vec is the equation of the curve of interest: here helix
r_vec=[x;y;z];
dx=diff(x);
dy=diff(y);
dz=diff(z);
tan_vec=[dx;dy;dz];
nor_vec=diff(tan_vec)/simplify(norm(diff(tan_vec)));
binor_vec=simplify((cross( tan_vec, nor_vec))/norm(tan_vec));
%parametrized for the next curve
syms a
assume(a,'real')
%Number of turns of the second helix
n2=12;
u=n2*t;
S= r_vec + a.*nor_vec.*cos(u)+a.*binor_vec.*sin(u);
% VARIABLES FOR THE SECOND LEVEL HELIX
% a= radius 
h_num=1; R_num=3; a_num=0.5;
S1=subs(S, [h, R, a], [h_num, R_num, a_num]);
% %%Plotting
% R=6; r=2; n=5;
t1 = 0:pi/100:40*pi;
% u1 = 0:pi/15:2*pi;
% x1=R.*cos(t)+r.*cos(t).*cos(n.*t);
% y1=R*sin(t)+r*cos(n.*t).*sin(t);
% z1=r*sin(n.*t);
% x1=S1(1)
% y1=S1(2)
% z1=S1(3)
x2=double(subs(S1(1),t ,t1)); y2=double(subs(S1(2),t ,t1)); z2=double(subs(S1(3),t ,t1));
% x3=double(subs(x2,u ,u1)); y3=double(subs(y2,u ,u1)); z3=double(subs(z2,u ,u1));

%x2=subs(x1,[t, u], [t1, u1])         
% %subs(x, t)

plot3(x2,y2,z2)
axis equal
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')


%%
%   HELIX ON TORUS DIRECT FORUMLA
%
t = 0:pi/500:40*pi;
%       VARIABLES
R=3;% MAJOR RADIUS
r=1;% MINOR RADIUS
n=6;% No. of loops

xt = (R + r*cos(n*t)).*cos(t);
yt = r*sin(n*t);
zt = (R + r*cos(n*t)).*sin(t);

plot3(xt,yt,zt)
axis equal
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)')