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
|
C Copyright 1981-2007 ECMWF
C
C Licensed under the GNU Lesser General Public License which
C incorporates the terms and conditions of version 3 of the GNU
C General Public License.
C See LICENSE and gpl-3.0.txt for details.
C
SUBROUTINE SETGPFA(TRIGS,N)
*
DIMENSION TRIGS(*)
DIMENSION NJ(3)
*
* DECOMPOSE N INTO FACTORS 2,3,5
* ------------------------------
NN = N
IFAC = 2
*
DO 30 LL = 1 , 3
KK = 0
10 CONTINUE
IF (MOD(NN,IFAC).NE.0) GO TO 20
KK = KK + 1
NN = NN / IFAC
GO TO 10
20 CONTINUE
NJ(LL) = KK
IFAC = IFAC + LL
30 CONTINUE
*
IF (NN.NE.1) THEN
WRITE(6,40) N
40 FORMAT(' *** WARNING!!!',I10,' IS NOT A LEGAL VALUE OF N ***')
RETURN
ENDIF
*
IP = NJ(1)
IQ = NJ(2)
IR = NJ(3)
*
* COMPUTE LIST OF ROTATED TWIDDLE FACTORS
* ---------------------------------------
NJ(1) = 2**IP
NJ(2) = 3**IQ
NJ(3) = 5**IR
*
TWOPI = 4.0 * ASIN(1.0)
I = 1
*
DO 60 LL = 1 , 3
NI = NJ(LL)
IF (NI.EQ.1) GO TO 60
*
DEL = TWOPI / FLOAT(NI)
IROT = N / NI
KINK = MOD(IROT,NI)
KK = 0
*
DO 50 K = 1 , NI
ANGLE = FLOAT(KK) * DEL
TRIGS(I) = COS(ANGLE)
TRIGS(I+1) = SIN(ANGLE)
I = I + 2
KK = KK + KINK
IF (KK.GT.NI) KK = KK - NI
50 CONTINUE
60 CONTINUE
*
RETURN
END
|
