Initial commit from Polytechnique Montreal

1 files changed, 381 insertions, 0 deletions

diff --git a/Dragon/src/LIBPTT.f b/Dragon/src/LIBPTT.f
new file mode 100644
index 0000000..689a4f7
--- /dev/null
+++ b/Dragon/src/LIBPTT.f
@@ -0,0 +1,381 @@
+*DECK LIBPTT
+      SUBROUTINE LIBPTT(IGRP,NDIL,NPART,DILUT,XSDIL,GOLD,HNAMIS,IMPX,
+     1 NOR,WEIGH,SIGX,SIGP)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Transform dilution dependent total and partial self-shielded cross
+* section into probability tables.
+*
+*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
+* IGRP    energy group index of the probability table.
+* NDIL    number of finite dilutions.
+* NPART   number of partial cross sections.
+* DILUT   dilutions with DILUT(NDIL+1)=1.e10.
+* XSDIL   dilution dependent self-shielded cross sections:
+*         XSDIL(I,1) total self-shielded cross sections;
+*         XSDIL(I,2) nu*fission self-shielded cross sections;
+*         XSDIL(I,3) P0 scattering cross sections;
+*         etc.
+*         XSDIL(NDIL+1,j) are the infinite dilution values.
+* GOLD    Goldstein-Cohen parameter.
+* HNAMIS  local name of the isotope:
+*         HNAMIS(1:8)  is the local isotope name;
+*         HNAMIS(9:12) is a suffix function of the mixture index.
+* IMPX    print parameter (equal to zero for no print).
+*
+*Parameters: output
+* NOR     order for the probability table.
+* WEIGH   quadrature weights for the probability table.
+* SIGX    base points for the total cross sections.
+* SIGP    base points for the partial cross sections.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      PARAMETER (MAXNOR=12)
+      INTEGER IGRP,NDIL,NPART,IMPX,NOR
+      REAL DILUT(NDIL+1),XSDIL(NDIL+1,NPART+1),GOLD,WEIGH(NOR),
+     1 SIGX(NOR),SIGP(MAXNOR,NPART)
+      CHARACTER HNAMIS*12
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER (TARGET=0.2D-3,EPSRID=5.0D-4)
+      REAL PRECA
+      DOUBLE PRECISION SIGXI,CC,DD,EE,DENOM
+      DOUBLE PRECISION DA(0:MAXNOR-1),DB(0:MAXNOR-1),DC(0:MAXNOR)
+      COMPLEX*16 SIGX0(MAXNOR),CCC,DCC,XCC
+      LOGICAL LCONV,LFAIL
+      CHARACTER HSMG*131
+      REAL, ALLOCATABLE, DIMENSION(:) :: WABS
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: TEST,SDDK
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: TOFIT
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(TEST(NPART+2),SDDK(NDIL),TOFIT(NDIL,2))
+*
+      IF(IMPX.GT.5) THEN
+         WRITE(6,'(/47H LIBPTT: DILUTION DEPENDANT CROSS SECTIONS OF I,
+     1   8HSOTOPE '',A12,10H'' IN GROUP,I4,1H:/9X,10HDILUTIONS=,1P,
+     2   9E12.4,:/(19X,9E12.4))') HNAMIS,IGRP,(DILUT(I),I=1,NDIL+1)
+         WRITE(6,'(10X,9HTOTAL XS=,1P,9E12.4,:/(19X,9E12.4))')
+     1   (XSDIL(I,1),I=1,NDIL+1)
+         IF(GOLD.NE.1.0) THEN
+            WRITE(6,'(6X,13HPRINCIPAL XS=,1P,9E12.4,:/(19X,9E12.4))')
+     1      (XSDIL(I,1)-(1.0-GOLD)*XSDIL(I,3),I=1,NDIL+1)
+         ENDIF
+         IF(IMPX.GT.6) THEN
+            DO 10 J=1,NPART
+            WRITE(6,'(4X,10HPARTIAL XS,I4,1H=,1P,9E12.4,:/(19X,
+     1      9E12.4))') J,(XSDIL(I,J+1),I=1,NDIL+1)
+   10       CONTINUE
+         ENDIF
+      ENDIF
+      IF(NPART.LT.2) CALL XABORT('LIBPTT: SCATTERING INFO MISSING.')
+*----
+*  CHECK IF THE ENERGY GROUP IS REALLY RESONANT.
+*----
+      IF(NDIL.EQ.0) THEN
+         NOR=1
+         WEIGH(1)=1.0
+         SIGX(1)=XSDIL(1,1)
+         DO 20 J=1,NPART
+         SIGP(1,J)=XSDIL(1,J+1)
+   20    CONTINUE
+         PREC0=0.0D0
+         PREC=0.0D0
+         PREC1=0.0D0
+         GO TO 480
+      ENDIF
+      PREC=0.0D0
+      PREC0=0.0D0
+      PREC1=0.0D0
+      LCONV=.FALSE.
+      DO 30 IDIL=1,NDIL
+      ERR=ABS(DBLE((XSDIL(NDIL+1,1))/(XSDIL(IDIL,1)))-1.0D0)
+      PREC=MAX(PREC,ERR)
+      ERR=ABS(DBLE((XSDIL(NDIL+1,1)-(1.0D0-GOLD)*XSDIL(NDIL+1,3))/
+     1 (XSDIL(IDIL,1)-(1.0-GOLD)*XSDIL(IDIL,3)))-1.0D0)
+      PREC0=MAX(PREC0,ERR)
+      ERR=ABS(DBLE((XSDIL(NDIL+1,1)-XSDIL(NDIL+1,3))/
+     1 (XSDIL(IDIL,1)-XSDIL(IDIL,3)))-1.0D0)
+      PREC1=MAX(PREC1,ERR)
+      LCONV=LCONV.OR.(ABS(XSDIL(IDIL,1)).EQ.ABS(XSDIL(IDIL+1,1)))
+   30 CONTINUE
+      IF(IMPX.GT.3) WRITE(6,'(/36H LIBPTT: ORDER  1 PROBABILITY TABLE ,
+     1 24HCALCULATION OF ISOTOPE '',A12,10H'' IN GROUP,I4,8H. ERROR=,1P,
+     2 3D11.3,1H.)') HNAMIS,IGRP,PREC0,PREC,PREC1
+      IF(PREC.LE.TARGET) THEN
+         NOR=1
+         WEIGH(1)=1.0
+         SIGX(1)=XSDIL(NDIL+1,1)
+         DO 40 J=1,NPART
+         SIGP(1,J)=XSDIL(NDIL+1,J+1)
+   40    CONTINUE
+         GO TO 360
+      ENDIF
+      IF(LCONV) THEN
+         WRITE(HSMG,'(45HLIBPTT: UNIFORM TOTAL XS IS NOT EXPECTED IN G,
+     1   4HROUP,I4,1H.)') IGRP
+         CALL XABORT(HSMG)
+      ENDIF
+*----
+*  FIND THE PADE APPROXIMATION FOR ABS+XGOLD*SIGS CROSS SECTION USING
+*  A PADE REGRESSION.
+*----
+      XGOLD=GOLD
+   45 ALLOCATE(WABS(NDIL+1))
+      DO 60 IDIL=1,NDIL+1
+      WABS(IDIL)=REAL(XSDIL(IDIL,1)-(1.0D0-XGOLD)*XSDIL(IDIL,3))
+   60 CONTINUE
+      CALL ALPLSF(3,NDIL+1,DILUT,WABS,EPSRID,.TRUE.,NOR,DA,DB,PRECA)
+      DEALLOCATE(WABS)
+      NOR=NOR+1
+      IF(NOR.GT.MAXNOR) CALL XABORT('LIBPTT: NOR IS TOO LARGE.')
+      PREC0=DBLE(PRECA)
+*
+*     FIND THE BASE POINTS IN ABS+XGOLD*SIGS CROSS SECTION.
+      SGN=1.0D0
+      DC(0)=DA(0)
+      DO 70 I=2,NOR
+      SGN=-SGN
+      DC(I-1)=SGN*(DB(I-2)+DA(I-1))
+   70 CONTINUE
+      DC(NOR)=-SGN
+      CALL ALROOT(DC,NOR,SIGX0,LFAIL)
+      IF(LFAIL) CALL XABORT('LIBPTT: POLYNOMIAL ROOT FINDING FAILURE.')
+*
+      DO 110 I=1,NOR
+*
+*     NEWTON IMPROVEMENT OF THE ROOTS.
+      CCC=0.0D0
+      XCC=1.0D0
+      DO 80 J=0,NOR
+      CCC=CCC+DC(J)*XCC
+      XCC=XCC*SIGX0(I)
+   80 CONTINUE
+      DCC=0.0D0
+      XCC=1.0D0
+      DO 85 J=1,NOR
+      DCC=DCC+DC(J)*XCC*REAL(J)
+      XCC=XCC*SIGX0(I)
+   85 CONTINUE
+      SIGX0(I)=SIGX0(I)-CCC/DCC
+*
+*     COMPUTE THE WEIGHTS.
+      IF(AIMAG(CMPLX(SIGX0(I))).NE.0.0) CALL XABORT('LIBPTT: COMPLEX '
+     1 //'ROOT.')
+      SIGXI=DBLE(SIGX0(I))
+      CC=1.0D0
+      DD=0.0D0
+      DO 90 J=0,NOR-1
+      DD=DD+DB(J)*CC
+      CC=-CC*SIGXI
+   90 CONTINUE
+      DO 100 J=1,NOR
+      IF(J.NE.I) DD=DD/(DBLE(SIGX0(J))-SIGXI)
+  100 CONTINUE
+      WEIGH(I)=REAL(DD)
+  110 CONTINUE
+*----
+*  PROCESS THE TOTAL CROSS SECTIONS.
+*----
+      DO 210 IDIL=1,NDIL
+      SCC=DA(NOR-1)
+      DO 200 I=NOR-2,0,-1
+      SCC=DA(I)+SCC*DILUT(IDIL)
+  200 CONTINUE
+      SDDK(IDIL)=(XSDIL(IDIL,1)-(1.0-XGOLD)*XSDIL(IDIL,3))/SCC
+      TOFIT(IDIL,1)=DILUT(IDIL)
+      TOFIT(IDIL,2)=XSDIL(IDIL,1)/SDDK(IDIL)
+      SDDK(IDIL)=SDDK(IDIL)*SDDK(IDIL)
+  210 CONTINUE
+      IF(XGOLD.NE.1.0) THEN
+         CALL ALDFIT(NDIL,NOR-1,TOFIT(1,1),TOFIT(1,2),SDDK,DA)
+      ENDIF
+      DO 220 I=0,NOR-1
+      DA(I)=DA(I)*XSDIL(NDIL+1,1)/DA(NOR-1)
+  220 CONTINUE
+*----
+*  COMPUTE THE BASE POINTS IN TOTAL CROSS SECTION.
+*----
+      DO 240 I=1,NOR
+      SIGXI=DBLE(SIGX0(I))
+      CC=1.0D0
+      DD=0.0D0
+      EE=0.0D0
+      DO 230 J=0,NOR-1
+      DD=DD+DA(J)*CC
+      EE=EE+DB(J)*CC
+      CC=-CC*SIGXI
+  230 CONTINUE
+      SIGX(I)=REAL(DD/EE)
+      IF(SIGX(I).LT.0.0) THEN
+         IF(XGOLD.EQ.1.0) CALL XABORT('LIBPTT: NEGATIVE BASE POINTS FO'
+     1   //'R THE TOTAL CROSS SECTION.')
+         XGOLD=MIN(1.0D0,XGOLD+0.1D0)
+         GO TO 45
+      ENDIF
+  240 CONTINUE
+*----
+*  PROCESS THE PARTIAL CROSS SECTIONS.
+*----
+      DO 300 IPART=1,NPART
+      IF(XSDIL(NDIL+1,IPART+1).EQ.0.0) THEN
+         DO 250 I=1,NOR
+         SIGP(I,IPART)=0.0
+  250    CONTINUE
+         GO TO 300
+      ENDIF
+      DO 260 IDIL=1,NDIL
+      TOFIT(IDIL,1)=DILUT(IDIL)
+      TOFIT(IDIL,2)=XSDIL(IDIL,IPART+1)/SQRT(SDDK(IDIL))
+  260 CONTINUE
+      CALL ALDFIT(NDIL,NOR-1,TOFIT(1,1),TOFIT(1,2),SDDK,DA)
+      IF(DA(NOR-1).EQ.0.0) THEN
+         DO 265 I=1,NOR
+         SIGP(I,IPART)=XSDIL(NDIL+1,IPART+1)
+  265    CONTINUE
+         GO TO 300
+      ENDIF
+      DO 270 I=0,NOR-1
+      DA(I)=DA(I)*XSDIL(NDIL+1,IPART+1)/DA(NOR-1)
+  270 CONTINUE
+*----
+*  COMPUTE THE BASE POINTS IN PARTIAL CROSS SECTION.
+*----
+      DO 290 I=1,NOR
+      SIGXI=DBLE(SIGX0(I))
+      CC=1.0D0
+      DD=0.0D0
+      EE=0.0D0
+      DO 280 J=0,NOR-1
+      DD=DD+DA(J)*CC
+      EE=EE+DB(J)*CC
+      CC=-CC*SIGXI
+  280 CONTINUE
+      SIGP(I,IPART)=REAL(DD/EE)
+  290 CONTINUE
+  300 CONTINUE
+*----
+*  REMOVING SMALL PROBABILITIES.
+*----
+      INOR=0
+  330 INOR=INOR+1
+      IF(INOR.GT.NOR) GO TO 360
+      IF(ABS(WEIGH(INOR)).LE.5.0E-7) THEN
+         DO 355 JNOR=INOR+1,NOR
+         WEIGH(JNOR-1)=WEIGH(JNOR)
+         SIGX(JNOR-1)=SIGX(JNOR)
+         DO 350 J=1,NPART
+         SIGP(JNOR-1,J)=SIGP(JNOR,J)
+  350    CONTINUE
+  355    CONTINUE
+         INOR=INOR-1
+         NOR=NOR-1
+      ENDIF
+      GO TO 330
+*----
+*  NORMALIZE THE PROBABILITY TABLE TO INFINITE DILUTION X-S.
+*----
+  360 CC=0.0D0
+      DO 390 I=1,NOR
+      CC=CC+WEIGH(I)
+  390 CONTINUE
+      DO 400 I=1,NOR
+      WEIGH(I)=WEIGH(I)/REAL(CC)
+  400 CONTINUE
+      CC=0.0D0
+      DO 410 I=1,NOR
+      CC=CC+WEIGH(I)*SIGX(I)
+  410 CONTINUE
+      IF(CC.NE.0.0) THEN
+         DO 420 I=1,NOR
+         SIGX(I)=SIGX(I)*XSDIL(NDIL+1,1)/REAL(CC)
+  420    CONTINUE
+      ENDIF
+      DO 450 J=1,NPART
+      CC=0.0D0
+      DO 430 I=1,NOR
+      CC=CC+WEIGH(I)*SIGP(I,J)
+  430 CONTINUE
+      IF(CC.NE.0.0) THEN
+         DO 440 I=1,NOR
+         SIGP(I,J)=SIGP(I,J)*XSDIL(NDIL+1,J+1)/REAL(CC)
+  440    CONTINUE
+      ENDIF
+  450 CONTINUE
+*----
+*  TEST THE ACCURACY OF THE PROBABILITY TABLE.
+*----
+      PREC=0.0D0
+      PREC1=0.0D0
+      DO 470 IDIL=1,NDIL+1
+         CC=0.0D0
+         DD=0.0D0
+         EE=0.0D0
+         DO 460 I=1,NOR
+         DENOM=SIGX(I)-(1.0-GOLD)*SIGP(I,2)+DILUT(IDIL)
+         CC=CC+WEIGH(I)/DENOM
+         DD=DD+WEIGH(I)*SIGX(I)/DENOM
+         EE=EE+WEIGH(I)*(SIGX(I)-SIGP(I,2))/DENOM
+  460    CONTINUE
+         PREC=MAX(PREC,ABS((DD/CC)/DBLE(XSDIL(IDIL,1))-1.0D0))
+         PREC1=MAX(PREC1,ABS((EE/CC)/DBLE(XSDIL(IDIL,1)-
+     1   XSDIL(IDIL,3))-1.0D0))
+  470 CONTINUE
+  480 IF((IMPX.GE.3).AND.(NOR.GT.1)) THEN
+         WRITE(6,'(14H LIBPTT: ORDER,I3,27H PROBABILITY TABLE CALCULAT,
+     1   16HION OF ISOTOPE '',A12,10H'' IN GROUP,I4,8H. ERROR=,1P,
+     2   3D11.3,1H.)') NOR,HNAMIS,IGRP,PREC0,PREC,PREC1
+      ENDIF
+      IF((IMPX.GE.2).AND.(NOR.GT.1)) THEN
+         WRITE(6,'(/34H LIBPTT: SIGT BASE POINTS IN GROUP,I4,1H:)') IGRP
+         WRITE(6,'(1X,1P,12E12.4)') SIGX(:NOR)
+      ENDIF
+*
+      IF(((IMPX.GT.4).AND.(NOR.GT.1)).OR.(IMPX.GT.5)) THEN
+         WRITE(6,'(/27H LIBPTT: PROBABILITY TABLE:/7X,11HPROBABILITY,
+     1   2X,10HTOTAL-----,2X,22HPARTIAL CROSS SECTIONS)')
+         DO 490 J=1,NPART+2
+         TEST(J)=0.0D0
+  490    CONTINUE
+         DO 510 INOR=1,NOR
+         TEST(1)=TEST(1)+WEIGH(INOR)
+         TEST(2)=TEST(2)+WEIGH(INOR)*SIGX(INOR)
+         DO 500 J=1,NPART
+         TEST(J+2)=TEST(J+2)+WEIGH(INOR)*SIGP(INOR,J)
+  500    CONTINUE
+         WRITE(6,'(1X,I5,1P,7E12.4,:/(30X,5E12.4))') INOR,
+     1   WEIGH(INOR),SIGX(INOR),(SIGP(INOR,J),J=1,NPART)
+  510    CONTINUE
+         WRITE(6,'(6H CHECK,1P,7E12.4,:/(30X,5E12.4))') (REAL(TEST(J)),
+     1   J=1,NPART+2)
+         TEST(1)=1.0D0
+         DO 520 J=1,NPART+1
+         TEST(J+1)=XSDIL(NDIL+1,J)
+  520    CONTINUE
+         WRITE(6,'(6H EXACT,1P,7E12.4,:/(30X,5E12.4))') (REAL(TEST(J)),
+     1   J=1,NPART+2)
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(TOFIT,SDDK,TEST)
+      RETURN
+      END