1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE quadrature2D(s, ndir)
   4 | 
   5 |         IMPLICIT NONE
   6 | 
   7 |         include 'dom_constants.h'
   8 | 
   9 | !       IN
  10 |         DOM_INT                     :: ndir
  11 | 
  12 | !       OUT
  13 |         DOM_REAL, DIMENSION(4,ndir) :: s
  14 | 
  15 | !       LOCAL
  16 |         DOM_INT                     :: i
  17 |         DOM_REAL                    :: a, w
  18 |         DOM_REAL                    :: zero_moment
  19 |         DOM_REAL                    :: first_moment(3)
  20 |         DOM_REAL                    :: second_moment(3)
  21 | 
  22 | !       -------------------------------------------!
  23 | !       Calculate the components of each direction !
  24 | !       -------------------------------------------!
  25 | 
  26 |         w = 4*pi / real(ndir)
  27 | 
  28 |         zero_moment   = 0.
  29 |         first_moment  = 0.
  30 |         second_moment = 0.
  31 | 
  32 |         DO i=1, ndir
  33 | 
  34 |           a = (i-1)*(2*pi/real(ndir))
  35 |           s(1,i) = 0.8165*cos(a)
  36 |           s(2,i) = 0.8165*sin(a)
  37 |           s(3,i) = 0.
  38 |           s(4,i) = w
  39 | 
  40 |           zero_moment = zero_moment + w
  41 | 
  42 |           first_moment(1) = first_moment(1) + w*s(1,i)
  43 |           first_moment(2) = first_moment(2) + w*s(2,i)
  44 |           first_moment(3) = first_moment(3) + w*s(3,i)
  45 | 
  46 |           second_moment(1) = second_moment(1) + w*s(1,i)*s(1,i)
  47 |           second_moment(2) = second_moment(2) + w*s(2,i)*s(2,i)
  48 |           second_moment(3) = second_moment(3) + w*s(3,i)*s(3,i)
  49 | 
  50 |         ENDDO
  51 | 
  52 |         WRITE(*,*) " ===================="
  53 |         WRITE(*,*) " + DOM moments:"
  54 |         WRITE(*,*)
  55 |         WRITE(*,'(4x,A,1x, F9.6)') "0th moment       : ", zero_moment
  56 |         WRITE(*,'(4x,A,1x,3F9.6)') "1st moment       : ", first_moment
  57 |         WRITE(*,'(4x,A,1x,3F9.6)') "2nd moment(trace): ", second_moment
  58 |         WRITE(*,*) " ===================="
  59 |         WRITE(*,*)
  60 | 
  61 |       END SUBROUTINE quadrature2D