unitVectors Subroutine

  subroutine unitVectors(params,Xo,F,b1,b2,b3,flag)
    !! @note Subrotuine that calculates an orthonormal basis using information 
    !! of the (local) magnetic field at position \(\mathbf{X}_0\). @endnote
    TYPE(KORC_PARAMS), INTENT(IN)                                      :: params
    !! Core KORC simulation parameters.
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(IN)                  :: Xo
    !! Array with the position of the simulated particles.
    TYPE(FIELDS), INTENT(IN)                                           :: F
    !! F An instance of the KORC derived type FIELDS.
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT)               :: b1
    !! Basis vector pointing along the local magnetic field, 
    !! that is, along \(\mathbf{b} = \mathbf{B}/B\).
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT)               :: b2
    !!  Basis vector perpendicular to b1
    REAL(rp), DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT)               :: b3
    !! Basis vector perpendicular to b1 and b2.
    INTEGER(is), DIMENSION(:), ALLOCATABLE, OPTIONAL, INTENT(INOUT)    :: flag
    !! Flag for each particle to decide whether it is being 
    !! followed (flag=T) or not (flag=F).
    TYPE(PARTICLES)                                                    :: vars
    !! A temporary instance of the KORC derived type PARTICLES.
    INTEGER                                                            :: ii
    !! Iterator.
    INTEGER                                                            :: ppp
    !! Number of particles.

    ppp = SIZE(Xo,1) ! Number of particles

    ALLOCATE( vars%X(ppp,3) )
    ALLOCATE( vars%Y(ppp,3) )
    ALLOCATE( vars%B(ppp,3) )
    ALLOCATE( vars%gradB(ppp,3) )
    ALLOCATE( vars%curlb(ppp,3) )
    ALLOCATE( vars%PSI_P(ppp) )
    ALLOCATE( vars%E(ppp,3) )
    ALLOCATE( vars%flag(ppp) )

    vars%X = Xo
    vars%flag = 1_idef
    vars%B=0._rp
    vars%PSI_P=0._rp


    call init_random_seed()

    call get_fields(params,vars,F)

    !    write(6,'("Bx: ",E17.10)') vars%B(:,1)
    !    write(6,'("By: ",E17.10)') vars%B(:,2)
    !    write(6,'("Bz: ",E17.10)') vars%B(:,3)

    do ii=1_idef,ppp
       if ( vars%flag(ii) .EQ. 1_idef ) then
          b1(ii,:) = vars%B(ii,:)/sqrt(vars%B(ii,1)*vars%B(ii,1)+ &
               vars%B(ii,2)*vars%B(ii,2)+vars%B(ii,3)*vars%B(ii,3))

          b2(ii,:) = cross(b1(ii,:),(/0.0_rp,0.0_rp,1.0_rp/))
          b2(ii,:) = b2(ii,:)/sqrt(b2(ii,1)*b2(ii,1)+b2(ii,2)*b2(ii,2)+ &
               b2(ii,3)*b2(ii,3))

          b3(ii,:) = cross(b1(ii,:),b2(ii,:))
          b3(ii,:) = b3(ii,:)/sqrt(b3(ii,1)*b3(ii,1)+b3(ii,2)*b3(ii,2)+ &
               b3(ii,3)*b3(ii,3))
       end if
    end do

    if (PRESENT(flag)) then
       flag = vars%flag
    end if

    DEALLOCATE( vars%X )
    DEALLOCATE( vars%Y )
    DEALLOCATE( vars%B )
    DEALLOCATE( vars%PSI_P )
    DEALLOCATE( vars%gradB )
    DEALLOCATE( vars%curlb )
    DEALLOCATE( vars%E )
    DEALLOCATE( vars%flag )
  end subroutine unitVectors