function helix_points = helix(h_freq, h_scale, h_range, h_inc)

  t = (h_range(1):h_inc:h_range(2))';
  helix_points = [t h_scale(1)*cos(t*h_freq) h_scale(2)*sin(t*h_freq)];

//end helix

function helix_rad_vecs = helix_rad(h_freq, h_scale, h_range, h_inc)

  t = (h_range(1):h_inc:h_range(2))';
  helix_rad_vecs = [zeros(size(t,1),1) ...
                    h_scale(1)*cos(t*h_freq) h_scale(2)*sin(t*h_freq)];

//end helix_rad

function helix_tan_vecs = helix_tan(h_freq, h_scale, h_range, h_inc)

  t = (h_range(1):h_inc:h_range(2))';
  helix_tan_vecs = [ones(size(t,1),1) ...
                    -h_scale(1)*h_freq*sin(t*h_freq) ...
                     h_scale(2)*h_freq*cos(t*h_freq)];

//end helix_tan

function cp_vec = cross_product(vec_a, vec_b)

  cp_vec = [vec_a(:,2).*vec_b(:,3)-vec_a(:,3).*vec_b(:,2) ...
            vec_a(:,3).*vec_b(:,1)-vec_a(:,1).*vec_b(:,3) ...
            vec_a(:,1).*vec_b(:,2)-vec_a(:,2).*vec_b(:,1)];

//end cross_product

function norm_vecs = normalize_rows(vecs)

  norm_vecs = zeros(vecs);

  for i=1:size(vecs,1)
    norm_vecs(i,:) = vecs(i,:)/norm(vecs(i,:));
  end //for i

//end normalize_rows

function helix_demo

  xselect(); xbasc();

  t            = (0:0.1:3.9*%pi)';
  helix_1      = helix(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_tan_1  = helix_tan(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_rad_1  = helix_rad(1., [1. 1.], [0. 3.9*%pi], 0.1)/15.;
  helix_perp_1 = normalize_rows(cross_product(helix_rad_1,helix_tan_1))/15.;

  tube_plot(helix_1, helix_rad_1, helix_perp_1, 8, [28 29 30 31 31 30 29 28]);

//end helix_demo

function sum_2_helix_demo

  xselect(); xbasc();

  t            = (0:0.1:3.9*%pi)';
  helix_1      = helix(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_tan_1  = helix_tan(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_rad_1  = helix_rad(1., [1. 1.], [0. 3.9*%pi], 0.1);

  helix_2      = helix(2., [0.6 0.6], [0. 3.9*%pi], 0.1);
  helix_tan_2  = helix_tan(2., [0.6 0.6], [0. 3.9*%pi], 0.1);
  helix_rad_2  = helix_rad(2., [0.6 0.6], [0. 3.9*%pi], 0.1);

  sum_helix      = helix_1 + helix_2;
  sum_helix_tan  = helix_tan_1 + helix_tan_2;
  sum_helix_rad  = normalize_rows(helix_rad_1 + helix_rad_2*2.)/15.;
  sum_helix_perp = ...
    normalize_rows(cross_product(sum_helix_rad,sum_helix_tan))/15.;

  tube_plot(sum_helix, sum_helix_rad, sum_helix_perp, ...
            8, [28 29 30 31 31 30 29 28]);

//end sum_2_helix_demo

function sum_3_helix_demo

  xselect(); xbasc();

  t            = (0:0.1:3.9*%pi)';
  helix_1      = helix(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_tan_1  = helix_tan(1., [1. 1.], [0. 3.9*%pi], 0.1);
  helix_rad_1  = helix_rad(1., [1. 1.], [0. 3.9*%pi], 0.1);

  helix_2      = helix(2., [0.6 0.6], [0. 3.9*%pi], 0.1);
  helix_tan_2  = helix_tan(2., [0.6 0.6], [0. 3.9*%pi], 0.1);
  helix_rad_2  = helix_rad(2., [0.6 0.6], [0. 3.9*%pi], 0.1);

  helix_3      = helix(3., [0.5 0.5], [0. 3.9*%pi], 0.1);
  helix_tan_3  = helix_tan(3., [0.5 0.5], [0. 3.9*%pi], 0.1);
  helix_rad_3  = helix_rad(3., [0.5 0.5], [%pi 3.9*%pi+%pi], 0.1);

  sum_helix      = helix_1 + helix_2 + helix_3;
  sum_helix_tan  = helix_tan_1 + helix_tan_2 + helix_tan_3;
  sum_helix_rad  = normalize_rows(helix_rad_1 + helix_rad_2*2. + ...
                                  helix_rad_3*3.)/15.;
  sum_helix_perp = ...
    normalize_rows(cross_product(sum_helix_rad,sum_helix_tan))/15.;

  tube_plot(sum_helix, sum_helix_rad, sum_helix_perp, ...
            8, [28 29 30 31 31 30 29 28]);

//end sum_3_helix_demo
