subroutine realtr(a,b,n,isn)
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 1.0/(2.0*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 with most fortran compilers the data can alternatively be
c stored in a single complex array a, then the magnitude of isn
c changed to two to give the correct indexing increment and a(2)
c used to 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, sept. 1968
dimension a(1),b(1)
real im
inc=iabs(isn)
nk=n*inc+2
nh=nk/2
sd=2.0*atan(1.0)/float(n)
cd=2.0*sin(sd)**2
sd=sin(sd+sd)
sn=0.0
if(isn .lt. 0) go to 30
cn=1.0
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
cn=2.0-(aa**2+sn**2)
sn=cn*sn
20 cn=cn*aa
return
30 cn=-1.0
sd=-sd
go to 10
end