unitVectors Subroutine

  subroutine unitVectors(params,Xo,F,b1,b2,b3,flag,cart,hint)
    !! @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(C_PTR), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)    :: hint
    !! 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.
    LOGICAL :: cart
    REAL(rp), DIMENSION(3) ::b1tmp,b2tmp,b3tmp,tmpvec

!    write(output_unit_write,*) 'in unitVector'
    
    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%flagCon(ppp) )
    ALLOCATE( vars%initLCFS(ppp) )

#ifdef FIO
    ALLOCATE( vars%hint(ppp) )
#endif
    
    vars%X = Xo
#ifdef FIO
    vars%hint = hint
#endif
    vars%flagCon = flag
    vars%initLCFS = 0_is
    vars%B=0._rp
    vars%PSI_P=0._rp
    vars%cart=.false.
    
    !write(output_unit_write,*) 'before init_random_seed'
    
    call init_random_seed()

   ! write(output_unit_write,*) 'before get_fields'

    !write(6,*) 'before first get fields'
    call get_fields(params,vars,F)
    !write(6,*) 'before second get fields'

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

        !write(output_unit_write,*) 'before b1,b2,b3 calculation'

    tmpvec=(/1.0_rp,1.0_rp,1.0_rp/)
    
    do ii=1_idef,ppp
       !write(6,*) 'ii',ii
       if ( vars%flagCon(ii) .EQ. 1_idef ) then
          b1tmp = 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))

          b2tmp = cross(b1tmp,tmpvec)
          b2tmp = b2tmp/sqrt(b2tmp(1)*b2tmp(1)+b2tmp(2)*b2tmp(2)+ &
               b2tmp(3)*b2tmp(3))

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

    !write(output_unit_write,*) 'before copying hint and flag'
#ifdef FIO
    hint = vars%hint
#endif
    
    if (PRESENT(flag)) then
       flag = vars%flagCon
    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%flagCon )
#ifdef FIO
    DEALLOCATE( vars%hint)
#endif
    
    !write(output_unit_write,*) 'out unitVectors'
    
  end subroutine unitVectors