Initial commit from Polytechnique Montreal

1 files changed, 181 insertions, 0 deletions

diff --git a/Donjon/src/DETCTL.f b/Donjon/src/DETCTL.f
new file mode 100644
index 0000000..111ca23
--- /dev/null
+++ b/Donjon/src/DETCTL.f
@@ -0,0 +1,181 @@
+*DECK DETCTL
+      SUBROUTINE DETCTL(NX,NY,NZ,NEL,VECT,RESP,NDET,XCT,YCT,ZCT,COR,
+     1           KEYF,IPRT)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Call the subroutines that perform the parabolic interpolation.
+*
+*Copyright:
+* Copyright (C) 2010 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* E. Varin, M. Guyot
+*
+*Parameters:
+* NX     number of x mesh-splitted elements 
+* NY     number of y mesh-splitted elements 
+* NZ     number of z mesh-splitted elements
+* NEL    number of finite elements
+* VECT
+* RESP   flux reads by the detector
+* NDET   number of detectors
+* XCT    center coordinates of each mesh-splitted elements for x
+* YCT    center coordinates of each mesh-splitted elements for y
+* ZCT    center coordinates of each mesh-splitted elements for z
+* COR    coordinates of the center of the detector
+* KEYF   keyflux recover from L_TRACK object
+* IPRT   printing index
+*
+*-----------------------------------------------------------------------                                                                    *
+ 
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER    NX,NY,NZ,NEL,NDET,KEYF(NEL),IPRT
+      REAL       VECT(*),COR(*),XCT,YCT,ZCT,RESP(NDET)
+*----
+*  LOCAL VARIABLES
+*----
+      REAL      D1,D2,D3,X1,X2,X3,XX1,XX2,XX3,Y1,Y2,
+     1          Y3,YY1,YY2,YY3,Z1,Z2,Z3,ZZ1,ZZ2,ZZ3,PD1,PD2,PD3,PPD1,
+     2          PPD2,PPD3,CE,BE,AH
+      INTEGER   I,III,NIJK,I1,I2,I3,IP1,IP2,IP3,J1,J2,J3,JP1,JP2,JP3,
+     1          K1,K2,K3,KP1,KP2,KP3,K0,JJJ
+
+      IF(IPRT.GT.4) WRITE(6,1000)
+
+      IF (NDET.LE.0) RETURN
+      NIJK = NX*NY
+
+      DO 10 III=1,NDET
+         I = (III-1)*3
+         D1 = COR(I+1)
+         D2 = COR(I+2)
+         D3 = COR(I+3)
+*----
+*  DETERMINE CENTER OF INTERPOLATE RANGE
+*----
+         CALL DETRTR(D1,XCT,NX,XX1,XX2,XX3,IP1,IP2,IP3)
+         X1 = XX1
+         X2 = XX2
+         X3 = XX3
+         I1 = IP1
+         I2 = IP2
+         I3 = IP3
+
+         CALL DETRTR(D2,YCT,NY,YY1,YY2,YY3,JP1,JP2,JP3)
+         Y1 = YY1
+         Y2 = YY2
+         Y3 = YY3
+         J1 = JP1
+         J2 = JP2
+         J3 = JP3
+
+         CALL DETRTR(D3,ZCT,NZ,ZZ1,ZZ2,ZZ3,KP1,KP2,KP3)
+         Z1 = ZZ1
+         Z2 = ZZ2
+         Z3 = ZZ3
+         K1 = KP1
+         K2 = KP2
+         K3 = KP3
+
+         IF (IPRT.GT.4) THEN
+           IF (MOD(III,25).EQ.0) WRITE(6,1000)
+         ENDIF
+
+         IF(IPRT.GT.4) THEN
+           WRITE(6,2000) III,D1,X1,X2,X3,D2,Y1,Y2,Y3,D3,Z1,Z2,Z3,
+     >                          I1,I2,I3,   J1,J2,J3,   K1,K2,K3
+         ENDIF
+*----
+*  INTERPOLATION IN X AT PLANE Z=K1
+*----
+         K0   = (K1-1)*NIJK
+*----
+*  INTERPOLATION IN X AT PLANE Y=J1,Z=K1
+*----
+         JJJ  = NX*(J1-1)
+         PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J2,Z=K1
+*----
+         JJJ  = NX*(J2-1)
+         PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J3,Z=K1
+*----
+         JJJ  = NX*(J3-1)
+         PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN Y AT PLANE Z=K1
+*----
+         CALL DETPAR(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
+         PPD1 = AH*D2*D2 + BE*D2 + CE
+*----
+*  INTERPOLATION IN X AT PLANE Z=K2
+*----
+         K0   = (K2-1)*NIJK
+*----
+*  INTERPOLATION IN X AT PLANE Y=J1,Z=K2
+*----
+         JJJ  = NX*(J1-1)
+         PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J2,Z=K2
+*----
+         JJJ  = NX*(J2-1)
+         PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J3,Z=K2
+*----
+         JJJ  = NX*(J3-1)
+         PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN Y AT PLANE Z=K2
+*----
+         CALL DETPAR2(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
+         PPD2 = AH*D2*D2 + BE*D2 + CE
+*----
+*  INTERPOLATION IN X AT PLANE Z=K3
+*----
+         K0   = (K3-1)*NIJK
+*----
+*  INTERPOLATION IN X AT PLANE Y=J1,Z=K3
+*----
+         JJJ  = NX*(J1-1)
+         PD1 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J2,Z=K3
+*----
+         JJJ  = NX*(J2-1)
+         PD2 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN X AT PLANE Y=J3,Z=K3
+*----
+         JJJ  = NX*(J3-1)
+         PD3 = DETPOL(VECT,KEYF,JJJ,K0,I1,I2,I3,X1,X2,X3,D1)
+*----
+*  INTERPOLATION IN Y AT PLANE Z=K3
+*----
+         CALL DETPAR2(Y1,Y2,Y3,PD1,PD2,PD3,AH,BE,CE)
+         PPD3 = AH*D2*D2 + BE*D2 + CE
+*----
+*  INTERPOLATION IN Z
+*----
+         CALL DETPAR2(Z1,Z2,Z3,PPD1,PPD2,PPD3,AH,BE,CE)
+         RESP(III) = AH*D3*D3 + BE*D3 + CE
+
+  10  CONTINUE
+
+      RETURN
+
+ 1000 FORMAT(//,57X,'BRACKETING PROCESS',
+     >        /,57X,'******************',
+     >       //,5X,'DET',4X,'X ',8X,'X1',8X,'X2',8X,'X3',4X,
+     >                   4X,'Y ',8X,'Y1',8X,'Y2',8X,'Y3',4X,
+     >                   4X,'Z ',8X,'Z1',8X,'Z2',8X,'Z3',4X,/)
+ 2000 FORMAT(5X,I3.3,12F10.4,/,5X,3X,3(10X,3(2X,I6.6,2X)))
+
+      END