Source code changes/additions to Parametric_Line.m
to form Parametric_3D_Line.m:

n_q = 5 ; % number of quadrature points required

Total_first_z = 0; Total_sec_zz = 0; Total_sec_xz = 0; Total_sec_zy = 0 ;

  xy_e (1:n_n, 3) = z(e_nodes)            ; % z coord at el nodes
  first_z = 0 ; sec_zz  = 0 ; sec_xz  = 0 ; sec_zy  = 0 ; % initialize

%
%   or /clear/www/htdocs/mech517/Matlab_Plots/mesh_3d_plot.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/mesh_3d_shrink.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/L3_3d_mesh_plot.m

    dz_dr = DLH_q * xy_e (:, 3)                 ;  % dz (r) / dr
    J_det = sqrt ( dx_dr^2 + dy_dr^2 + dz_dr^2) ;  % dL (r) / dr
    first_z = first_z + xy_q (3) * J_det * w_q (q)                ; % 1st moment
    sec_zz  = sec_zz  + (xy_q (2)^2 + xy_q (1)^2) * J_det * w_q(q); % 2nd moment
    sec_xz  = sec_xz  - xy_q (1) * xy_q (3) * J_det * w_q(q)     ; % 2nd product
    sec_zy  = sec_zy  - xy_q (3) * xy_q (2) * J_det * w_q(q)     ; % 2nd product

  Total_first_z = Total_first_z + first_z  ; % update system sums
  Total_sec_zz  = Total_sec_zz  + sec_zz   ; % update system sums
  Total_sec_xz  = Total_sec_xz  + sec_xz   ; % update system sums
  Total_sec_zy  = Total_sec_zy  + sec_zy   ; % update system sums

z_bar = Total_first_z / Total_length ; % centroid coordinate
fprintf ('Z-centroid point = %g \n', z_bar)  
fprintf ('Global I_zz      = %g \n', Total_sec_zz)
fprintf ('Global I_xz      = %g \n', Total_sec_xz)
fprintf ('Global I_zy      = %g \n', Total_sec_zy)

OPTIONAL SOURCE CHANGES:
%  find second moments wrt the centroid (parallel axis theorem)

CG_sec_xx = Total_sec_xx - (y_bar^2 + z_bar^2) * Total_length ;
CG_sec_yy = Total_sec_yy - (x_bar^2 + z_bar^2) * Total_length ;
CG_sec_zz = Total_sec_zz - (x_bar^2 + y_bar^2) * Total_length ;
CG_sec_xy = Total_sec_xy - y_bar * x_bar * Total_length ;
CG_sec_xz = Total_sec_xz - z_bar * x_bar * Total_length ;
CG_sec_zy = Total_sec_zy - y_bar * z_bar * Total_length ;
fprintf ('\n')
fprintf ('Centroidal I_xx  = %g \n', CG_sec_xx)
fprintf ('Centroidal I_yy  = %g \n', CG_sec_yy)
fprintf ('Centroidal I_zz  = %g \n', CG_sec_zz)
fprintf ('Centroidal I_xy  = %g \n', CG_sec_xy)
fprintf ('Centroidal I_xz  = %g \n', CG_sec_xz)
fprintf ('Centroidal I_zy  = %g \n', CG_sec_zy)

RESULTS FROM SUPPLIED DATA:
>Parametric_Line_3D

(Echo of msh_remarks.tmp) 
============== Begin Application Remarks ====================
Numerical evaluation of the geometric properties of a
planar or space curve. Here a variable helix along the z-axis
is defined by a circular arc (R=4) coupled with z = 2x + 3y.
Nodal dof: x-,  y-, and z-coordinate 
Element type = 1 (a parametric line)
Element connection: two or three or four nodes per element
Element properties: none (later add)
Space Curve Properties for 9 node nesh:
Mesh type          Four L3  Eight L2
Length           = 11.3564  11.3094
X-centroid point = 2.76228  2.75691
Y-centroid point = 2.20195  2.19818
Z-centroid point = 12.1304  12.1084
=============== End Application Remarks =======================

User control variables 
n_g number of DOF per node        = 1 
n_n number of nodes per element   = 2 
n_p dimension of parametric space = 1 
n_q number of quadrature points   = 5 
n_s dimension of physical space   = 1 
Read 9 nodes. 
(Echo node, file msh_bc_xyz.tmp) 
node, bc_flags, 3 coordinates 
1, 0    4 0 8 
2, 0    3.9231 0.7804 10.1874 
3, 0    3.6955 1.5307 11.9831 
4, 0    3.3259 2.2223 13.3187 
5, 0    2.8284 2.8284 14.142 
6, 0    2.2223 3.3259 14.4223 
7, 0    1.5307 3.6955 14.1479 
8, 0    0.7804 3.9231 13.3301
9, 0    0 4 12 
WARNING: Read 3 coordinates, but expected 1 

Read 4 elements with type & 3 nodes each. 
(Echo element number, file msh_typ_nodes.tmp) 
1, 1 1 2 3    
2, 1 3 4 5    
3, 1 5 6 7    
4, 1 7 8 9    
WARNING Read 3 nodes per element but expected 2 

Space Curve Properties: 
Length           = 11.3564 
X-centroid point = 2.76228 
Y-centroid point = 2.20195 
Z-centroid point = 12.1304 
Global I_xx      = 1787.27 
Global I_yy      = 1817.79 
Global I_zz      = 181.692 
Global I_xy      = -50.5838 
Global I_xz      = -363.971 
Global I_zy      = -327.914 

Centroidal I_xx  = 61.1541 
Centroidal I_yy  = 60.0927 
Centroidal I_zz  = 39.9786 
Centroidal I_xy  = -119.658 
Centroidal I_xz  = -744.495 
Centroidal I_zy  = -631.248