scale_vector Module Reference

module vector More...

Functions/Subroutines
subroutine, public	vectr_xyz2latlon (x, y, z, lat, lon)
subroutine, public	vectr_latlon2xyz (lat, lon, x, y, z, radius)
subroutine, public	vectr_cross (nv, a, b, c, d)
	exterior product of vector a->b and c->d More...
subroutine, public	vectr_dot (l, a, b, c, d)
	interior product of vector a->b and c->d More...
subroutine, public	vectr_abs (l, a)
	length of vector o->a More...
subroutine, public	vectr_angle (angle, a, b, c)
	calc angle between two vector(b->a,b->c) More...
subroutine, public	vectr_intersec (ifcross, p, a, b, c, d)
	judge intersection of two vector More...
subroutine, public	vectr_anticlockwise (vertex, nvert)
	bubble sort anticlockwise by angle More...
real(rp) function, public	vectr_triangle (a, b, c, polygon_type, radius)
	calc triangle area More...
real(rp) function, public	vectr_triangle_plane (a, b, c)
	calc triangle area on plane More...
subroutine, public	vectr_rotation (a, angle, iaxis)
	Apply rotation matrix. More...
subroutine, public	vectr_distance (r, lon1, lat1, lon2, lat2, dist)
	Get horizontal distance on the sphere. More...

Variables
integer, parameter, public	i_xaxis = 1
integer, parameter, public	i_yaxis = 2
integer, parameter, public	i_zaxis = 3

Detailed Description

module vector

Description

module for 3D vector on the sphere

Author

NICAM developers, Team SCALE

Function/Subroutine Documentation

vectr_xyz2latlon()

subroutine, public scale_vector::vectr_xyz2latlon	(	real(rp), intent(in)	x,
		real(rp), intent(in)	y,
		real(rp), intent(in)	z,
		real(rp), intent(out)	lat,
		real(rp), intent(out)	lon
	)

Definition at line 64 of file scale_vector.F90.

    use scale_const, only: &
       eps => const_eps
    implicit none
 
    real(RP), intent(in)  :: x
    real(RP), intent(in)  :: y
    real(RP), intent(in)  :: z
    real(RP), intent(out) :: lat
    real(RP), intent(out) :: lon
 
    real(RP) :: length, length_h
    !---------------------------------------------------------------------------
 
    length = sqrt( x*x + y*y + z*z )
 
    if ( length < eps ) then ! 3D vector length is
       lat = 0.0_rp
       lon = 0.0_rp
       return
    endif
 
    if    ( z / length >= 1.0_rp ) then ! vector is parallele to z axis.
       lat = asin(1.0_rp)
       lon = 0.0_rp
       return
    elseif( z / length <= -1.0_rp ) then ! vector is parallele to z axis.
       lat = asin(-1.0_rp)
       lon = 0.0_rp
       return
    else
       lat = asin( z / length )
    endif
 
    length_h = sqrt( x*x + y*y )
 
    if ( length_h < eps ) then
       lon = 0.0_rp
       return
    endif
 
    if    ( x / length_h >= 1.0_rp ) then
       lon = acos(1.0_rp)
    elseif( x / length_h <= -1.0_rp ) then
       lon = acos(-1.0_rp)
    else
       lon = acos( x / length_h )
    endif
 
    if( y < 0.0_rp ) lon = -lon
 
    return

References scale_const::const_eps.

vectr_latlon2xyz()

subroutine, public scale_vector::vectr_latlon2xyz	(	real(rp), intent(in)	lat,
		real(rp), intent(in)	lon,
		real(rp), intent(out)	x,
		real(rp), intent(out)	y,
		real(rp), intent(out)	z,
		real(rp), intent(in)	radius
	)

Definition at line 125 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(in)  :: lat
    real(RP), intent(in)  :: lon
    real(RP), intent(out) :: x
    real(RP), intent(out) :: y
    real(RP), intent(out) :: z
    real(RP), intent(in)  :: radius
    !---------------------------------------------------------------------------
 
    x = radius * cos(lat) * cos(lon)
    y = radius * cos(lat) * sin(lon)
    z = radius * sin(lat)
 
    return

vectr_cross()

subroutine, public scale_vector::vectr_cross	(	real(rp), dimension(3), intent(out)	nv,
		real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c,
		real(rp), dimension(3), intent(in)	d
	)

exterior product of vector a->b and c->d

