home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Frostbyte's 1980s DOS Shareware Collection
/
floppyshareware.zip
/
floppyshareware
/
DOOG
/
PCSSP1.ZIP
/
ITRPAPSM.ZIP
/
DRHARM.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
4KB
|
128 lines
C
C ..................................................................
C
C SUBROUTINE DRHARM
C
C PURPOSE
C FINDS THE FOURIER COEFFICIENTS OF ONE DIMENSIONAL DOUBLE
C PRECISION REAL DATA
C
C USAGE
C CALL DRHARM(A,M,INV,S,IFERR)
C
C DESCRIPTION OF PARAMETERS
C A - A DOUBLE PRECISION VECTOR
C AS INPUT, CONTAINS ONE DIMENSIONAL REAL DATA. A IS
C 2*N+4 CORE LOCATIONS, WHERE N = 2**M. 2*N REAL
C NUMBERS ARE PUT INTO THE FIRST 2*N CORE LOCATIONS
C OF A
C AS OUTPUT, A CONTAINS THE FOURIER COEFFICIENTS
C A0/2,B0=0,A1,B1,A2,B2,...,AN/2,BN=0 RESPECTIVELY IN
C THE FIRST 2N+2 CORE LOCATIONS OF A
C M - AN INTEGER WHICH DETERMINES THE SIZE OF THE VECTOR
C A. THE SIZE OF A IS 2*(2**M) + 4
C INV - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION OF
C DIMENSION ONE EIGHTH THE NUMBER OF REAL INPUT, VIZ.,
C (1/8)*2*(2**M)
C S - A DOUBLE PRECISION VECTOR WORK AREA FOR SINE TABLES
C WITH DIMENSION THE SAME AS INV
C IFERR - A RETURNED VALUE OF 1 MEANS THAT M IS LESS THAN 3 OR
C GREATER THAN 20. OTHERWISE IFERR IS SET = 0
C
C REMARKS
C THIS SUBROUTINE GIVES THE FOURIER COEFFICIENTS OF 2*(2**M)
C REAL POINTS. SEE SUBROUTINE DHARM FOR THREE DIMENSIONAL,
C DOUBLE PRECISION, COMPLEX FOURIER TRANSFORMS.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C DHARM
C
C METHOD
C THE FOURIER COEFFICIENTS A0,B0=0,A1,B1,...,AN,BN=0 ARE
C OBTAINED FOR INPUT XJ, J=0,1,2,...,2N-1 FOR THE FOLLOWING
C EQUATION (PI = 3.14159...)
C
C N-1 J
C XJ=(1/2)A0+SUM (AK*COS(PI*J*K/N)+BK*SIN(PI*J*K/N))+(1/2)AN(-1)
C K=1
C
C SEE REFERENCE UNDER SUBROUTINE DHARM
C
C ..................................................................
C
SUBROUTINE DRHARM(A,M,INV,S,IFERR)
DIMENSION A(1),L(3),INV(1),S(1)
DOUBLE PRECISION A,SI,AP1IM,FN,CO,CIRE,AP2IM,S,SS,DEL,CIIM,AP1RE,
1 CNIRE,SC,SIS,AP2RE,CNIIM
IFSET=1
L(1)=M
L(2)=0
L(3)=0
NTOT=2**M
NTOT2 = 2*NTOT
FN = NTOT
DO 3 I = 2,NTOT2,2
3 A(I) = -A(I)
DO 6 I = 1,NTOT2
6 A(I) = A(I)/FN
CALL DHARM(A,L,INV,S,IFSET,IFERR)
C
C MOVE LAST HALF OF A(J)S DOWN ONE SLOT AND ADD A(N) AT BOTTOM TO
C GIVE ARRAY FOR A1PRIME AND A2PRIME CALCULATION
C
21 DO 52 I=1,NTOT,2
J0=NTOT2+2-I
A(J0)=A(J0-2)
52 A(J0+1)=A(J0-1)
A(NTOT2+3)=A(1)
A(NTOT2+4)=A(2)
C
C CALCULATE A1PRIMES AND STORE IN FIRST N SLOTS
C CALCULATE A2PRIMES AND STORE IN SECOND N SLOTS IN REVERSE ORDER
K0=NTOT+1
DO 104 I=1,K0,2
K1=NTOT2-I+4
AP1RE=.5*(A(I)+A(K1))
AP2RE=-.5*(A(I+1)+A(K1+1))
AP1IM=.5*(-A(I+1)+A(K1+1))
AP2IM=-.5*(A(I)-A(K1))
A(I)=AP1RE
A(I+1)=AP1IM
A(K1)=AP2RE
104 A(K1+1)=AP2IM
NTO = NTOT/2
110 NT=NTO+1
DEL=3.141592653589793/DFLOAT(NTOT)
SS=DSIN(DEL)
SC=DCOS(DEL)
SI=0.0
CO=1.0
C
C COMPUTE C(J)S FOR J=0 THRU J=N
114 DO 116 I=1,NT
K6=NTOT2-2*I+5
AP2RE=A(K6)*CO+A(K6+1)*SI
AP2IM=-A(K6)*SI+A(K6+1)*CO
CIRE=.5*(A(2*I-1)+AP2RE)
CIIM=.5*(A(2*I)+AP2IM)
CNIRE=.5*(A(2*I-1)-AP2RE)
CNIIM=.5*(A(2*I)-AP2IM)
A(2*I-1)=CIRE
A(2*I)=CIIM
A(K6)=CNIRE
A(K6+1)=-CNIIM
SIS=SI
SI=SI*SC+CO*SS
116 CO=CO*SC-SIS*SS
C
C SHIFT C(J)S FOR J=N/2+1 TO J=N UP ONE SLOT
DO 117 I=1,NTOT,2
K8=NTOT+4+I
A(K8-2)=A(K8)
117 A(K8-1)=A(K8+1)
DO 500 I=3,NTOT2,2
A(I) = 2. * A(I)
500 A(I + 1) = -2. * A(I + 1)
RETURN
END