/* realtr.f -- translated by f2c (version of 23 April 1993  18:34:30).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Subroutine */ int realtr_(real *a, real *b, integer *n, integer *isn)
{
    /* System generated locals */
    integer i__1, i__2;
    real r__1, r__2;

    /* Builtin functions */
    double atan(doublereal), sin(doublereal);

    /* Local variables */
    integer j, k;
    real aa, ab, ba, bb, cd, cn, im;
    integer nh;
    real sd;
    integer nk;
    real re, sn;
    integer inc;

/*  if isn=1, this subroutine completes the fourier transform */
/*    of 2*n real data values, where the original data values are */
/*    stored alternately in arrays a and b, and are first */
/*    transformed by a complex fourier transform of dimension n. */
/*    the cosine coefficients are in a(1),a(2),...,a(n+1) and */
/*    the sine coefficients are in b(1),b(2),...,b(n+1). */
/*    a typical calling sequence is */
/*      call fft(a,b,n,n,n,1) */
/*      call realtr(a,b,n,1) */
/*    the results should be multiplied by 1.0/(2.0*n) to give the */
/*    usual scaling of coefficients. */
/*  if isn=-1, the inverse transformation is done, the first step */
/*    in evaluating a real fourier series. */
/*    a typical calling sequence is */
/*      call realtr(a,b,n,-1) */
/*      call fft(a,b,n,n,n,-1) */
/*    the results should be multiplied by 0.5 to give the usual */
/*    scaling, and the time domain results alternate in arrays a */
/*    and b, i.e. a(1),b(1),a(2),b(2),...,a(n),b(n). */
/*  with most fortran compilers the data can alternatively be */
/*    stored in a single complex array a, then the magnitude of isn */
/*    changed to two to give the correct indexing increment and a(2) */
/*    used to pass the initial address for the sequence of imaginary */
/*    values, e.g. */
/*      call fft(a,a(2),n,n,n,2) */
/*      call realtr(a,a(2),n,2) */
/*    in this case, the cosine and sine coefficients alternate in a. */
/*  by r. c. singleton, stanford research institute, sept. 1968 */
    /* Parameter adjustments */
    --b;
    --a;

    /* Function Body */
    inc = abs(*isn);
    nk = *n * inc + 2;
    nh = nk / 2;
    sd = atan(1.f) * 2.f / (real) (*n);
/* Computing 2nd power */
    r__1 = sin(sd);
    cd = r__1 * r__1 * 2.f;
    sd = sin(sd + sd);
    sn = 0.f;
    if (*isn < 0) {
	goto L30;
    }
    cn = 1.f;
    a[nk - 1] = a[1];
    b[nk - 1] = b[1];
L10:
    i__1 = nh;
    i__2 = inc;
    for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
	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;
/* Computing 2nd power */
	r__1 = aa;
/* Computing 2nd power */
	r__2 = sn;
	cn = 2.f - (r__1 * r__1 + r__2 * r__2);
	sn = cn * sn;
/* L20: */
	cn *= aa;
    }
    return 0;
L30:
    cn = -1.f;
    sd = -(doublereal)sd;
    goto L10;
} /* realtr_ */