Definition at line 145 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(out) :: nv(3)                  ! normal vector
    real(RP), intent(in)  :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
    !---------------------------------------------------------------------------
 
    nv(1) = ( b(2)-a(2) ) * ( d(3)-c(3) ) &
          - ( b(3)-a(3) ) * ( d(2)-c(2) )
    nv(2) = ( b(3)-a(3) ) * ( d(1)-c(1) ) &
          - ( b(1)-a(1) ) * ( d(3)-c(3) )
    nv(3) = ( b(1)-a(1) ) * ( d(2)-c(2) ) &
          - ( b(2)-a(2) ) * ( d(1)-c(1) )
 
    return

Referenced by vectr_angle(), vectr_anticlockwise(), vectr_intersec(), and vectr_triangle().

Here is the caller graph for this function:

vectr_dot()

subroutine, public scale_vector::vectr_dot	(	real(rp), intent(out)	l,
		real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c,
		real(rp), dimension(3), intent(in)	d
	)

interior product of vector a->b and c->d

Definition at line 164 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(out) :: l
    real(RP), intent(in)  :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
    !---------------------------------------------------------------------------
    ! if a=c=zero-vector and b=d, result is abs|a|^2
 
    l = ( b(1)-a(1) ) * ( d(1)-c(1) ) &
      + ( b(2)-a(2) ) * ( d(2)-c(2) ) &
      + ( b(3)-a(3) ) * ( d(3)-c(3) )
 
    return

Referenced by vectr_angle(), vectr_anticlockwise(), vectr_intersec(), and vectr_triangle_plane().

Here is the caller graph for this function:

vectr_abs()

subroutine, public scale_vector::vectr_abs	(	real(rp), intent(out)	l,
		real(rp), dimension(3), intent(in)	a
	)

length of vector o->a

Definition at line 181 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(out) :: l
    real(RP), intent(in)  :: a(3) ! x,y,z(cartesian)
    !---------------------------------------------------------------------------
 
    l = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
    l = sqrt(l)
 
    return

Referenced by vectr_angle(), vectr_intersec(), and vectr_triangle().

Here is the caller graph for this function:

vectr_angle()

subroutine, public scale_vector::vectr_angle	(	real(rp), intent(out)	angle,
		real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c
	)

calc angle between two vector(b->a,b->c)

Definition at line 196 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(out) :: angle
    real(RP), intent(in)  :: a(3), b(3), c(3)
 
    real(RP) :: nv(3), nvlenS, nvlenC
    !---------------------------------------------------------------------
 
    call vectr_dot  ( nvlenc, b, a, b, c )
    call vectr_cross( nv(:),  b, a, b, c )
    call vectr_abs  ( nvlens, nv(:) )
    angle = atan2( nvlens, nvlenc )
 
    return

References vectr_abs(), vectr_cross(), and vectr_dot().

Referenced by vectr_anticlockwise(), vectr_intersec(), and vectr_triangle().

Here is the call graph for this function:

Here is the caller graph for this function:

vectr_intersec()

subroutine, public scale_vector::vectr_intersec	(	logical, intent(out)	ifcross,
		real(rp), dimension(3), intent(out)	p,
		real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c,
		real(rp), dimension(3), intent(in)	d
	)

judge intersection of two vector

Definition at line 215 of file scale_vector.F90.

    use scale_const, only: &
       eps => const_eps
    implicit none
 
    logical,  intent(out) :: ifcross
    ! .true. : line a->b and c->d intersect
    ! .false.: line a->b and c->d do not intersect and p = (0,0)
    real(RP), intent(out) :: p(3) ! intersection point
    real(RP), intent(in)  :: a(3), b(3), c(3), d(3)
 
    real(RP), parameter :: o(3) = 0.0_rp
 
    real(RP) :: oaob(3), ocod(3), cdab(3)
    real(RP) :: ip, length
    real(RP) :: angle_aop, angle_pob, angle_aob
    real(RP) :: angle_cop, angle_pod, angle_cod
    !---------------------------------------------------------------------
 
    call vectr_cross( oaob, o, a, o, b )
    call vectr_cross( ocod, o, c, o, d )
    call vectr_cross( cdab, o, ocod, o, oaob )
 
    call vectr_abs  ( length, cdab )
    call vectr_dot  ( ip, o, cdab, o, a )
 
    p(:) = cdab(:) / sign(length,ip)
