From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 Dragon/src/NXTIAA.f | 188 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 188 insertions(+)
 create mode 100644 Dragon/src/NXTIAA.f

(limited to 'Dragon/src/NXTIAA.f')

diff --git a/Dragon/src/NXTIAA.f b/Dragon/src/NXTIAA.f
new file mode 100644
index 0000000..a962ece
--- /dev/null
+++ b/Dragon/src/NXTIAA.f
@@ -0,0 +1,188 @@
+*DECK NXTIAA
+      FUNCTION NXTIAA(POSANN ,POSPIN,VOLINT)
+*
+*----------
+*
+*Purpose:
+* Compute the volume of intersection between
+* a 2--D annular region and an annular pin.
+*
+*Copyright:
+* Copyright (C) 2005 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): G. Marleau.
+*
+*Parameters: input
+* POSANN  spatial description of the annular region with
+*         POSANN(0) the radius, POSANN(1) the $X$ position
+*         of center and POSANN(2) the $Y$ position
+*         of center.
+* POSPIN  spatial description of the annular pin region with
+*         POSPIN(0) the radius, POSPIN(1) the $X$ position
+*         of center and POSPIN(2) the $Y$ position
+*         of center.
+*
+*Parameters: output
+* NXTIAA  type of intersection between annular region and
+*         annular pin, where:
+*         = 0 means that there is no intersection
+*         between the two regions;
+*         = 1 means that the annular region
+*         is all located inside the annular pin;
+*         = 2 means that the annular pin
+*         is all located inside the annular region;
+*         =-1 means that the intersection between
+*         the annular region and the annular pin is partial.
+* VOLINT  2-D volume of intersection (area) between annular region and
+*         annular pin.
+*
+*Reference:
+*  G. Marleau,
+*  New Geometries Processing in DRAGON: The NXT: Module,
+*  Report IGE-260, Polytechnique Montreal,
+*  Montreal, 2005.
+*
+*----
+*
+      IMPLICIT         NONE
+*----
+*  Subroutine arguments
+*----
+      INTEGER          NXTIAA
+      DOUBLE PRECISION POSANN(0:2),POSPIN(0:2)
+      DOUBLE PRECISION VOLINT
+*----
+*  Local parameters
+*----
+      INTEGER          IOUT
+      CHARACTER        NAMSBR*6
+      PARAMETER       (IOUT=6,NAMSBR='NXTIAA')
+      INTEGER          IPRINT
+      PARAMETER       (IPRINT=100)
+      DOUBLE PRECISION DCUTOF
+      PARAMETER       (DCUTOF=1.0D-8)
+      DOUBLE PRECISION DZERO
+      PARAMETER       (DZERO=0.0D0)
+*----
+*  Functions
+*----
+      DOUBLE PRECISION XDRCST,PI
+*----
+*  Local variables
+*----
+      INTEGER          IFACE
+      DOUBLE PRECISION PX,PY,PX2,PY2,RA2,RP2,DAP2,DAP,XINT,YINT,
+     >                 ANGLEP,ANGLEA,SP,SA,VOLANN,VOLPIN,ACARG
+      DOUBLE PRECISION DT1,DT2,DT3
+*----
+*  Initialize NXTIAA and VOLINT
+*----
+      IF(IPRINT .GE. 200) THEN
+        WRITE(IOUT,6000) NAMSBR
+        WRITE(IOUT,6010) (POSANN(IFACE),IFACE=0,2)
+        WRITE(IOUT,6011) (POSPIN(IFACE),IFACE=0,2)
+      ENDIF
+      PI=XDRCST('Pi',' ')
+      NXTIAA=0
+      VOLINT=DZERO
+*----
+*  Find distance from center of annulus to pin center
+*  $d=\sqrt(y^{2}+y^{2})$
+*  annular region/pin intersection if $d<r_{p}+r_{a}$
+*  where $r_{a}$ is annular region radius
+*  and $r_{p}$ is pin radius
+*  Annular region in pin if $r_{a}<r_{p}$ $d<r_{p}-r_{a}$
+*  Pin in annular region if $r_{p}<r_{a}$ $d<r_{a}-r_{p}$
+*----
+      PX=POSANN(1)-POSPIN(1)
+      PY=POSANN(2)-POSPIN(2)
+      PX2=PX*PX
+      PY2=PY*PY
+      RA2=POSANN(0)*POSANN(0)
+      RP2=POSPIN(0)*POSPIN(0)
+      VOLANN=PI*RA2
+      VOLPIN=PI*RP2
+      DAP2=PX2+PY2
+      DAP=SQRT(DAP2)
+      IF(DAP .LE. ABS(POSANN(0)-POSPIN(0)) ) THEN
+        IF(POSANN(0) .LE. POSPIN(0)) THEN
+          VOLINT=VOLANN
+          NXTIAA=1
+        ELSE
+          VOLINT=VOLPIN
+          NXTIAA=2
+        ENDIF
+      ELSE IF(DAP .LT. POSANN(0)+POSPIN(0)) THEN
+*----
+*  Find intersection points $(x_{i},y_{i})$ assuming
+*  the annular region is rotated in such a way that its center
+*  is located at $(d,0)$.
+*  Equation for annular region is :
+*    $(x-d)^{2}+y^{2}=r^{2}_{a}$
+*  Equation for pin region is :
+*    $x^{2}+y^{2}=r^{2}_{p}$
+*  Intersection point is :
+*    $x_{i}=(d^{2}+r^{2}_{p}-r^{2}_{s})/(2*d)$
+*    $y_{i}=\sqrt{r^{2}_{p}-x^{2}_{i}}$
+*----
+        XINT=(DAP2+RP2-RA2)/(2*DAP)
+        YINT=SQRT(RP2-XINT*XINT)
+*----
+*  Find angular opening for pin and annulus
+*  Pin surface to the right of $x=x_{i}$ plane
+*  $S_{p}=r_{p}^{2}\arccos(x_{i}/r_{p})-x_{i}y_{i}$
+*  Annular surface to the left of $x=x_{i}$ plane
+*  $S_{a}=r_{a}^{2}\arccos((d-x_{i})/r_{a})-(d-x_{i})y_{i}$
+*----
+        ACARG=XINT/POSPIN(0)
+        IF(ACARG .GE. 1.0D0) THEN
+          ANGLEP=ACOS(1.0D0)
+        ELSE IF(ACARG .LE. -1.0D0) THEN
+          ANGLEP=ACOS(-1.0D0)
+        ELSE
+          ANGLEP=ACOS(ACARG)
+        ENDIF
+        ACARG=(DAP-XINT)/POSANN(0)
+        IF(ACARG .GE. 1.0D0) THEN
+          ANGLEA=ACOS(1.0D0)
+        ELSE IF(ACARG .LE. -1.0D0) THEN
+          ANGLEA=ACOS(-1.0D0)
+        ELSE
+          ANGLEA=ACOS(ACARG)
+        ENDIF
+        SP=RP2*ANGLEP-XINT*YINT
+        SA=RA2*ANGLEA-(DAP-XINT)*YINT
+        VOLINT=SP+SA
+        DT1=ABS(VOLINT-VOLPIN)
+        DT2=ABS(VOLINT-VOLANN)
+        DT3=ABS(VOLINT)
+        IF(DT1 .LT. DCUTOF) THEN
+          VOLINT=VOLPIN
+          NXTIAA=2
+        ELSE IF(DT2 .LT. DCUTOF) THEN
+          VOLINT=VOLANN
+          NXTIAA=1
+        ELSE IF(DT3 .GT. DCUTOF) THEN
+          NXTIAA=-1
+        ELSE
+          VOLINT=DZERO
+        ENDIF
+      ENDIF
+      IF(IPRINT .GE. 200) THEN
+        WRITE(IOUT,6012) NAMSBR,NXTIAA,VOLINT
+        WRITE(IOUT,6001) NAMSBR
+      ENDIF
+      RETURN
+*----
+*  Output formats
+*----
+ 6000 FORMAT('(* Output from --',A6,'-- follows ')
+ 6001 FORMAT('   Output from --',A6,'-- completed *)')
+ 6010 FORMAT('POSANN={',2(F20.10,','),F20.10,'};')
+ 6011 FORMAT('POSPIN={',2(F20.10,','),F20.10,'};')
+ 6012 FORMAT(A6,'={',I5,',',F20.10,'};')
+      END
-- 
cgit v1.2.3