1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
| program test
implicit none
integer :: k
double precision, dimension(3) :: u,v,w
interface
function vect_product(u,v)
double precision, dimension(3) :: vect_product
double precision,dimension(3),intent(in) :: u,v
end function vect_product
end interface
u(:) = (/ 1.d0, 0.d0, 0.d0 /)
v(:) = (/ 0.d0, 1.d0, 0.d0 /)
w(:) = (/ 0.d0, 0.d0, 0.d0 /)
w(:) = vect_product(u(:),v(:))
write(*,"(3F5.2)")(w(k),k=1,3)
end program test
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
function vect_product( u, v)
! calcule le produit vectoriel ( u x v )
! retourne un vecteur norme
implicit none
! la fonction :
double precision, dimension(3) :: vect_product
! arguments :
double precision, dimension(3), intent(in) :: u, v
! locales :
double precision, dimension(3) :: w
double precision :: norme
w(1) = u(2)*v(3) - u(3)*v(2)
w(2) = u(3)*v(1) - u(1)*v(3)
w(3) = u(1)*v(2) - u(2)*v(1)
norme = sqrt( sum( (w(:))**2 ) )
if( norme > 1.d-15 ) then
vect_product(:) = w(:) / norme
else
vect_product(:) = 0.d0
write(*,"('ATTENTION : norme du produit vectoriel nulle')")
end if
return
end function vect_product |
Partager