!   LOG_INFO("VECTR_intersec",*), "p:", p(:)
 
    call vectr_angle( angle_aop, a, o, p )
    call vectr_angle( angle_pob, p, o, b )
    call vectr_angle( angle_aob, a, o, b )
!   LOG_INFO("VECTR_intersec",*), "angle a-p-b:", angle_aop, angle_pob, angle_aob
 
    call vectr_angle( angle_cop, c, o, p )
    call vectr_angle( angle_pod, p, o, d )
    call vectr_angle( angle_cod, c, o, d )
!   LOG_INFO("VECTR_intersec",*), "angle c-p-d:", angle_cop, angle_pod, angle_cod
 
!   LOG_INFO("VECTR_intersec",*), "judge:", angle_aob-(angle_aop+angle_pob), angle_cod-(angle_cop+angle_pod)
 
    ! --- judge intersection
    if (       abs(angle_aob-(angle_aop+angle_pob)) < eps &
         .AND. abs(angle_cod-(angle_cop+angle_pod)) < eps &
         .AND. abs(angle_aop) > eps                       &
         .AND. abs(angle_pob) > eps                       &
         .AND. abs(angle_cop) > eps                       &
         .AND. abs(angle_pod) > eps                       ) then
       ifcross = .true.
    else
       ifcross = .false.
       p(:) = 0.0_rp
    endif
 
    return

References scale_const::const_eps, vectr_abs(), vectr_angle(), vectr_cross(), and vectr_dot().

Here is the call graph for this function:

vectr_anticlockwise()

subroutine, public scale_vector::vectr_anticlockwise	(	real(rp), dimension(nvert,3), intent(inout)	vertex,
		integer, intent(in)	nvert
	)

bubble sort anticlockwise by angle

Definition at line 274 of file scale_vector.F90.

    use scale_const, only: &
       eps => const_eps
    implicit none
 
    integer,  intent(in)    :: nvert
    real(RP), intent(inout) :: vertex(nvert,3)
 
    real(RP), parameter :: o(3) = 0.0_rp
 
    real(RP) :: v1(3), v2(3), v3(3)
    real(RP) :: xp(3), ip
    real(RP) :: angle1, angle2
 
    integer :: i, j
    !---------------------------------------------------------------------
 
    do j = 2  , nvert-1
    do i = j+1, nvert
       v1(:) = vertex(1,:)
       v2(:) = vertex(j,:)
       v3(:) = vertex(i,:)
 
       call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
       call vectr_dot  ( ip, o(:), v1(:), o(:), xp(:) )
 
       if ( ip < -eps ) then ! right hand : exchange
!         LOG_INFO("VECTR_anticlockwise",*) 'exchange by ip', i, '<->',j
          vertex(i,:) = v2(:)
          vertex(j,:) = v3(:)
       endif
 
    enddo
    enddo
 
    v1(:) = vertex(1,:)
    v2(:) = vertex(2,:)
    v3(:) = vertex(3,:)
    ! if 1->2->3 is on the line
    call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
    call vectr_dot  ( ip, o(:), v1(:), o(:), xp(:) )
    call vectr_angle( angle1, v1(:), o, v2(:) )
    call vectr_angle( angle2, v1(:), o, v3(:) )
