SCALE-RM
scale_vector.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
11 module scale_vector
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
19  !-----------------------------------------------------------------------------
20  implicit none
21  private
22  !-----------------------------------------------------------------------------
23  !
24  !++ Public procedure
25  !
26  public :: vectr_xyz2latlon
27  public :: vectr_latlon2xyz
28  public :: vectr_cross
29  public :: vectr_dot
30  public :: vectr_abs
31  public :: vectr_angle
32  public :: vectr_intersec
33  public :: vectr_anticlockwise
34  public :: vectr_triangle
35  public :: vectr_triangle_plane
36  public :: vectr_rotation
37  public :: vectr_distance
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Public parameters & variables
42  !
43  integer, public, parameter :: i_xaxis = 1
44  integer, public, parameter :: i_yaxis = 2
45  integer, public, parameter :: i_zaxis = 3
46 
47  !-----------------------------------------------------------------------------
48  !
49  !++ Private procedure
50  !
51  !-----------------------------------------------------------------------------
52  !
53  !++ Private parameters & variables
54  !
55  !-----------------------------------------------------------------------------
56 contains
57  !-----------------------------------------------------------------------------
58  subroutine vectr_xyz2latlon( &
59  x, &
60  y, &
61  z, &
62  lat, &
63  lon )
64  use scale_const, only: &
65  eps => const_eps
66  implicit none
67 
68  real(rp), intent(in) :: x
69  real(rp), intent(in) :: y
70  real(rp), intent(in) :: z
71  real(rp), intent(out) :: lat
72  real(rp), intent(out) :: lon
73 
74  real(rp) :: length, length_h
75  !---------------------------------------------------------------------------
76 
77  length = sqrt( x*x + y*y + z*z )
78 
79  if ( length < eps ) then ! 3D vector length is
80  lat = 0.0_rp
81  lon = 0.0_rp
82  return
83  endif
84 
85  if ( z / length >= 1.0_rp ) then ! vector is parallele to z axis.
86  lat = asin(1.0_rp)
87  lon = 0.0_rp
88  return
89  elseif( z / length <= -1.0_rp ) then ! vector is parallele to z axis.
90  lat = asin(-1.0_rp)
91  lon = 0.0_rp
92  return
93  else
94  lat = asin( z / length )
95  endif
96 
97  length_h = sqrt( x*x + y*y )
98 
99  if ( length_h < eps ) then
100  lon = 0.0_rp
101  return
102  endif
103 
104  if ( x / length_h >= 1.0_rp ) then
105  lon = acos(1.0_rp)
106  elseif( x / length_h <= -1.0_rp ) then
107  lon = acos(-1.0_rp)
108  else
109  lon = acos( x / length_h )
110  endif
111 
112  if( y < 0.0_rp ) lon = -lon
113 
114  return
115  end subroutine vectr_xyz2latlon
116 
117  !-----------------------------------------------------------------------------
118  subroutine vectr_latlon2xyz( &
119  lat, &
120  lon, &
121  x, &
122  y, &
123  z, &
124  radius )
125  implicit none
126 
127  real(rp), intent(in) :: lat
128  real(rp), intent(in) :: lon
129  real(rp), intent(out) :: x
130  real(rp), intent(out) :: y
131  real(rp), intent(out) :: z
132  real(rp), intent(in) :: radius
133  !---------------------------------------------------------------------------
134 
135  x = radius * cos(lat) * cos(lon)
136  y = radius * cos(lat) * sin(lon)
137  z = radius * sin(lat)
138 
139  return
140  end subroutine vectr_latlon2xyz
141 
142  !-----------------------------------------------------------------------------
144  subroutine vectr_cross( nv, a, b, c, d )
145  implicit none
146 
147  real(rp), intent(out) :: nv(3) ! normal vector
148  real(rp), intent(in) :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
149  !---------------------------------------------------------------------------
150 
151  nv(1) = ( b(2)-a(2) ) * ( d(3)-c(3) ) &
152  - ( b(3)-a(3) ) * ( d(2)-c(2) )
153  nv(2) = ( b(3)-a(3) ) * ( d(1)-c(1) ) &
154  - ( b(1)-a(1) ) * ( d(3)-c(3) )
155  nv(3) = ( b(1)-a(1) ) * ( d(2)-c(2) ) &
156  - ( b(2)-a(2) ) * ( d(1)-c(1) )
157 
158  return
159  end subroutine vectr_cross
160 
161  !-----------------------------------------------------------------------------
163  subroutine vectr_dot( l, a, b, c, d )
164  implicit none
165 
166  real(rp), intent(out) :: l
167  real(rp), intent(in) :: a(3), b(3), c(3), d(3) ! x,y,z(cartesian)
168  !---------------------------------------------------------------------------
169  ! if a=c=zero-vector and b=d, result is abs|a|^2
170 
171  l = ( b(1)-a(1) ) * ( d(1)-c(1) ) &
172  + ( b(2)-a(2) ) * ( d(2)-c(2) ) &
173  + ( b(3)-a(3) ) * ( d(3)-c(3) )
174 
175  return
176  end subroutine vectr_dot
177 
178  !-----------------------------------------------------------------------------
180  subroutine vectr_abs( l, a )
181  implicit none
182 
183  real(rp), intent(out) :: l
184  real(rp), intent(in) :: a(3) ! x,y,z(cartesian)
185  !---------------------------------------------------------------------------
186 
187  l = a(1)*a(1) + a(2)*a(2) + a(3)*a(3)
188  l = sqrt(l)
189 
190  return
191  end subroutine vectr_abs
192 
193  !---------------------------------------------------------------------
195  subroutine vectr_angle( angle, a, b, c )
196  implicit none
197 
198  real(rp), intent(out) :: angle
199  real(rp), intent(in) :: a(3), b(3), c(3)
200 
201  real(rp) :: nv(3), nvlens, nvlenc
202  !---------------------------------------------------------------------
203 
204  call vectr_dot ( nvlenc, b, a, b, c )
205  call vectr_cross( nv(:), b, a, b, c )
206  call vectr_abs ( nvlens, nv(:) )
207  angle = atan2( nvlens, nvlenc )
208 
209  return
210  end subroutine vectr_angle
211 
212  !-----------------------------------------------------------------------------
214  subroutine vectr_intersec( ifcross, p, a, b, c, d )
215  use scale_const, only: &
216  eps => const_eps
217  implicit none
218 
219  logical, intent(out) :: ifcross
220  ! .true. : line a->b and c->d intersect
221  ! .false.: line a->b and c->d do not intersect and p = (0,0)
222  real(rp), intent(out) :: p(3) ! intersection point
223  real(rp), intent(in) :: a(3), b(3), c(3), d(3)
224 
225  real(rp), parameter :: o(3) = 0.0_rp
226 
227  real(rp) :: oaob(3), ocod(3), cdab(3)
228  real(rp) :: ip, length
229  real(rp) :: angle_aop, angle_pob, angle_aob
230  real(rp) :: angle_cop, angle_pod, angle_cod
231  !---------------------------------------------------------------------
232 
233  call vectr_cross( oaob, o, a, o, b )
234  call vectr_cross( ocod, o, c, o, d )
235  call vectr_cross( cdab, o, ocod, o, oaob )
236 
237  call vectr_abs ( length, cdab )
238  call vectr_dot ( ip, o, cdab, o, a )
239 
240  p(:) = cdab(:) / sign(length,ip)
241 ! LOG_INFO("VECTR_intersec",*), "p:", p(:)
242 
243  call vectr_angle( angle_aop, a, o, p )
244  call vectr_angle( angle_pob, p, o, b )
245  call vectr_angle( angle_aob, a, o, b )
246 ! LOG_INFO("VECTR_intersec",*), "angle a-p-b:", angle_aop, angle_pob, angle_aob
247 
248  call vectr_angle( angle_cop, c, o, p )
249  call vectr_angle( angle_pod, p, o, d )
250  call vectr_angle( angle_cod, c, o, d )
251 ! LOG_INFO("VECTR_intersec",*), "angle c-p-d:", angle_cop, angle_pod, angle_cod
252 
253 ! LOG_INFO("VECTR_intersec",*), "judge:", angle_aob-(angle_aop+angle_pob), angle_cod-(angle_cop+angle_pod)
254 
255  ! --- judge intersection
256  if ( abs(angle_aob-(angle_aop+angle_pob)) < eps &
257  .AND. abs(angle_cod-(angle_cop+angle_pod)) < eps &
258  .AND. abs(angle_aop) > eps &
259  .AND. abs(angle_pob) > eps &
260  .AND. abs(angle_cop) > eps &
261  .AND. abs(angle_pod) > eps ) then
262  ifcross = .true.
263  else
264  ifcross = .false.
265  p(:) = 0.0_rp
266  endif
267 
268  return
269  end subroutine vectr_intersec
270 
271  !---------------------------------------------------------------------
273  subroutine vectr_anticlockwise( vertex, nvert )
274  use scale_const, only: &
275  eps => const_eps
276  implicit none
277 
278  integer, intent(in) :: nvert
279  real(rp), intent(inout) :: vertex(nvert,3)
280 
281  real(rp), parameter :: o(3) = 0.0_rp
282 
283  real(rp) :: v1(3), v2(3), v3(3)
284  real(rp) :: xp(3), ip
285  real(rp) :: angle1, angle2
286 
287  integer :: i, j
288  !---------------------------------------------------------------------
289 
290  do j = 2 , nvert-1
291  do i = j+1, nvert
292  v1(:) = vertex(1,:)
293  v2(:) = vertex(j,:)
294  v3(:) = vertex(i,:)
295 
296  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
297  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
298 
299  if ( ip < -eps ) then ! right hand : exchange
300 ! LOG_INFO("VECTR_anticlockwise",*) 'exchange by ip', i, '<->',j
301  vertex(i,:) = v2(:)
302  vertex(j,:) = v3(:)
303  endif
304 
305  enddo
306  enddo
307 
308  v1(:) = vertex(1,:)
309  v2(:) = vertex(2,:)
310  v3(:) = vertex(3,:)
311  ! if 1->2->3 is on the line
312  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
313  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
314  call vectr_angle( angle1, v1(:), o, v2(:) )
315  call vectr_angle( angle2, v1(:), o, v3(:) )
316 ! LOG_INFO("VECTR_anticlockwise",*) ip, angle1, angle2, abs(angle1)-abs(angle2)
317 
318  if ( abs(ip) < eps & ! on the same line
319  .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
320 ! LOG_INFO("VECTR_anticlockwise",*) 'exchange by angle', 2, '<->', 3
321  vertex(2,:) = v3(:)
322  vertex(3,:) = v2(:)
323  endif
324 
325  v2(:) = vertex(nvert ,:)
326  v3(:) = vertex(nvert-1,:)
327  ! if 1->nvert->nvert-1 is on the line
328  call vectr_cross( xp(:), v1(:), v2(:), v1(:), v3(:) )
329  call vectr_dot ( ip, o(:), v1(:), o(:), xp(:) )
330  call vectr_angle( angle1, v1(:), o, v2(:) )
331  call vectr_angle( angle2, v1(:), o, v3(:) )
332 ! LOG_INFO("VECTR_anticlockwise",*) ip, angle1, angle2, abs(angle1)-abs(angle2)
333 
334  if ( abs(ip) < eps & ! on the same line
335  .AND. abs(angle2)-abs(angle1) < 0.0_rp ) then ! which is far?
336 ! LOG_INFO("VECTR_anticlockwise",*) 'exchange by angle', nvert, '<->', nvert-1
337  vertex(nvert, :) = v3(:)
338  vertex(nvert-1,:) = v2(:)
339  endif
340 
341  return
342  end subroutine vectr_anticlockwise
343 
344  !-----------------------------------------------------------------------------
347  function vectr_triangle( &
348  a, b, c, &
349  polygon_type, &
350  radius ) &
351  result(area)
352  use scale_const, only: &
353  pi => const_pi, &
354  eps => const_eps
355  implicit none
356 
357  real(rp), intent(in) :: a(3), b(3), c(3)
358  character(len=*), intent(in) :: polygon_type
359  real(rp), intent(in) :: radius
360  real(rp) :: area
361 
362  real(rp), parameter :: o(3) = 0.0_rp
363 
364  ! ON_PLANE
365  real(rp) :: abc(3)
366  real(rp) :: prd, r
367 
368  ! ON_SPHERE
369  real(rp) :: angle(3)
370  real(rp) :: oaob(3), oaoc(3)
371  real(rp) :: oboc(3), oboa(3)
372  real(rp) :: ocoa(3), ocob(3)
373  real(rp) :: abab, acac
374  real(rp) :: bcbc, baba
375  real(rp) :: caca, cbcb
376  !---------------------------------------------------------------------------
377 
378  area = 0.0_rp
379 
380  if ( polygon_type == 'ON_PLANE' ) then ! Note : On a plane, area = | ourter product of two vectors |.
381 
382  call vectr_cross( abc(:), a(:), b(:), a(:), c(:) )
383  call vectr_abs( prd, abc(:) )
384  call vectr_abs( r , a(:) )
385 
386  prd = 0.5_rp * prd !! triangle area
387  if ( r < eps ) then
388  print *, "zero length?", a(:)
389  else
390  r = 1.0_rp / r !! 1 / length
391  endif
392 
393  area = prd * r*r * radius*radius
394 
395  elseif( polygon_type == 'ON_SPHERE' ) then ! On a unit sphere, area = sum of angles - pi.
396 
397  ! angle 1
398  call vectr_cross( oaob(:), o(:), a(:), o(:), b(:) )
399  call vectr_cross( oaoc(:), o(:), a(:), o(:), c(:) )
400  call vectr_abs( abab, oaob(:) )
401  call vectr_abs( acac, oaoc(:) )
402 
403  if ( abab < eps .OR. acac < eps ) then
404  !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length abab or acac:", abab, acac
405  return
406  endif
407 
408  call vectr_angle( angle(1), oaob(:), o(:), oaoc(:) )
409 
410  ! angle 2
411  call vectr_cross( oboc(:), o(:), b(:), o(:), c(:) )
412  oboa(:) = -oaob(:)
413  call vectr_abs( bcbc, oboc(:) )
414  baba = abab
415 
416  if ( bcbc < eps .OR. baba < eps ) then
417  !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length bcbc or baba:", bcbc, baba
418  return
419  endif
420 
421  call vectr_angle( angle(2), oboc(:), o(:), oboa(:) )
422 
423  ! angle 3
424  ocoa(:) = -oaoc(:)
425  ocob(:) = -oboc(:)
426  caca = acac
427  cbcb = bcbc
428 
429  if ( caca < eps .OR. cbcb < eps ) then
430  !LOG_WARN("VECTR_triangle",'(A,3(ES20.10))') "zero length caca or cbcb:", caca, cbcb
431  return
432  endif
433 
434  call vectr_angle( angle(3), ocoa(:), o(:), ocob(:) )
435 
436  ! calc area
437  area = ( angle(1)+angle(2)+angle(3) - pi ) * radius*radius
438 
439  endif
440 
441  return
442  end function vectr_triangle
443 
444  !-----------------------------------------------------------------------------
447  function vectr_triangle_plane( &
448  a, b, c ) &
449  result(area)
450  implicit none
451 
452  real(rp), intent(in) :: a(3), b(3), c(3)
453  real(rp) :: area
454  !
455  real(rp) :: len_ab, len_ac, prd
456  !---------------------------------------------------------------------------
457 
458  call vectr_dot( len_ab, a, b, a, b )
459  call vectr_dot( len_ac, a, c, a, c )
460  call vectr_dot( prd , a, b, a, c )
461 
462  area = 0.5_rp * sqrt( len_ab * len_ac - prd * prd )
463 
464  end function vectr_triangle_plane
465 
466  !-----------------------------------------------------------------------------
468  subroutine vectr_rotation( &
469  a, &
470  angle, &
471  iaxis )
472  implicit none
473 
474  real(rp), intent(inout) :: a(3)
475  real(rp), intent(in) :: angle
476  integer, intent(in) :: iaxis
477 
478  real(rp) :: m(3,3), b(3)
479  !---------------------------------------------------------------------------
480 
481  if ( iaxis == i_xaxis ) then
482  m(1,1) = 1.0_rp
483  m(1,2) = 0.0_rp
484  m(1,3) = 0.0_rp
485 
486  m(2,1) = 0.0_rp
487  m(2,2) = cos(angle)
488  m(2,3) = sin(angle)
489 
490  m(3,1) = 0.0_rp
491  m(3,2) = -sin(angle)
492  m(3,3) = cos(angle)
493  elseif( iaxis == i_yaxis ) then
494  m(1,1) = cos(angle)
495  m(1,2) = 0.0_rp
496  m(1,3) = -sin(angle)
497 
498  m(2,1) = 0.0_rp
499  m(2,2) = 1.0_rp
500  m(2,3) = 0.0_rp
501 
502  m(3,1) = sin(angle)
503  m(3,2) = 0.0_rp
504  m(3,3) = cos(angle)
505  elseif( iaxis == i_zaxis ) then
506  m(1,1) = cos(angle)
507  m(1,2) = sin(angle)
508  m(1,3) = 0.0_rp
509 
510  m(2,1) = -sin(angle)
511  m(2,2) = cos(angle)
512  m(2,3) = 0.0_rp
513 
514  m(3,1) = 0.0_rp
515  m(3,2) = 0.0_rp
516  m(3,3) = 1.0_rp
517  else
518  return
519  endif
520 
521  b(1) = m(1,1) * a(1) + m(1,2) * a(2) + m(1,3) * a(3)
522  b(2) = m(2,1) * a(1) + m(2,2) * a(2) + m(2,3) * a(3)
523  b(3) = m(3,1) * a(1) + m(3,2) * a(2) + m(3,3) * a(3)
524 
525  a(:) = b(:)
526 
527  return
528  end subroutine vectr_rotation
529 
530  !-----------------------------------------------------------------------
532  subroutine vectr_distance( &
533  r, &
534  lon1, &
535  lat1, &
536  lon2, &
537  lat2, &
538  dist )
539  implicit none
540 
541  real(rp), intent(in) :: r ! radius in meter
542  real(rp), intent(in) :: lon1, lat1 ! in radian
543  real(rp), intent(in) :: lon2, lat2 ! in radian
544  real(rp), intent(out) :: dist ! distance of the two points in meter
545 
546  real(rp) :: gmm, gno_x, gno_y
547  !-----------------------------------------------------------------------
548 
549  gmm = sin(lat1) * sin(lat2) &
550  + cos(lat1) * cos(lat2) * cos(lon2-lon1)
551 
552  gno_x = ( cos(lat2) * sin(lon2-lon1) )
553  gno_y = ( cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(lon2-lon1) )
554 
555  dist = r * atan2( sqrt(gno_x*gno_x+gno_y*gno_y), gmm )
556 
557  return
558  end subroutine vectr_distance
559 
560 end module scale_vector
scale_vector::i_zaxis
integer, parameter, public i_zaxis
Definition: scale_vector.F90:45
scale_vector::vectr_triangle_plane
real(rp) function, public vectr_triangle_plane(a, b, c)
calc triangle area on plane
Definition: scale_vector.F90:450
scale_vector::vectr_anticlockwise
subroutine, public vectr_anticlockwise(vertex, nvert)
bubble sort anticlockwise by angle
Definition: scale_vector.F90:274
scale_vector::i_xaxis
integer, parameter, public i_xaxis
Definition: scale_vector.F90:43
scale_vector::vectr_latlon2xyz
subroutine, public vectr_latlon2xyz(lat, lon, x, y, z, radius)
Definition: scale_vector.F90:125
scale_vector::vectr_dot
subroutine, public vectr_dot(l, a, b, c, d)
interior product of vector a->b and c->d
Definition: scale_vector.F90:164
scale_vector::vectr_abs
subroutine, public vectr_abs(l, a)
length of vector o->a
Definition: scale_vector.F90:181
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_vector::vectr_rotation
subroutine, public vectr_rotation(a, angle, iaxis)
Apply rotation matrix.
Definition: scale_vector.F90:472
scale_vector::vectr_triangle
real(rp) function, public vectr_triangle(a, b, c, polygon_type, radius)
calc triangle area
Definition: scale_vector.F90:352
scale_vector::i_yaxis
integer, parameter, public i_yaxis
Definition: scale_vector.F90:44
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_vector::vectr_distance
subroutine, public vectr_distance(r, lon1, lat1, lon2, lat2, dist)
Get horizontal distance on the sphere.
Definition: scale_vector.F90:539
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_vector::vectr_cross
subroutine, public vectr_cross(nv, a, b, c, d)
exterior product of vector a->b and c->d
Definition: scale_vector.F90:145
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_vector::vectr_xyz2latlon
subroutine, public vectr_xyz2latlon(x, y, z, lat, lon)
Definition: scale_vector.F90:64
scale_vector
module vector
Definition: scale_vector.F90:11
scale_vector::vectr_angle
subroutine, public vectr_angle(angle, a, b, c)
calc angle between two vector(b->a,b->c)
Definition: scale_vector.F90:196
scale_vector::vectr_intersec
subroutine, public vectr_intersec(ifcross, p, a, b, c, d)
judge intersection of two vector
Definition: scale_vector.F90:215
Generated by   doxygen 1.8.17