1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
*DECK LIBMPA
SUBROUTINE LIBMPA(NOR,JINI,WEIGHT,BASEPT,DEMP,SP)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute a set of partial base points preserving NOR partial moments
* of a function using the modified Ribon approach.
*
*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
* NOR number of partial moments to preserve.
* JINI minimum order of the partial moment we want to preserve.
* We must have 1-NOR <= JINI <= 0 (order 0 moment is always
* preserved).
* WEIGHT weights of the probability table.
* BASEPT base points of the probability table.
* DEMP partial moments.
*
*Parameters: output
* SP base points for the partial cross section.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NOR,JINI
DOUBLE PRECISION DEMP(JINI:NOR+JINI-1)
REAL WEIGHT(NOR),BASEPT(NOR),SP(NOR)
*----
* LOCAL VARIABLES
*----
PARAMETER (MAXNOR=20)
DOUBLE PRECISION DDA(0:MAXNOR),DD,DSIGX
*
IF(NOR.GT.MAXNOR) CALL XABORT('LIBMPA: STORAGE OVERFLOW.')
IF(NOR.LE.0) CALL XABORT('LIBMPA: NEGATIVE OR ZERO VALUE OF NOR.')
IF((1-NOR.GT.JINI).OR.(JINI.GT.0)) CALL XABORT('LIBMPA: INCONSIST'
1 //'ENT VALUE OF JINI.')
*
DO 50 I=1,NOR
DSIGX=DBLE(BASEPT(I))
DDA(0)=1.0D0
J0=0
DO 20 J=1,NOR
IF(J.EQ.I) GO TO 20
J0=J0+1
DDA(J0)=DDA(J0-1)
DO 10 K=1,J0-1
DDA(J0-K)=DDA(J0-K-1)-DDA(J0-K)*DBLE(BASEPT(J))
10 CONTINUE
DDA(0)=-DDA(0)*DBLE(BASEPT(J))
20 CONTINUE
DD=0.0D0
DO 30 J=0,NOR-1
DD=DD+DDA(J)*DEMP(J+JINI)
30 CONTINUE
DO 40 J=1,NOR
IF(J.NE.I) DD=DD/(DBLE(BASEPT(J))-DSIGX)
40 CONTINUE
SP(I)=REAL(((-1.0D0)**(NOR-1))*DD*DSIGX**(-JINI))/WEIGHT(I)
50 CONTINUE
RETURN
END
|
