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