add_cylinder Subroutine

    subroutine add_cylinder(me,v1,v2,radius,num_points,initial_cap,final_cap,&
                            initial_normal,final_normal,final_normal_used,initial_vector,final_initial_vector_used)

    implicit none

    class(stl_file),intent(inout)              :: me
    real(wp),dimension(3),intent(in)           :: v1                        !! coordinates of initial point
    real(wp),dimension(3),intent(in)           :: v2                        !! coordinates of final point
    real(wp),intent(in)                        :: radius                    !! radius of the cylinder
    integer,intent(in)                         :: num_points                !! number of point on the circle (>=3)
    logical,intent(in)                         :: initial_cap               !! add a cap plate to the initial point
    logical,intent(in)                         :: final_cap                 !! add a cap plate to the final point
    real(wp),dimension(3),intent(in),optional  :: initial_normal            !! outward normal vector for initial circle
    real(wp),dimension(3),intent(in),optional  :: final_normal              !! outward normal vector for final circle
    real(wp),dimension(3),intent(out),optional :: final_normal_used         !! outward normal vector for final circle
                                                                            !! actually used
    real(wp),dimension(3),intent(in),optional  :: initial_vector            !! vector to use to generate the initial
                                                                            !! circle (x_unit by default)
    real(wp),dimension(3),intent(out),optional :: final_initial_vector_used !! the initial vector used for the final
                                                                            !! cap to generate the points

    integer :: i  !! counter
    integer :: nc !! number of points on the circle
    real(wp),dimension(3) :: n0 !! normal vector for initial circle
    real(wp),dimension(3) :: nf !! normal vector for final circle
    real(wp),dimension(:,:),allocatable :: n0_cap_points !! points for the initial cap
    real(wp),dimension(:,:),allocatable :: nf_cap_points !! points for the final cap

    nc = max(3, num_points)

    ! compute the end unit vectors
    !
    !        1 _________2
    !        |          |
    !  n0 <--*----------*--> nf
    !        |          |
    !         ----------

    if (present(initial_normal)) then
        n0 = unit(initial_normal)
    else
        n0 = unit(v1-v2)
    end if
    if (present(final_normal)) then
        nf = unit(final_normal)
    else
        nf = unit(v2-v1)
    end if
    if (present(final_normal_used)) final_normal_used = nf ! return if necessary

    ! create the points on the initial cap (optionally add the plate)
    call me%generate_circle(v1,radius,n0,nc,initial_cap,n0_cap_points,initial_vector=initial_vector)

    ! create the points on the final cap (optionally add the plate)
    ! [use the same initial vector to sure that the plate will form a good cylinder]
    call me%generate_circle(v2,radius,nf,nc,final_cap,nf_cap_points,&
                            initial_vector=unit(n0_cap_points(:,1)-v1),cw=.true.)
    if (present(final_initial_vector_used)) final_initial_vector_used = unit(nf_cap_points(:,1)-v2)

    ! now connect the points to form the cylinder:
    !   1----2  nf
    !   |  / |
    !   | /  |
    !   1----2  n0
    do i = 1, nc-1
        call me%add_plate(n0_cap_points(:,i),n0_cap_points(:,i+1),nf_cap_points(:,i+1))
        call me%add_plate(n0_cap_points(:,i),nf_cap_points(:,i+1),nf_cap_points(:,i))
    end do
    ! last one:
    !   n----1  nf
    !   |  / |
    !   | /  |
    !   n----1  n0
   call me%add_plate(n0_cap_points(:,nc),n0_cap_points(:,1),nf_cap_points(:,1))
   call me%add_plate(n0_cap_points(:,nc),nf_cap_points(:,1),nf_cap_points(:,nc))

    end subroutine add_cylinder