/*-------------------------------------------------------------------- * $Id: fourt.c,v 1.2.4.2 2002/02/27 20:23:56 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* Translation of old fourt.f FORTRAN code to C using the automatic * translator f2c written by S.I. Feldman, David M. Gay, Mark W. Maimone, * and N.L. Schryer. Translated version provided by Andrew MacRae * at the University of Calgary. I've cleaned up the resulting code * since much of the f2c.h include stuff was unneccesary for this * function. * * P. Wessel *--------------------------------------------------------------------*/ #include "gmt.h" typedef int integer; typedef float real; typedef double doublereal; /* C-callable wrapper for fourt_ */ void GMT_fourt (float *data, int *nn, int ndim, int ksign, int iform, float *work) /* Data array */ /* Dimension array */ /* Number of dimensions */ /* Forward(-1) or Inverse(+1) */ /* Real(0) or complex(1) data */ /* Work array */ { int fourt_ (real *data, integer *nn, integer *ndim, integer *ksign, integer *iform, real *work); (void) fourt_ (data, nn, &ndim, &ksign, &iform, work); } int fourt_ (real *data, integer *nn, integer *ndim, integer *ksign, integer *iform, real *work) { /* System generated locals */ integer i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9, i__10, i__11, i__12; /* Builtin functions */ double sin(double); /* Local variables */ static integer j1rg2, idiv, irem, ipar, kmin, imin, jmin, lmax, mmax, imax, jmax; static integer ntwo, j1cnj, np1hf, np2hf, j1min, i1max, i1rng, j1rng, j2min, j3max; static integer j1max, j2max, i2max, non2t, j2stp, i, j, k, l, m, n, icase, ifact[1]; static integer nhalf, krang, kconj, kdif, idim, ntot, kstep, k2, k3, k4, nprev, iquot; static integer i2, i1, i3, j3, k1, j2, j1, if_, np0, np1, np2, ifp1, ifp2, non2; static real theta, oldsi, tempi, oldsr, sinth, difi, difr, sumi, sumr, tempr, twopi; static real wstpi, wstpr, twowr, wi, wr, u1i, w2i, w3i, u2i, u3i, u4i, t2i, u1r; static real u2r, u3r, w2r, w3r, u4r, t2r, t3r, t3i, t4r, t4i; static doublereal wrd, wid; /*--------------------------------------------------------------------------- ARGUMENTS : DATA - COMPLEX ARRAY, LENGTH NN NN - ARRAY OF NUMBER OF POINTS IN EACH DIMENSION NDIM - NUMBER OF DIMENSIONS (FOR OUR PURPOSES, NDIM=1) KSIGN - +1 FOR INVERSE TRANSFORM (FREQ TO TIME DOMAIN) -1 FOR FORWARD TRANSFORM (TIME TO FREQ DOMAIN) IFORM - 0 REAL DATA +1 COMPLEX DATA WORK - 0 IF ALL DIMENSIONS ARE RADIX 2 COMPLEX ARRAY, LARGE AS LARGEST NON-RADIX 2 DIMENSI0N PROGRAM BY NORMAN BRENNER FROM THE BASIC PROGRAM BY CHARLES RADER. RALPH ALTER SUGGESTED THE IDEA FOR THE DIGIT REVERSAL. MIT LINCOLN LABORATORY, AUGUST 1967. ---------------------------------------------------------------------------*/ /* Parameter adjustments */ --work; --nn; --data; /* Function Body */ wr = wi = wstpr = wstpi = (real)0.0; twopi = (real)6.283185307; if (*ndim - 1 >= 0) { goto L1; } else { goto L920; } L1: ntot = 2; i__1 = *ndim; for (idim = 1; idim <= i__1; ++idim) { if (nn[idim] <= 0) { goto L920; } else { goto L2; } L2: ntot *= nn[idim]; } np1 = 2; i__1 = *ndim; for (idim = 1; idim <= i__1; ++idim) { n = nn[idim]; np2 = np1 * n; if ((i__2 = n - 1) < 0) { goto L920; } else if (i__2 == 0) { goto L900; } else { goto L5; } L5: m = n; ntwo = np1; if_ = 1; idiv = 2; L10: iquot = m / idiv; irem = m - idiv * iquot; if (iquot - idiv >= 0) { goto L11; } else { goto L50; } L11: if (irem != 0) { goto L20; } else { goto L12; } L12: ntwo += ntwo; m = iquot; goto L10; L20: idiv = 3; L30: iquot = m / idiv; irem = m - idiv * iquot; if (iquot - idiv >= 0) { goto L31; } else { goto L60; } L31: if (irem != 0) { goto L40; } else { goto L32; } L32: ifact[if_ - 1] = idiv; ++if_; m = iquot; goto L30; L40: idiv += 2; goto L30; L50: if (irem != 0) { goto L60; } else { goto L51; } L51: ntwo += ntwo; goto L70; L60: ifact[if_ - 1] = m; L70: non2 = np1 * (np2 / ntwo); icase = 1; if (idim - 4 >= 0) { goto L90; } else { goto L71; } L71: if (*iform <= 0) { goto L72; } else { goto L90; } L72: icase = 2; if (idim - 1 <= 0) { goto L73; } else { goto L90; } L73: icase = 3; if (ntwo - np1 <= 0) { goto L90; } else { goto L74; } L74: icase = 4; ntwo /= 2; n /= 2; np2 /= 2; ntot /= 2; i = 3; i__2 = ntot; for (j = 2; j <= i__2; ++j) { data[j] = data[i]; /* L80: */ i += 2; } L90: i1rng = np1; if (icase - 2 != 0) { goto L100; } else { goto L95; } L95: i1rng = np0 * (nprev / 2 + 1); L100: if (ntwo - np1 <= 0) { goto L600; } else { goto L110; } L110: np2hf = np2 / 2; j = 1; i__2 = np2; i__3 = non2; for (i2 = 1; i__3 < 0 ? i2 >= i__2 : i2 <= i__2; i2 += i__3) { if (j - i2 >= 0) { goto L130; } else { goto L120; } L120: i1max = i2 + non2 - 2; i__4 = i1max; for (i1 = i2; i1 <= i__4; i1 += 2) { i__5 = ntot; i__6 = np2; for (i3 = i1; i__6 < 0 ? i3 >= i__5 : i3 <= i__5; i3 += i__6) { j3 = j + i3 - i2; tempr = data[i3]; tempi = data[i3 + 1]; data[i3] = data[j3]; data[i3 + 1] = data[j3 + 1]; data[j3] = tempr; /* L125: */ data[j3 + 1] = tempi; } } L130: m = np2hf; L140: if (j - m <= 0) { goto L150; } else { goto L145; } L145: j -= m; m /= 2; if (m - non2 >= 0) { goto L140; } else { goto L150; } L150: j += m; } non2t = non2 + non2; ipar = ntwo / np1; L310: if ((i__3 = ipar - 2) < 0) { goto L350; } else if (i__3 == 0) { goto L330; } else { goto L320; } L320: ipar /= 4; goto L310; L330: i__3 = i1rng; for (i1 = 1; i1 <= i__3; i1 += 2) { i__2 = non2; i__6 = np1; for (j3 = i1; i__6 < 0 ? j3 >= i__2 : j3 <= i__2; j3 += i__6) { i__5 = ntot; i__4 = non2t; for (k1 = j3; i__4 < 0 ? k1 >= i__5 : k1 <= i__5; k1 += i__4) { k2 = k1 + non2; tempr = data[k2]; tempi = data[k2 + 1]; data[k2] = data[k1] - tempr; data[k2 + 1] = data[k1 + 1] - tempi; data[k1] += tempr; /* L340: */ data[k1 + 1] += tempi; } } } L350: mmax = non2; L360: if (mmax - np2hf >= 0) { goto L600; } else { goto L370; } L370: /* Computing MAX */ i__4 = non2t, i__5 = mmax / 2; lmax = MAX(i__4,i__5); if (mmax - non2 <= 0) { goto L405; } else { goto L380; } L380: theta = -twopi * (real) non2 / (real) (mmax << 2); if (*ksign >= 0) { goto L390; } else { goto L400; } L390: theta = -theta; L400: sincos ((doublereal)theta, &wid, &wrd); wr = (real)wrd; wi = (real)wid; wstpr = (real)-2.0 * wi * wi; wstpi = (real)2.0 * wr * wi; L405: i__4 = lmax; i__5 = non2t; for (l = non2; i__5 < 0 ? l >= i__4 : l <= i__4; l += i__5) { m = l; if (mmax - non2 <= 0) { goto L420; } else { goto L410; } L410: w2r = wr * wr - wi * wi; w2i = (real)(wr * 2.0 * wi); w3r = w2r * wr - w2i * wi; w3i = w2r * wi + w2i * wr; L420: i__6 = i1rng; for (i1 = 1; i1 <= i__6; i1 += 2) { i__2 = non2; i__3 = np1; for (j3 = i1; i__3 < 0 ? j3 >= i__2 : j3 <= i__2; j3 += i__3) { kmin = j3 + ipar * m; if (mmax - non2 <= 0) { goto L430; } else { goto L440; } L430: kmin = j3; L440: kdif = ipar * mmax; L450: kstep = kdif << 2; i__7 = ntot; i__8 = kstep; for (k1 = kmin; i__8 < 0 ? k1 >= i__7 : k1 <= i__7; k1 += i__8) { k2 = k1 + kdif; k3 = k2 + kdif; k4 = k3 + kdif; if (mmax - non2 <= 0) { goto L460; } else { goto L480; } L460: u1r = data[k1] + data[k2]; u1i = data[k1 + 1] + data[k2 + 1]; u2r = data[k3] + data[k4]; u2i = data[k3 + 1] + data[k4 + 1]; u3r = data[k1] - data[k2]; u3i = data[k1 + 1] - data[k2 + 1]; if (*ksign >= 0) { goto L475; } else { goto L470; } L470: u4r = data[k3 + 1] - data[k4 + 1]; u4i = data[k4] - data[k3]; goto L510; L475: u4r = data[k4 + 1] - data[k3 + 1]; u4i = data[k3] - data[k4]; goto L510; L480: t2r = w2r * data[k2] - w2i * data[k2 + 1]; t2i = w2r * data[k2 + 1] + w2i * data[k2]; t3r = wr * data[k3] - wi * data[k3 + 1]; t3i = wr * data[k3 + 1] + wi * data[k3]; t4r = w3r * data[k4] - w3i * data[k4 + 1]; t4i = w3r * data[k4 + 1] + w3i * data[k4]; u1r = data[k1] + t2r; u1i = data[k1 + 1] + t2i; u2r = t3r + t4r; u2i = t3i + t4i; u3r = data[k1] - t2r; u3i = data[k1 + 1] - t2i; if (*ksign >= 0) { goto L500; } else { goto L490; } L490: u4r = t3i - t4i; u4i = t4r - t3r; goto L510; L500: u4r = t4i - t3i; u4i = t3r - t4r; L510: data[k1] = u1r + u2r; data[k1 + 1] = u1i + u2i; data[k2] = u3r + u4r; data[k2 + 1] = u3i + u4i; data[k3] = u1r - u2r; data[k3 + 1] = u1i - u2i; data[k4] = u3r - u4r; /* L520: */ data[k4 + 1] = u3i - u4i; } kmin = ((kmin - j3) << 2) + j3; kdif = kstep; if (kdif - np2 >= 0) { goto L530; } else { goto L450; } L530: ; } } m = mmax - m; if (*ksign >= 0) { goto L550; } else { goto L540; } L540: tempr = wr; wr = -wi; wi = -tempr; goto L560; L550: tempr = wr; wr = wi; wi = tempr; L560: if (m - lmax <= 0) { goto L565; } else { goto L410; } L565: tempr = wr; wr = wr * wstpr - wi * wstpi + wr; /* L570: */ wi = wi * wstpr + tempr * wstpi + wi; } ipar = 3 - ipar; mmax += mmax; goto L360; L600: if (ntwo - np2 >= 0) { goto L700; } else { goto L605; } L605: ifp1 = non2; if_ = 1; np1hf = np1 / 2; L610: ifp2 = ifp1 / ifact[if_ - 1]; j1rng = np2; if (icase - 3 != 0) { goto L612; } else { goto L611; } L611: j1rng = (np2 + ifp1) / 2; j2stp = np2 / ifact[if_ - 1]; j1rg2 = (j2stp + ifp2) / 2; L612: j2min = ifp2 + 1; if (ifp1 - np2 >= 0) { goto L640; } else { goto L615; } L615: i__5 = ifp1; i__4 = ifp2; for (j2 = j2min; i__4 < 0 ? j2 >= i__5 : j2 <= i__5; j2 += i__4) { theta = -twopi * (real) (j2 - 1) / (real) np2; if (*ksign >= 0) { goto L620; } else { goto L625; } L620: theta = -theta; L625: sinth = (real)sin((double)(0.5 * theta)); wstpr = sinth * (real)(-2. * sinth); wstpi = (real)sin((double)theta); wr = wstpr + (real)1.0; wi = wstpi; j1min = j2 + ifp1; i__3 = j1rng; i__2 = ifp1; for (j1 = j1min; i__2 < 0 ? j1 >= i__3 : j1 <= i__3; j1 += i__2) { i1max = j1 + i1rng - 2; i__6 = i1max; for (i1 = j1; i1 <= i__6; i1 += 2) { i__8 = ntot; i__7 = np2; for (i3 = i1; i__7 < 0 ? i3 >= i__8 : i3 <= i__8; i3 += i__7) { j3max = i3 + ifp2 - np1; i__9 = j3max; i__10 = np1; for (j3 = i3; i__10 < 0 ? j3 >= i__9 : j3 <= i__9; j3 += i__10) { tempr = data[j3]; data[j3] = data[j3] * wr - data[j3 + 1] * wi; /* L630: */ data[j3 + 1] = tempr * wi + data[j3 + 1] * wr; } } } tempr = wr; wr = wr * wstpr - wi * wstpi + wr; /* L635: */ wi = tempr * wstpi + wi * wstpr + wi; } } L640: theta = -twopi / (real) ifact[if_ - 1]; if (*ksign >= 0) { goto L645; } else { goto L650; } L645: theta = -theta; L650: sinth = (real)sin((double)(0.5 * theta)); wstpr = sinth * (real)(-2. * sinth); wstpi = (real)sin((double)theta); kstep = (n << 1) / ifact[if_ - 1]; krang = kstep * (ifact[if_ - 1] / 2) + 1; i__2 = i1rng; for (i1 = 1; i1 <= i__2; i1 += 2) { i__3 = ntot; i__4 = np2; for (i3 = i1; i__4 < 0 ? i3 >= i__3 : i3 <= i__3; i3 += i__4) { i__5 = krang; i__10 = kstep; for (kmin = 1; i__10 < 0 ? kmin >= i__5 : kmin <= i__5; kmin += i__10) { j1max = i3 + j1rng - ifp1; i__9 = j1max; i__7 = ifp1; for (j1 = i3; i__7 < 0 ? j1 >= i__9 : j1 <= i__9; j1 += i__7) { j3max = j1 + ifp2 - np1; i__8 = j3max; i__6 = np1; for (j3 = j1; i__6 < 0 ? j3 >= i__8 : j3 <= i__8; j3 += i__6) { j2max = j3 + ifp1 - ifp2; k = kmin + (j3 - j1 + (j1 - i3) / ifact[if_ - 1]) / np1hf; if (kmin - 1 <= 0) { goto L655; } else { goto L665; } L655: sumr = (real)0.0; sumi = (real)0.0; i__11 = j2max; i__12 = ifp2; for (j2 = j3; i__12 < 0 ? j2 >= i__11 : j2 <= i__11; j2 += i__12) { sumr += data[j2]; /* L660: */ sumi += data[j2 + 1]; } work[k] = sumr; work[k + 1] = sumi; goto L680; L665: kconj = k + ((n - kmin + 1) << 1); j2 = j2max; sumr = data[j2]; sumi = data[j2 + 1]; oldsr = (real)0.0; oldsi = (real)0.0; j2 -= ifp2; L670: tempr = sumr; tempi = sumi; sumr = twowr * sumr - oldsr + data[j2]; sumi = twowr * sumi - oldsi + data[j2 + 1]; oldsr = tempr; oldsi = tempi; j2 -= ifp2; if (j2 - j3 <= 0) { goto L675; } else { goto L670; } L675: tempr = wr * sumr - oldsr + data[j2]; tempi = wi * sumi; work[k] = tempr - tempi; work[kconj] = tempr + tempi; tempr = wr * sumi - oldsi + data[j2 + 1]; tempi = wi * sumr; work[k + 1] = tempr + tempi; work[kconj + 1] = tempr - tempi; L680: ; } } if (kmin - 1 <= 0) { goto L685; } else { goto L686; } L685: wr = wstpr + (real)1.0; wi = wstpi; goto L690; L686: tempr = wr; wr = wr * wstpr - wi * wstpi + wr; wi = tempr * wstpi + wi * wstpr + wi; L690: twowr = wr + wr; } if (icase - 3 != 0) { goto L692; } else { goto L691; } L691: if (ifp1 - np2 >= 0) { goto L692; } else { goto L695; } L692: k = 1; i2max = i3 + np2 - np1; i__10 = i2max; i__5 = np1; for (i2 = i3; i__5 < 0 ? i2 >= i__10 : i2 <= i__10; i2 += i__5) { data[i2] = work[k]; data[i2 + 1] = work[k + 1]; /* L693: */ k += 2; } goto L698; L695: j3max = i3 + ifp2 - np1; i__5 = j3max; i__10 = np1; for (j3 = i3; i__10 < 0 ? j3 >= i__5 : j3 <= i__5; j3 += i__10) { j2max = j3 + np2 - j2stp; i__6 = j2max; i__8 = j2stp; for (j2 = j3; i__8 < 0 ? j2 >= i__6 : j2 <= i__6; j2 += i__8) { j1max = j2 + j1rg2 - ifp2; j1cnj = j3 + j2max + j2stp - j2; i__7 = j1max; i__9 = ifp2; for (j1 = j2; i__9 < 0 ? j1 >= i__7 : j1 <= i__7; j1 += i__9) { k = j1 + 1 - i3; data[j1] = work[k]; data[j1 + 1] = work[k + 1]; if (j1 - j2 <= 0) { goto L697; } else { goto L696; } L696: data[j1cnj] = work[k]; data[j1cnj + 1] = -work[k + 1]; L697: j1cnj -= ifp2; } } } L698: ; } } ++if_; ifp1 = ifp2; if (ifp1 - np1 <= 0) { goto L700; } else { goto L610; } L700: switch ((int)icase) { case 1: goto L900; case 2: goto L800; case 3: goto L900; case 4: goto L701; } L701: nhalf = n; n += n; theta = -twopi / (real) n; if (*ksign >= 0) { goto L702; } else { goto L703; } L702: theta = -theta; L703: sinth = (real)sin((double)(0.5 * theta)); wstpr = sinth * (real)(-2. * sinth); wstpi = (real)sin((double)theta); wr = wstpr + (real)1.0; wi = wstpi; imin = 3; jmin = (nhalf << 1) - 1; goto L725; L710: j = jmin; i__4 = ntot; i__3 = np2; for (i = imin; i__3 < 0 ? i >= i__4 : i <= i__4; i += i__3) { sumr = (real)0.5 * (data[i] + data[j]); sumi = (real)0.5 * (data[i + 1] + data[j + 1]); difr = (real)0.5 * (data[i] - data[j]); difi = (real)0.5 * (data[i + 1] - data[j + 1]); tempr = wr * sumi + wi * difr; tempi = wi * sumi - wr * difr; data[i] = sumr + tempr; data[i + 1] = difi + tempi; data[j] = sumr - tempr; data[j + 1] = -difi + tempi; /* L720: */ j += np2; } imin += 2; jmin += -2; tempr = wr; wr = wr * wstpr - wi * wstpi + wr; wi = tempr * wstpi + wi * wstpr + wi; L725: if ((i__3 = imin - jmin) < 0) { goto L710; } else if (i__3 == 0) { goto L730; } else { goto L740; } L730: if (*ksign >= 0) { goto L740; } else { goto L731; } L731: i__3 = ntot; i__4 = np2; for (i = imin; i__4 < 0 ? i >= i__3 : i <= i__3; i += i__4) { /* L735: */ data[i + 1] = -data[i + 1]; } L740: np2 += np2; ntot += ntot; j = ntot + 1; imax = ntot / 2 + 1; L745: imin = imax - (nhalf << 1); i = imin; goto L755; L750: data[j] = data[i]; data[j + 1] = -data[i + 1]; L755: i += 2; j += -2; if (i - imax >= 0) { goto L760; } else { goto L750; } L760: data[j] = data[imin] - data[imin + 1]; data[j + 1] = (real)0.0; if (i - j >= 0) { goto L780; } else { goto L770; } L765: data[j] = data[i]; data[j + 1] = data[i + 1]; L770: i += -2; j += -2; if (i - imin <= 0) { goto L775; } else { goto L765; } L775: data[j] = data[imin] + data[imin + 1]; data[j + 1] = (real)0.0; imax = imin; goto L745; L780: data[1] += data[2]; data[2] = (real)0.0; goto L900; L800: if (i1rng - np1 >= 0) { goto L900; } else { goto L805; } L805: i__4 = ntot; i__3 = np2; for (i3 = 1; i__3 < 0 ? i3 >= i__4 : i3 <= i__4; i3 += i__3) { i2max = i3 + np2 - np1; i__2 = i2max; i__9 = np1; for (i2 = i3; i__9 < 0 ? i2 >= i__2 : i2 <= i__2; i2 += i__9) { imin = i2 + i1rng; imax = i2 + np1 - 2; jmax = (i3 << 1) + np1 - imin; if (i2 - i3 <= 0) { goto L820; } else { goto L810; } L810: jmax += np2; L820: if (idim - 2 <= 0) { goto L850; } else { goto L830; } L830: j = jmax + np0; i__7 = imax; for (i = imin; i <= i__7; i += 2) { data[i] = data[j]; data[i + 1] = -data[j + 1]; /* L840: */ j += -2; } L850: j = jmax; i__7 = imax; i__8 = np0; for (i = imin; i__8 < 0 ? i >= i__7 : i <= i__7; i += i__8) { data[i] = data[j]; data[i + 1] = -data[j + 1]; /* L860: */ j -= np0; } } } L900: np0 = np1; np1 = np2; /* L910: */ nprev = n; } L920: return 0; } /* fourt_ */