1 files changed, 142 insertions, 0 deletions
diff --git a/Dragon/src/SYB4VO.f b/Dragon/src/SYB4VO.f
new file mode 100644
index 0000000..e92c3d4
--- /dev/null
+++ b/Dragon/src/SYB4VO.f
@@ -0,0 +1,142 @@
+*DECK SYB4VO
+      SUBROUTINE SYB4VO(NSECT,NR,DX,DY,RAD,VOLINT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the volumes of a Cartesian sectorized cell.
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input
+* NSECT   number of sectors (multiple of 4).
+* NR      one plus the number of tubes.
+* DX      X-oriented side.
+* DY      Y-oriented side.
+* RAD     radius of the tubes.
+*
+*Parameters: output
+* VOLINT  volumes.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER   NSECT,NR
+      REAL      DX,DY,RAD(NR-1),VOLINT(NSECT,NR)
+*
+      IF(MOD(4+NSECT,4).NE.0) CALL XABORT('SYB4VO: INVALID NSECT.')
+      DDX = DX / 2.
+      DDY = DY / 2.
+      RADDD = SQRT(DDX*DDX + DDY*DDY)
+      TETA0 = 3.14159265358979 * 2. / REAL(NSECT)
+*----
+*  ITER=1 IS FOR (0, PI/4) AND ITER=2 IS FOR (PI/2, PI/4)
+*----
+      DO 15 IR=1,NR
+      DO 10 IS=1,NSECT
+      VOLINT(IS,IR)=0.0
+   10 CONTINUE
+   15 CONTINUE
+      DO 60 ITER=1,2
+      IF(ITER.EQ.1) THEN
+         ISA = 1
+         ISB = 0
+         DDXY=DDX
+      ELSE
+         ISA = - 1
+         ISB = NSECT/4 + 1
+         DDXY=DDY
+      ENDIF
+*
+      TETAC = ACOS(DDXY / RADDD)
+      TETA2   = 0.
+      DO 20 IR = 1, NR-1
+      IF (RAD(IR) .GT. DDXY) THEN
+         IRNEXT = IR
+         GOTO 30
+      ENDIF
+   20 CONTINUE
+      IRNEXT = NR
+*----
+*  IRNEXT : NEXT RADIUS INTERCEPTING A SIDE
+*  DDY2   : LAST PROCESSED COORDINATE
+*  ISNEXT : NEXT SECTOR
+*----
+   30 DDY2   = 0.
+      ISNEXT = 1
+*----
+*  NEXT SECTOR
+*----
+   40 IS = ISNEXT
+      IF (IS .NE. 0) THEN
+         ISV   = ISA * IS + ISB
+         IR0   = IRNEXT
+         DDY1  = DDY2
+         TETA1 = TETA2
+         TETA2 = IS * TETA0
+*
+*        THE ANGLE IS LIMITED BY THE DIAGONAL.
+         IF (TETA2 .GE. (TETAC - 1.E-6)) THEN
+            ISNEXT = 0
+            TETA2  = TETAC
+            RAD2 = RADDD
+         ELSE
+            ISNEXT = IS + 1
+            RAD2 = DDXY / COS(TETA2)
+         ENDIF
+*
+*        THE NEXT RADIUS IS INTERCEPTING THE SECTOR.
+         IF (IR0 .LT. NR) THEN
+            RADIR = RAD(IR0)
+            IF (RADIR .LT. (RAD2 * (1. - 1.E-6))) THEN
+               RAD2 = RADIR
+               TETA2 = ACOS(DDXY / RAD2)
+               IRNEXT = IR0 + 1
+               ISNEXT = IS
+            ELSE IF (RADIR .LE. (RAD2 * (1. + 1.E-6))) THEN
+*              THE NEXT RADIUS IS EQUAL TO THE SECTOR.
+               IRNEXT = IR0 + 1
+            ENDIF
+         ENDIF
+*
+*        DDY2 IS THE NEXT COORDINATE AND DT IS HALF THE ANGLE
+*        INCREMENT FOR THE SECTOR.
+         DDY2 = RAD2 * SIN(TETA2)
+         DT = (TETA2 - TETA1) * 0.5
+*
+*        COMPLETE TUBES.
+         R1 = 0.
+         DO 50 IR = 1, IR0 - 1
+         R0 = R1
+         R1 = RAD(IR)
+         DR = (R1 - R0) * (R1 + R0)
+         VOLINT(ISV, IR) = VOLINT(ISV, IR) + DT * DR
+   50    CONTINUE
+*
+*        LAST SIDE-INTERCEPTED TUBE.
+         VOLMAX = DDXY * 0.5 * (DDY2 - DDY1)
+         VOLMAX = VOLMAX - DT * R1 * R1
+         VOLINT(ISV, IR0) = VOLINT(ISV, IR0) + VOLMAX
+         GOTO 40
+      ENDIF
+   60 CONTINUE
+*
+      DO 90 IR=1,NR
+      DO 70 IS=NSECT/4+1,NSECT/2
+      VOLINT(IS,IR)=VOLINT(NSECT/2-IS+1,IR)
+   70 CONTINUE
+      DO 80 IS=NSECT/2+1,NSECT
+      VOLINT(IS,IR)=VOLINT(NSECT-IS+1,IR)
+   80 CONTINUE
+   90 CONTINUE
+      RETURN
+      END