!   LOG_INFO("VECTR_anticlockwise",*) ip, angle1, angle2, abs(angle1)-abs(angle2)
 
    if (       abs(ip)                 < eps  &      ! on the same line
         .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
!      LOG_INFO("VECTR_anticlockwise",*) 'exchange by angle', 2, '<->', 3
       vertex(2,:) = v3(:)
       vertex(3,:) = v2(:)
    endif
 
    v2(:) = vertex(nvert  ,:)
    v3(:) = vertex(nvert-1,:)
    ! if 1->nvert->nvert-1 is on the line
    call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
    call vectr_dot  ( ip, o(:), v1(:), o(:), xp(:) )
    call vectr_angle( angle1, v1(:), o, v2(:) )
    call vectr_angle( angle2, v1(:), o, v3(:) )
!   LOG_INFO("VECTR_anticlockwise",*) ip, angle1, angle2, abs(angle1)-abs(angle2)
 
    if (       abs(ip)                 < eps  &      ! on the same line
         .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
!      LOG_INFO("VECTR_anticlockwise",*) 'exchange by angle', nvert, '<->', nvert-1
       vertex(nvert,  :) = v3(:)
       vertex(nvert-1,:) = v2(:)
    endif
 
    return

References scale_const::const_eps, vectr_angle(), vectr_cross(), and vectr_dot().

Here is the call graph for this function:

vectr_triangle()

real(rp) function, public scale_vector::vectr_triangle	(	real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c,
		character(len=*), intent(in)	polygon_type,
		real(rp), intent(in)	radius
	)

calc triangle area

Returns

area

Definition at line 352 of file scale_vector.F90.

    use scale_const, only: &
       pi  => const_pi, &
       eps => const_eps
    implicit none
 
    real(RP),         intent(in) :: a(3), b(3), c(3)
    character(len=*), intent(in) :: polygon_type
    real(RP),         intent(in) :: radius
    real(RP)                     :: area
 
    real(RP), parameter :: o(3) = 0.0_rp
 
    ! ON_PLANE
    real(RP) :: abc(3)
    real(RP) :: prd, r
 
    ! ON_SPHERE
    real(RP) :: angle(3)
    real(RP) :: oaob(3), oaoc(3)
    real(RP) :: oboc(3), oboa(3)
    real(RP) :: ocoa(3), ocob(3)
    real(RP) :: abab, acac
    real(RP) :: bcbc, baba
    real(RP) :: caca, cbcb
    !---------------------------------------------------------------------------
 
    area = 0.0_rp
 
    if ( polygon_type == 'ON_PLANE' ) then ! Note : On a plane, area = | ourter product of two vectors |.
 
       call vectr_cross( abc(:), a(:), b(:), a(:), c(:) )
       call vectr_abs( prd, abc(:) )
       call vectr_abs( r  , a(:)   )
 
       prd = 0.5_rp * prd !! triangle area
       if ( r < eps ) then
          print *, "zero length?", a(:)
       else
          r = 1.0_rp / r   !! 1 / length
       endif
 
       area = prd * r*r * radius*radius
 
    elseif( polygon_type == 'ON_SPHERE' ) then ! On a unit sphere, area = sum of angles - pi.
 
       ! angle 1
       call vectr_cross( oaob(:), o(:), a(:), o(:), b(:) )
       call vectr_cross( oaoc(:), o(:), a(:), o(:), c(:) )
       call vectr_abs( abab, oaob(:) )
       call vectr_abs( acac, oaoc(:) )
 
       if ( abab < eps .OR. acac < eps ) then
          !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length abab or acac:", abab, acac
          return
       endif
 
       call vectr_angle( angle(1), oaob(:), o(:), oaoc(:) )
 
       ! angle 2
       call vectr_cross( oboc(:), o(:), b(:), o(:), c(:) )
       oboa(:) = -oaob(:)
       call vectr_abs( bcbc, oboc(:) )
       baba = abab
 
       if ( bcbc < eps .OR. baba < eps ) then
          !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length bcbc or baba:", bcbc, baba
          return
       endif
 
       call vectr_angle( angle(2), oboc(:), o(:), oboa(:) )
 
       ! angle 3
       ocoa(:) = -oaoc(:)
       ocob(:) = -oboc(:)
       caca = acac
       cbcb = bcbc
 
       if ( caca < eps .OR. cbcb < eps ) then
          !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length caca or cbcb:", caca, cbcb
          return
       endif
 
       call vectr_angle( angle(3), ocoa(:), o(:), ocob(:) )
 
       ! calc area
       area = ( angle(1)+angle(2)+angle(3) - pi ) * radius*radius
 
    endif
 
    return

References scale_const::const_eps, scale_const::const_pi, vectr_abs(), vectr_angle(), and vectr_cross().

Here is the call graph for this function:

vectr_triangle_plane()

real(rp) function, public scale_vector::vectr_triangle_plane	(	real(rp), dimension(3), intent(in)	a,
		real(rp), dimension(3), intent(in)	b,
		real(rp), dimension(3), intent(in)	c
	)

calc triangle area on plane

Returns

area

Definition at line 450 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(in) :: a(3), b(3), c(3)
    real(RP)             :: area
    !
    real(RP) :: len_ab, len_ac, prd
    !---------------------------------------------------------------------------
 
    call vectr_dot( len_ab, a, b, a, b )
    call vectr_dot( len_ac, a, c, a, c )
    call vectr_dot( prd   , a, b, a, c )
 
    area = 0.5_rp * sqrt( len_ab * len_ac - prd * prd )
 

References vectr_dot().

Here is the call graph for this function:

vectr_rotation()

subroutine, public scale_vector::vectr_rotation	(	real(rp), dimension(3), intent(inout)	a,
		real(rp), intent(in)	angle,
		integer, intent(in)	iaxis
	)

Apply rotation matrix.

Definition at line 472 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(inout) :: a(3)
    real(RP), intent(in)    :: angle
    integer,  intent(in)    :: iaxis
 
    real(RP) :: m(3,3), b(3)
    !---------------------------------------------------------------------------
 
    if ( iaxis == i_xaxis ) then
       m(1,1) =        1.0_rp
       m(1,2) =        0.0_rp
       m(1,3) =        0.0_rp
 
       m(2,1) =        0.0_rp
       m(2,2) =  cos(angle)
       m(2,3) =  sin(angle)
 
       m(3,1) =        0.0_rp
       m(3,2) = -sin(angle)
       m(3,3) =  cos(angle)
    elseif( iaxis == i_yaxis ) then
       m(1,1) =  cos(angle)
       m(1,2) =        0.0_rp
       m(1,3) = -sin(angle)
 
       m(2,1) =        0.0_rp
       m(2,2) =        1.0_rp
       m(2,3) =        0.0_rp
 
       m(3,1) =  sin(angle)
       m(3,2) =        0.0_rp
       m(3,3) =  cos(angle)
    elseif( iaxis == i_zaxis ) then
       m(1,1) =  cos(angle)
       m(1,2) =  sin(angle)
       m(1,3) =        0.0_rp
 
       m(2,1) = -sin(angle)
       m(2,2) =  cos(angle)
       m(2,3) =        0.0_rp
 
       m(3,1) =        0.0_rp
       m(3,2) =        0.0_rp
       m(3,3) =        1.0_rp
    else
       return
    endif
 
    b(1) = m(1,1) * a(1) + m(1,2) * a(2) + m(1,3) * a(3)
    b(2) = m(2,1) * a(1) + m(2,2) * a(2) + m(2,3) * a(3)
    b(3) = m(3,1) * a(1) + m(3,2) * a(2) + m(3,3) * a(3)
 
    a(:) = b(:)
 
    return

References i_xaxis, i_yaxis, and i_zaxis.

vectr_distance()

subroutine, public scale_vector::vectr_distance	(	real(rp), intent(in)	r,
		real(rp), intent(in)	lon1,
		real(rp), intent(in)	lat1,
		real(rp), intent(in)	lon2,
		real(rp), intent(in)	lat2,
		real(rp), intent(out)	dist
	)

Get horizontal distance on the sphere.

Definition at line 539 of file scale_vector.F90.

    implicit none
 
    real(RP), intent(in)  :: r          ! radius in meter
    real(RP), intent(in)  :: lon1, lat1 ! in radian
    real(RP), intent(in)  :: lon2, lat2 ! in radian
    real(RP), intent(out) :: dist       ! distance of the two points in meter
 
    real(RP) :: gmm, gno_x, gno_y
    !-----------------------------------------------------------------------
 
    gmm = sin(lat1) * sin(lat2) &
        + cos(lat1) * cos(lat2) * cos(lon2-lon1)
 
    gno_x = ( cos(lat2) * sin(lon2-lon1) )
    gno_y = ( cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2-lon1) )
 
    dist = r * atan2( sqrt(gno_x*gno_x+gno_y*gno_y), gmm )
 
    return

Referenced by scale_atmos_grid_cartesc_metric::atmos_grid_cartesc_metric_rotcoef().

Here is the caller graph for this function:

Variable Documentation

i_xaxis

integer, parameter, public scale_vector::i_xaxis = 1

Definition at line 43 of file scale_vector.F90.

  integer, public, parameter :: I_Xaxis = 1

Referenced by vectr_rotation().

i_yaxis

integer, parameter, public scale_vector::i_yaxis = 2

Definition at line 44 of file scale_vector.F90.

  integer, public, parameter :: I_Yaxis = 2

Referenced by vectr_rotation().

i_zaxis

integer, parameter, public scale_vector::i_zaxis = 3

Definition at line 45 of file scale_vector.F90.

  integer, public, parameter :: I_Zaxis = 3

Referenced by vectr_rotation().