diff options
Diffstat (limited to 'Utilib/src/ALCACT.f')
| -rw-r--r-- | Utilib/src/ALCACT.f | 164 |
1 files changed, 164 insertions, 0 deletions
diff --git a/Utilib/src/ALCACT.f b/Utilib/src/ALCACT.f new file mode 100644 index 0000000..9318176 --- /dev/null +++ b/Utilib/src/ALCACT.f @@ -0,0 +1,164 @@ +*DECK ALCACT + SUBROUTINE ALCACT(LCACT,NG,XG,WG) +* +*----------------------------------------------------------------------- +* +*Purpose: +* set the weights and base points for the different polar quadrature: +* - "Cactus" (Halsall) +* - "optimized" (Leonard and extension by Le Tellier) +* +*Copyright: +* Copyright (C) 1999 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): R. Roy and R. Le Tellier +* +*Reference: +* A. Leonard and C.T. Mc Daniel, "Optimal polar angles and weights for +* the characteristic method", Trans. Am. Nucl. Soc., 73, 172 (1995). +* +*Parameters: input +* LCACT type of quadrature (=1,2: values used by Halsall in +* Cactus (1980); =3: values proposed by Mc Daniel ng=2 only, +* extended to 3 and 4 with the same approach by R. Le Tellier +* in 06/2005) +* NG number of weights/base points. +* +*Parameters: output +* XG base points. +* WG weights. +* +*----------------------------------------------------------------------- +* + IMPLICIT NONE +*---- +* SUBROUTINE ARGUMENTS +*---- + INTEGER LCACT,NG + REAL XG(NG),WG(NG) +*---- +* LOCAL VARIABLES +*---- + INTEGER IG + DOUBLE PRECISION PI, PHIM, PHIP, DELPHI, WTEST, XTEST, YTEST, + > ZTEST, ZERO, ONE, HALF, QUART + PARAMETER( PI=3.14159265358979D0, ZERO=0.D0, ONE=1.D0, + > HALF=0.5D0, QUART=0.25D0 ) +* For "Optimized" Quadrature + INTEGER I, I2, I3, I4, I5, IDEP, IFIN + PARAMETER ( I2= 1, I3= 3, I4= 6, I5= 10) + REAL XOP0(I5-1), WOP0(I5-1) ! Quadrature which minimizes the error on Ki2 without constraints. + REAL XOP1(I5-1), WOP1(I5-1) ! Quadrature which minimizes the error on Ki2 with P1 constraints. + REAL XGOP(I5-1), WGOP(I5-1) ! Gauss optimized quadrature. + SAVE XOP0, WOP0, XOP1, WOP1, XGOP, WGOP +* NG=2 + DATA (XOP0(I),I=I2,I3-1) / 0.273658, 0.865714/ + DATA (WOP0(I),I=I2,I3-1) / 0.139473, 0.860527/ + DATA (XOP1(I),I=I2,I3-1) / 0.340183, 0.894215/ + DATA (WOP1(I),I=I2,I3-1) / 0.194406, 0.805594/ + DATA (XGOP(I),I=I2,I3-1) / 0.399374, 0.914448/ + DATA (WGOP(I),I=I2,I3-1) / 0.250547, 0.749453/ +* NG=3 + DATA (XOP0(I),I=I3,I4-1) / 0.891439, 0.395534, 0.099812/ + DATA (WOP0(I),I=I3,I4-1) / 0.793820, 0.188560, 0.017620/ + DATA (XOP1(I),I=I3,I4-1) / 0.131209, 0.478170, 0.920079/ + DATA (WOP1(I),I=I3,I4-1) / 0.029991, 0.250860, 0.719149/ + DATA (XGOP(I),I=I3,I4-1) / 0.231156, 0.639973, 0.954497/ + DATA (WGOP(I),I=I3,I4-1) / 0.085302, 0.341456, 0.573242/ +* NG=4 + DATA (XOP0(I),I=I4,I5-1) / 0.464167, 0.908274, 0.166004, 0.042181/ + DATA (WOP0(I),I=I4,I5-1) / 0.218331, 0.746430, 0.032141, 0.003098/ + DATA (XOP1(I),I=I4,I5-1) / 0.054819, 0.212313, 0.546065, 0.932318/ + DATA (WOP1(I),I=I4,I5-1) / 0.005225, 0.051270, 0.272789, 0.670716/ + DATA (XGOP(I),I=I4,I5-1) / 0.152641, 0.450820, 0.769181, 0.972320/ + DATA (WGOP(I),I=I4,I5-1) / 0.037508, 0.167623, 0.338496, 0.456373/ +! DATA (XGOP(I),I=I5,I6-1) / 0.159153, 0.450941, 0.724750, 0.868405, +! 1 0.980414/ +! DATA (WGOP(I),I=I5,I6-1) / 0.039977, 0.156595, 0.220669, 0.204037, +! 1 0.378722/ +* + IF( LCACT.EQ.1 )THEN +*--- +* CACTUS 1: +*--- +* Equal weight quadrature + DELPHI= ONE/DBLE(NG) + PHIM = ZERO + XTEST = ZERO + DO 10 IG= 1, NG + PHIP = DBLE(IG) * DELPHI + WTEST = PHIP - PHIM + WG(IG)= REAL( WTEST ) + IF(IG.EQ.NG) THEN + YTEST = PI/2.0 + ELSE + YTEST = SQRT(ONE - PHIP * PHIP) * PHIP + ASIN(PHIP) + ENDIF + ZTEST = HALF * (YTEST - XTEST) / WTEST + XG(IG)= REAL( SQRT(ONE - ZTEST*ZTEST) ) + PHIM = PHIP + XTEST = YTEST + 10 CONTINUE + ELSEIF( LCACT.EQ.2 )THEN +*--- +* CACTUS 2: +*--- +* Uniformly distributed angle quadrature + DELPHI= PI/DBLE(2*NG) + PHIM = ZERO + DO 20 IG= 1, NG + PHIP = DBLE(IG) * DELPHI + WTEST = SIN(PHIP) - SIN(PHIM) + WG(IG)= REAL( WTEST ) + XTEST = HALF * (PHIP - PHIM) + YTEST = QUART * (SIN(PHIP+PHIP) - SIN(PHIM+PHIM)) + ZTEST = (XTEST + YTEST) / WTEST + XG(IG)= REAL( SQRT(ONE - ZTEST*ZTEST) ) + PHIM= PHIP + 20 CONTINUE + ELSEIF(( LCACT.GE.3 ).AND.( LCACT.LE.5 ))THEN +*--- +* OPTIMIZED: +*--- + IDEP=0 + IFIN=0 + IF( NG.EQ.2 ) THEN + IDEP=I2 + IFIN=I3-1 + ELSE IF( NG.EQ.3 ) THEN + IDEP=I3 + IFIN=I4-1 + ELSE IF( NG.EQ.4 ) THEN + IDEP=I4 + IFIN=I5-1 + ELSE + CALL XABORT('ALCACT: LCACA=3 => NG= 2, 3 OR 4') + ENDIF + IF (LCACT.EQ.3) THEN +* Quadrature which minimizes the error on Ki2 without constraints. + DO 30 IG= IDEP, IFIN + XG(IG-IDEP+1)= REAL(SQRT(ONE - XOP0(IG)*XOP0(IG))) + WG(IG-IDEP+1)= WOP0(IG) + 30 CONTINUE + ELSEIF (LCACT.EQ.4) THEN +* Quadrature which minimizes the error on Ki2 with P1 constraints. + DO 40 IG= IDEP, IFIN + XG(IG-IDEP+1)= REAL(SQRT(ONE - XOP1(IG)*XOP1(IG))) + WG(IG-IDEP+1)= WOP1(IG) + 40 CONTINUE + ELSEIF (LCACT.EQ.5) THEN +* Gauss optimized quadrature. + DO 50 IG= IDEP, IFIN + XG(IG-IDEP+1)= REAL(SQRT(ONE - XGOP(IG)*XGOP(IG))) + WG(IG-IDEP+1)= WGOP(IG) + 50 CONTINUE + ENDIF + ELSE + CALL XABORT('ALCACT: *LCACT* IN [1,5]') + ENDIF + RETURN + END |