C....*...1.........2.........3.........4.........5.........6.........7.*.......8
      SUBROUTINE REALTR(A,B,N,ISN)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
C     IF ISN=1, THIS SUBROUTINE COMPLETES THE FOURIER TRANSFORM
C       OF 2*N REAL DATA VALUES, WHERE THE ORIGINAL DATA VALUES ARE
C       STORED ALTERNATELY IN ARRAYS A AND B, AND ARE FIRST
C       TRANSFORMED BY A COMPLEX FOURIER TRANSFORM OF DIMENSION N.
C       THE COSINE COEFFICIENTS ARE IN A(1),A(2),...A(N+1) AND
C       THE SINE COEFFICIENTS ARE IN B(1),B(2),...B(N+1).
C       A TYPICAL CALLING SEQUENCE IS
C          CALL FFT(A,B,N,N,N,1)
C          CALL REALTR(A,B,N,1)
C       THE RESULTS SHOULD BE MULTIPLIED BY 0.5/N TO GIVE THE
C       USUAL SCALING OF COEFFICIENTS.
C     IF ISN=-1, THE INVERSE TRANSFORMATION IS DONE, THE FIRST STEP
C       IN EVALUATING A REAL FOURIER SERIES.
C       A TYPICAL CALLING SEQUENCE IS
C         CALL REALTR(A,B,N,-1)
C         CALL FFT(A,B,N,N,N,-1)
C       THE RESULTS SHOULD BE MULTIPLIED BY 0.5 TO GIVE THE USUAL
C       SCALING, AND THE TIME DOMAIN RESULTS ALTERNATE IN ARRAYS A
C       AND B, I.E. A(1),B(1),A(2),B(2),...A(N),B(N).
C     THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
C       ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO
C       GIVE THE CORRECT INDEXING INCREMENT AND A(2) USED TO
C       PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF IMAGINARY
C       VALUES, E.G.
C          CALL FFT(A,A(2),N,N,N,2)
C          CALL REALTR(A,A(2),N,2)
C       IN THIS CASE, THE COSINE AND SINE COEFFICIENTS ALTERNATE IN A.
C     BY R. C. SINGLETON, STANFORD RESEARCH INSTITUTE, OCT. 1968
      DIMENSION A(1),B(1)
      REAL*8 IM
      INC=IABS(ISN)
      NK=N*INC+2
      NH=NK/2
      SD=2.D0*DATAN(1.D0)/DFLOAT(N)
      CD=2.D0*DSIN(SD)**2
      SD=DSIN(SD+SD)
      SN=0.D0
      IF(ISN.LT.0) GO TO 30
      CN=1.D0
      A(NK-1)=A(1)
      B(NK-1)=B(1)
   10 DO 20 J=1,NH,INC
      K=NK-J
      AA=A(J)+A(K)
      AB=A(J)-A(K)
      BA=B(J)+B(K)
      BB=B(J)-B(K)
      RE=CN*BA+SN*AB
      IM=SN*BA-CN*AB
      B(K)=IM-BB
      B(J)=IM+BB
      A(K)=AA-RE
      A(J)=AA+RE
      AA=CN-(CD*CN+SD*SN)
      SN=(SD*CN-CD*SN)+SN
C     THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C       ERROR.  IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE
   20 CN=AA
C     CN=0.5D0/(AA**2+SN**2)+0.5D0
C     SN=CN*SN
C  20 CN=CN*AA
      RETURN
   30 CN=-1.D0
      SD=-SD
      GO TO 10
      END