PROGRAM TRUESPIKE
REAL BESSMOD(30,-10:10) !Absolute values of Jn(x) for 10 arguements.
REAL RBESS(30,-700:700),IBESS(30,-700:700)
REAL RUNIT(30,-700:700),IUNIT(30,-700:700)
REAL A(30),PHI(30) !.......................!Amplit. and phase of
!30 Fourier Coeffs.
INTEGER GAMMA,GMAX,SEED
PI=4.0*ATAN(1.0)
C.... GMAX is highest spectral position reached, taken to be 700 here.
C.....GAMMA denotes which delta--fn of unit comb is point of observation.
GMAX=700
CALL RANDOM(CLOW,CHIGH,A,PHI,NMAX) !.......!Select Fourier Coeffs.
C.....Each comb has to begin with a complex value at each delta-fn due to
C.....the i^Mn multiplier but also the exp(Mn * PHI) from the phases of
C.....each Fourier Coefficient.
CALL BC1(A,PHI,RBESS,IBESS,BESSMOD)!.......!Initialise the Bessel
!.......!combs with their REAL
!.......!and IMAG. parts.
C.....First initialise the fundamental comb with its REAL (RUNIT) and
C.....IMAGINARY (IUNIT) parts:
DO 100 GAMMA=-6,6,1 !............!Move along comb orders.
RUNIT(1,GAMMA)=RBESS(1,GAMMA)
IUNIT(1,GAMMA)=IBESS(1,GAMMA)
100 CONTINUE
C... Now compute convolution of unit-spaced comb with the Mn comb.
DO 500 N=2,NMAX !..................!Convolve with next comb:
DO 450 GAMMA=-GMAX,GMAX,1 !........!Move point of observation:
DO 400 Mn=-6,6,1 !...........!Orders of N-comb:
M1=GAMMA-(Mn*N) !...........!Order of unit comb where
!...........!N-comb sits.
IF (ABS(M1).GE.GMAX) GOTO 400!.....!Cannot sit at frequency not
!.....!yet created by process:
AA=RUNIT(N-1,M1) !....!Real part of unit comb amplitude.
BB=IUNIT(N-1,M1) !........!Imag. part of unit comb amplitude.
XX=RBESS(N,Mn) !....!Real of N-comb amplitude
YY=IBESS(N,Mn) !........!Imag. part of N-comb amplitude.
BIT=((AA*XX)-(BB*YY)) !...........!Real part of complex
!...........!multiplication.
RUNIT(N,GAMMA)=RUNIT(N,GAMMA)+BIT !Add REAL additional amp.
!to total REAL amplitude there.
GIT=((BB*XX)+(AA*YY)) !...........!Imag. part of complex
IUNIT(N,GAMMA)=IUNIT(N,GAMMA)+GIT !Add IMAG additional amp.
!to total amplitude there.
400 CONTINUE
450 CONTINUE
500 CONTINUE
CALL SPECTLIST(A,RUNIT,IUNIT,NMAX) !.......!List Spectrum to file.
STOP
END
C.....**************************
SUBROUTINE ETEST(I,EVEN)
C.....Finds out if 'I' is an even or odd number.
EVEN=0
L=REAL(I)
A=MOD(L,2)
IF (A.EQ.0) EVEN=1
RETURN
END
C.....************************
SUBROUTINE RANDOM(CLOW,CHIGH,A,PHI,NMAX)
C.....Creates Fourier Coefficients of the phase object series.
REAL A(30),ZMIN(3),ZMAX(3),PHI(30)
INTEGER SEED
DATA(ZMIN(K),K=1,3)/0.0,0.4,1.1/
DATA(ZMAX(K),K=1,3)/0.4,1.1,1.6/
C.....This data defines three ranges of Fourier Coeff. values.
PI=4.0*ATAN(1.0)
WRITE(6,5)
5 FORMAT($,T2,'$Enter MINIMUM and MAXIMUM limits of An : ')
READ*,CLOW,CHIGH
WRITE(6,8)
8 FORMAT($,T2,'$Enter number of Fourier Coefficients : ')
READ*,NMAX
SEED=3112 !...............!Seed for Random Number Generator.
10 CONTINUE !...............!Marker to jump to if NEG is too big/small.
NEG=0 !...............!% of Fourier Coeffs. with phase<0.
SIGMA=16.0 !...............!Gaussian distribution of coeffs.
DO 50 I=1,NMAX
X=RAN(SEED)
ARG=(1.0*I/SIGMA)**2 !..!Gaussian frequency envelope .
A(I)=CHIGH*EXP(-0.5*ARG) !..!Amplitude of Coeff.
PHI(I)=-(0.95*PI)+(X*1.9*PI) !..!Phase of Coeff.
IF (PHI(I).LE.0) NEG=NEG+1
50 CONTINUE
C.....Check % of data with phase<zero.
ZNEG=(100.0*NEG)/(1.0*NMAX)
IF (ZNEG.GT.65.OR.ZNEG.LT.35) GOTO 10
PRINT*,'% of data set negative= ',ZNEG
RETURN
END
C.....************************
SUBROUTINE BC1(Z,PHI,RBESS,IBESS,BESSMOD)
C.....Creates the Bessel combs, each delta-fn. location is a complex number.
REAL Z(30),PHI(30),BESSMOD(30,-10:10)
REAL RBESS(30,-700:700),IBESS(30,-700:700)
C.....RBESS is 'n' of Zn, and order of function, real part of term.
C.....IBESS is same but Imaginary part of term.
C.....BESSMOD is modulus of Bessel Function.
DO 200 I=1,30 !.................!Step Fourier Coeffs. Comb no.
DO 100 M=6,-6,-1 !...........!Step Bessel Orders on each comb.
AMP=Z(I) !...........!I'th Fourier Coefficient.
CALL BC2(AMP,BESSMOD,I,M) !Find straight Bessel function Jm(Zi)
ALPHA=BESSMOD(I,M)*COS(M*PHI(I)) !...!REAL part of amplitude.
BETA=BESSMOD(I,M)*SIN(M*PHI(I)) !...!IMAG part of amplitude.
CALL ETEST(M,EVEN) !..........!Test if 'm' is odd or even no.
IF (EVEN.EQ.1) THEN !..........!then we have REAL part of term
ALPHA=ALPHA*((-1)**INT(M/2.0))
RBESS(I,M)=ALPHA
ELSE IF (EVEN.EQ.0) THEN !......!then we have IMAG part.
ALPHA=ALPHA*((-1)**INT((M-1.0)/2.0))
IBESS(I,M)=ALPHA
ENDIF
M2=M+1
CALL ETEST(M2,EVEN) !Test if 'm2' is odd or even no.
IF (EVEN.EQ.1) THEN !...........!then we have REAL part of term.
BETA=BETA*((-1)**INT(M2/2.0))
RBESS(I,M)=BETA
ELSE IF (EVEN.EQ.0) THEN
BETA=BETA*((-1)**INT((M2-1.0)/2.0))
IBESS(I,M)=BETA
ENDIF
100 CONTINUE
200 CONTINUE
RETURN
END
C.....************************
SUBROUTINE BC2(AMP,BESSMOD,I,MOO)
C.....Computes the Bessel function modulus given order and arguement.
REAL BESSMOD(30,-10:10),FAC(0:18)
c.....'I' is which of 30 coeffs. and 'MOO' is order of Bessel Fn.
C.....BESSMOD is amplitude of Jn, Bess('n' in An,Order).
DATA (FAC(J),J=0,18)/1,1,2,6,24,120,720,5040,
*40320,362880,3628800,39916800,479001600,
*6.2270208E9,8.71782912E10,1.3076E12,
*2.092278E13,3.5568E14,6.40237E15/
C... The factorial will be too large if J is much larger.
NEG=0
IF (MOO.LT.0) NEG=1
M=ABS(MOO) !Do for positive M and alter for -ve M at end of routine.
ASUM=0
P=+1
NMAX=18-M
DO 40 NLOOP=0,NMAX
DEN1=FAC(NLOOP)
DEN2=FAC(M+NLOOP)
NPOW=M+2*NLOOP
IF (NPOW.NE.0) THEN
A=(AMP/2.0)**(NPOW)
ELSE IF (NPOW.EQ.0) THEN
A=1.0
ENDIF
A=A/(DEN1*DEN2)
P=-1*P
A=A*P
ASUM=ASUM+A
CALL CRASH(A,NLOOP,NMAX)
40 CONTINUE
IF (NEG.EQ.0) THEN
BESSMOD(I,M)=ASUM
ELSE IF (NEG.EQ.1) THEN
CALL NEGBESS(BESSMOD,I,M)
ENDIF
C.....M is order of function.
RETURN
END
C.....*******************
SUBROUTINE CRASH(A,NLOOP,NMAX)
C... Stops numbers too small from crashing program.
V=ABS(A)
IF (V.LE.1E-8) NLOOP=NMAX
RETURN
END
C.....***********************
SUBROUTINE NEGBESS(BESSMOD,I,M)
C.....Takes care of negative order Bessel functions.
REAL BESSMOD(30,-10:10)
CALL ETEST(M,EVEN)
IF (EVEN.EQ.0) THEN
BESSMOD(I,-M)=-1*BESSMOD(I,M)
ELSE IF (EVEN.EQ.1) THEN
BESSMOD(I,-M)=BESSMOD(I,M)
ENDIF
RETURN
END
C.....********************
SUBROUTINE SPECTLIST(A,RUNIT,IUNIT,NMAX)
REAL RUNIT(30,-700:700),IUNIT(30,-700:700)
REAL A(30)
CHARACTER*8 NAME
NAME='SPIKE'
C.....Displays the spectrum, output suitable for the TOPDRAWER
C.....graphics package.
OPEN(UNIT=1,STATUS='UNKNOWN',FORM='FORMATTED',FILE=NAME)
WRITE(1,20)
DO 10 I=-70,70
VAL=SQRT(RUNIT(NMAX,I)**2+IUNIT(NMAX,I)**2)
IF (I.EQ.0) VAL=0.0
XA=I-0.5
XB=1.0*I
WRITE(1,*)XA,0.0
WRITE(1,*)XB,0.0
WRITE(1,*)XB,VAL
WRITE(1,*)XB,0.0
10 CONTINUE
20 FORMAT(1X,'SET WINDOW X 1 OF 1 Y 3 OF 3',/,
1 1X,'TITLE TOP"TRUESPIKE Results"',/,
2 1X,'TITLE BOTTOM"FREQUENCY"',/,
3 1X,'TITLE LEFT"Amplitude"')
WRITE(1,30)
30 FORMAT(1X,'JOIN 1')
CLOSE(1)
C.....Compute rms difference of spectrum from ideal.
DIFF=0.0
DO 40 I=-NMAX,NMAX,1
IF (I.EQ.0) GOTO 40
VAL=SQRT(RUNIT(NMAX,I)**2+IUNIT(NMAX,I)**2)
DIFF=(0.5*A(ABS(I))-VAL)**2
40 CONTINUE
DIFF=SQRT(DIFF/(2*NMAX-1.0))
PRINT*,'RMS spectrum difference from ideal is',DIFF
RETURN
END