/******************************************************************************/
/*                                                                            */
/*                   ---- NIH NMR Software System ----                        */
/*                        Copyright 1992 and 1993                             */
/*                             Frank Delaglio                                 */
/*                   NIH Laboratory of Chemical Physics                       */
/*                                                                            */
/*               This software is not for distribution without                */
/*                  the written permission of the author.                     */
/*                                                                            */
/******************************************************************************/

/***/
/* userProc: User 1D Process Functions; process and write a 1D slice.
/***/

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>

#define USE_CHAN_DEFS

#include "cmndargs.h"
#include "specunit.h"
#include "memory.h"
#include "dataio.h"
#include "dimloc.h"
#include "fdatap.h"

#include "namelist.h"
#include "nmrserver.h"
#include "nmrpipe.h"
#include "userproc.h"
#include "prec.h"

#include "fdm.h"

#define ISQUAD( D )  (1 != (int)getQuad( dataInfo->fdata, NDQUADFLAG, (D) ))

#define FPR    (void)fprintf

#define ARGC   dataInfo->argc 
#define ARGV   dataInfo->argv 
#define HDR    dataInfo->fdata
#define DI     dataInfo

static int     vAlt(), vNeg();
static int     rSet(), rMove();
static int     slowFFT(), slowIFT(), nextPower2();

/***/
/* getFuncList: returns pointer to UPROC function table.
/***/

struct ProcFuncInfo *getFuncList()
{
    return( uprocList );
}

int PASTDATE()
{
   return( 0 );
}

int lnType()
{
   return( 0 );
}

/***/
/* uNULL: null 1D processing function.
/***/

int uNULL( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    int error;

    error = 0;

/***/
/* Process, Write:
/***/

    if (dataInfo->sliceCode > 0)
       {
        error = putSlice( dataInfo );
       }

/***/
/* Command Line:
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        if (flagLoc( ARGC, ARGV, "-help" ))
           {
            FPR( stderr, "NULL: No Change to Data.\n" );
            return( 1 );
           }
       }

/***/
/* Initialization and Header:
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
       }

/***/
/* Shutdown:
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
       }

/***/
/* Unknown mode:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* uFT: 1D Fourier Transform.
/***/

int uFT( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    static int invFlag, realFlag, altFlag, negFlag;
    int    tdSize, error;

    error = 0;

/***/
/* Process, Write:
/*   Sign-alternate if needed for forward transform.
/*   Zero imags for real transform if needed.
/*   Perform FFT.
/*   Extract right half of result for real transform.
/*   Sign-alternate if needed for inverse transform.
/*   Write the result.
/***/

    if (dataInfo->sliceCode > 0)
       {
        if (altFlag && !invFlag)
           {
            (void) vAlt( dataInfo->rdata, dataInfo->size );
            (void) vAlt( dataInfo->idata, dataInfo->size );
           }

        if (negFlag && !invFlag)
           {
            (void) vNeg( dataInfo->idata, dataInfo->outSize );
           }

        if (realFlag && dataInfo->quadState != 1) 
           {
            (void) rSet( dataInfo->idata, dataInfo->size, (float) 0.0 );
           }

        if (invFlag)
           {
            (void) slowIFT( dataInfo->rdata, dataInfo->idata, 
                            dataInfo->size );
           }
        else
           {
            (void) slowFFT( dataInfo->rdata, dataInfo->idata, 
                            dataInfo->size );
           }

        if (realFlag)
           {
            (void) rMove( &dataInfo->rdata[dataInfo->size/2],
                          dataInfo->rdata,
                          dataInfo->outSize );

            (void) rMove( &dataInfo->idata[dataInfo->size/2],
                          dataInfo->idata,
                          dataInfo->outSize );
           }

        if (altFlag && invFlag)
           {
            (void) vAlt( dataInfo->rdata, dataInfo->outSize );
            (void) vAlt( dataInfo->idata, dataInfo->outSize );
           }

        if (negFlag && invFlag)
           {
            (void) vNeg( dataInfo->idata, dataInfo->outSize );
           }

        error = putSlice( dataInfo );
       }

/***/
/* Command line:
/*   Extract FT modes.
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        if (flagLoc( ARGC, ARGV, "-help" ))
           {
            FPR( stderr, "FT: Complex Fourier Transform.\n" );
            FPR( stderr, " -real  Transform of Real-Only Data.\n" );
            FPR( stderr, " -inv   Perform Inverse Transform.\n" );
            FPR( stderr, " -alt   Use Sign Alternation.\n" );
            FPR( stderr, " -neg   Negate Imaginaries.\n" );
            return( 1 );
           }

        realFlag = 0;
        invFlag  = 0;
        altFlag  = 0;
        negFlag  = 0;
        
        if (flagLoc( ARGC, ARGV, "-real" ))   realFlag = 1;
        if (flagLoc( ARGC, ARGV, "-inv" ))    invFlag  = 1;
        if (flagLoc( ARGC, ARGV, "-alt" ))    altFlag  = 1;
        if (flagLoc( ARGC, ARGV, "-neg" ))    negFlag  = 1;
       }

/***/
/* Initialization, Header:
/*   Data will be complex after FT.
/*   Adjust FT state of dimension.
/*   Reduce output size in case of real transform.
/*   Allocate input, output and trig buffers.
/*
/* DMX handling:
/*  Auto mode, test whether DMX adjustment should be performed.
/*  If adjustment is needed, initialize phase buffers. 
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
        if (dataInfo->ftFlag)
           dataInfo->ftFlag = 0;
        else
           dataInfo->ftFlag = 1;

        (void) setParm( HDR, NDFTFLAG, (float)dataInfo->ftFlag, CUR_XDIM );

        if (dataInfo->ftFlag)
           {
            (void) setParm( HDR, NDAQSIGN, (float)ALT_NONE,       CUR_XDIM );
            (void) setParm( HDR, NDFTSIZE, (float)dataInfo->size, CUR_XDIM );
           }

        (void) setParm( HDR, FDQUADFLAG, (float)0.0, CUR_XDIM );
        (void) setParm( HDR, NDQUADFLAG, (float)0.0, CUR_XDIM );

        dataInfo->outQuadState = 2;

        if (realFlag)
           {
            tdSize = getParm( HDR, NDTDSIZE, CUR_XDIM );

            dataInfo->outSize  = dataInfo->size/2;
            dataInfo->timeSize = dataInfo->timeSize/2;
            tdSize             = tdSize/2;
          
            (void)setParm( HDR, NDSIZE,   (float)dataInfo->outSize,  CUR_XDIM );
            (void)setParm( HDR, NDAPOD,   (float)dataInfo->timeSize, CUR_XDIM );
            (void)setParm( HDR, NDTDSIZE, (float)tdSize,             CUR_XDIM );
           }
       }

/***/
/* Shutdown: 
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
        error = 0;
       }

/***/
/* Unknown code:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* uZF: 1D Zero Fill
/***/

int uZF( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    int   i, mid, ix1, ixn, izf, fSize, error;
    float orig, sw, obs, car, *rPtr;
 
    static int zfSize, zfCount, padCount, midFlag, invFlag, autoFlag;

    error = 0;

/***/
/* Process, Write:
/*  Process only for actual zero fill, as opposed to extraction.
/***/

    if (dataInfo->sliceCode > 0)
       {
        if (dataInfo->outSize > dataInfo->size)
           {
            (void) rSet( dataInfo->rdata + DI->size, zfSize - DI->size, 0.0 );
            (void) rSet( dataInfo->idata + DI->size, zfSize - DI->size, 0.0 );
           }
 
        error = putSlice( dataInfo );
       }

/***/
/* Command Line:
/*  Extract desired size after ZF; default is 2*size.
/*  Explicit Final Size:      -size finalSize
/*  Extra Zeros by Pad Count: -pad  padCount 
/*  Final Size by ZF Count:   -zf   zfCount 
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        if (flagLoc( ARGC, ARGV, "-help" ))
           {
            FPR( stderr, "ZF: Extend By Zero Filling;\n" );
            FPR( stderr, "Use only one of the following:\n" );
            FPR( stderr, " -zf   zfCnt   [1] Number of Times to\n" );
            FPR( stderr, "                   Double the Size.\n" );
            FPR( stderr, " -pad  padCnt      Zeros to Add.\n" );
            FPR( stderr, " -size finSize     Final Size.\n" );
            FPR( stderr, "Other Flags:\n" );
            FPR( stderr, " -auto Round Final Size to Power of 2.\n");
            FPR( stderr, "Removing Previous Zero Filling:\n" );
            FPR( stderr, " -inv  Extract Original Time Domain.\n" );
            return( 1 );
           }

        zfCount  = 1;
        padCount = 0;
        zfSize   = 0;
        invFlag  = 0;
        autoFlag = 0;

        if (flagLoc( ARGC, ARGV, "-inv" ))  invFlag  = 1;
        if (flagLoc( ARGC, ARGV, "-auto" )) autoFlag = 1;

        (void) intArgD( ARGC, ARGV, "-zf",   &zfCount );
        (void) intArgD( ARGC, ARGV, "-pad",  &padCount );
        (void) intArgD( ARGC, ARGV, "-size", &zfSize );

        if (invFlag && !flagLoc( ARGC, ARGV, "-zf" )) zfCount = 0;
       }

/***/
/* Initialization and Header:
/*
/*  Inverse Zero Fill (Extraction):
/*     If Zero Fill Count is supplied, undo that amount of zero fill.
/*     If Pad Count is supplied, undo that amount of zero fill.
/*     Otherwise, extract apparent time domain size.
/*  Forward Mode (Zero Filling:
/*     If Explicit Final Size is set, use that.
/*     If Pad Count is set, use that.
/*     Otherwise, double data according to zero fill count.
/*
/*     Round final size to next power of two if needed.
/* 
/*  Find final size after Zero Fill or Extraction.
/*  Update output size.
/*  Update zero frequency location for ppm calibration if needed.
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
        if (invFlag)
           {
            if (zfCount > 0)
               {
                for( i = 0; i < zfCount; i++ ) dataInfo->outSize /= 2;
                zfSize = dataInfo->outSize > 0 ? dataInfo->outSize : 1;
                dataInfo->outSize = zfSize;
               }
            else if (padCount > 0)
               {
                zfSize = dataInfo->size - padCount;
                zfSize = zfSize < 1 ? 1 : zfSize;
                dataInfo->outSize = zfSize;
               } 
            else
               {
                zfSize = dataInfo->timeSize;
                dataInfo->outSize = zfSize;
               }
           }
        else
           {
            if (zfSize)
               {
                dataInfo->outSize = zfSize;
               }
            else if (padCount)
               {
                zfSize = dataInfo->size + padCount;
                dataInfo->outSize = zfSize;
               }
            else
               {
                if (zfCount < 0) zfCount = 0;
                for( i = 0; i < zfCount; i++ ) dataInfo->outSize *= 2;
                zfSize = dataInfo->outSize;
               }

            if (autoFlag)
               {
                zfSize = nextPower2( zfSize );
                dataInfo->outSize = zfSize;
               }
           }

        mid   = getParm( HDR, NDCENTER, CUR_XDIM );
        orig  = getParm( HDR, NDORIG,   CUR_XDIM );
        car   = getParm( HDR, NDCAR,    CUR_XDIM );
        obs   = getParm( HDR, NDOBS,    CUR_XDIM );
        sw    = getParm( HDR, NDSW,     CUR_XDIM );
        ix1   = getParm( HDR, NDX1,     CUR_XDIM );
        ixn   = getParm( HDR, NDXN,     CUR_XDIM );
        izf   = getParm( HDR, NDZF,     CUR_XDIM );

        if (dataInfo->ftFlag)
           {
            mid += dataInfo->outSize - dataInfo->size;
            (void) setParm( HDR, NDCENTER, (float)mid, CUR_XDIM );
           }
        else
           { 
            if (dataInfo->quadState == 1)
               fSize = dataInfo->outSize/2;
            else 
               fSize = dataInfo->outSize;

            if (ix1 || ixn)
               {
                /* Fix this case. */
               }
            else
               {
                mid  = fSize/2 + 1;
                orig = obs*car - sw*(fSize - mid)/fSize;
               }


            (void)setParm( HDR, NDZF,     -(float)dataInfo->outSize, CUR_XDIM );
            (void)setParm( HDR, NDCENTER,  (float)mid,               CUR_XDIM );
            (void)setParm( HDR, NDX1,      (float)ix1,               CUR_XDIM );
            (void)setParm( HDR, NDXN,      (float)ixn,               CUR_XDIM );
            (void)setParm( HDR, NDORIG,    orig,                     CUR_XDIM );
           }

        (void) setParm( HDR, NDSIZE, (float)dataInfo->outSize, CUR_XDIM );

        if (dataInfo->outSize > dataInfo->maxSize)
           {
            dataInfo->maxSize = dataInfo->outSize;
           }
       }

/***/
/* Shutdown:
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
       }


/***/
/* Unknown mode:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* uXY2YX: 2D Transpose Function XY->YX.
/***/

int uXY2YX( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    float *srcPtr1, *srcPtr2, *destPtr;
    int   i, xID, yID, error;

    static int sliceCount = 0, autoFlag, hyperFlag = 0, length = 0, jump = 0;
    static int maxTries, timeOut;

    static float *work = (float *) NULL;

    error = 0;

/***/
/* Process:
/*   Allocate transpose buffer if this is the first slice.
/*   Disperse the 1D vectors into the 2D matrix:
/*     Hypercomplex data gets dispersed as alternating R/I points.
/*     Hypercomplex odd slices go to the real part of result.
/*     Hypercomplex even slices go the the imag part of result.
/*     Normal data dispersed as vector of reals, then vector of imags if any.
/*   Write results if this slice completes a plane.
/***/
      
    if (dataInfo->sliceCode > 0)
       {
        sliceCount++;

        if (dataInfo->sliceCode == 1)
           {
            if (!(work = fltAlloc( "tp", length ))) return( 1 );
           }

        if (hyperFlag)
           {
            if (sliceCount % 2)
               destPtr = work + sliceCount/2;
            else
               destPtr = work + dataInfo->outSize + sliceCount/2 - 1;
               
            srcPtr1 = dataInfo->rdata;
            srcPtr2 = dataInfo->idata;

            for( i = 0; i < dataInfo->size; i++ )
               {
                *destPtr  = *srcPtr1++;
                destPtr  += jump;
                *destPtr  = *srcPtr2++;
                destPtr  += jump;
               }
           }
        else
           {
            destPtr = work + sliceCount - 1;
            srcPtr1 = dataInfo->rdata;

            for( i = 0; i < dataInfo->size; i++ )
               {
                *destPtr = *srcPtr1++;
                destPtr += jump;
               }
       
            if (dataInfo->quadState != 1)
               { 
                destPtr = work + dataInfo->outSize + sliceCount - 1;
                srcPtr1 = dataInfo->idata;
    
                for( i = 0; i < dataInfo->size; i++ )
                   {
                    *destPtr = *srcPtr1++;
                    destPtr += jump;
                   }
               }
           }

        if (sliceCount >= dataInfo->specnum)
           {
            sliceCount = 0;

            error = dataWriteB( dataInfo->outUnit,
                                work, sizeof(float)*length,
                                maxTries, timeOut );

            if (error)
               {
                FPR( stderr, "YTP Error writing plane.\n" );
                return( error );
               }
           }
       }
 
/***/
/* Command line: test for hypercomplex and N-dimensional XY TP. 
/*  Test for blocking I/O paramaters, set only for pipe I/O.
/*  Abort if any data type adjustments are set.
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        work = (float *) NULL;

        if (flagLoc( ARGC, ARGV, "-help" ))
           {
            FPR( stderr, "TP: 2D Plane Transpose. [YTP XY2YX]\n" );
            FPR( stderr, " -hyper       Hypercomplex Transpose.\n" );
            FPR( stderr, " -nohyper     Suppress Hypercomplex Mode.\n" );
            FPR( stderr, " -auto        Choose Mode Automatically.\n" );
            return( 1 );
           }

        hyperFlag = 0;
        autoFlag  = 0;
        maxTries  = -1;
        timeOut   = getDefTimeout();

        (void) intArgD( ARGC, ARGV, "-mt",   &maxTries );
        (void) intArgD( ARGC, ARGV, "-to",   &timeOut );
        (void) logArgD( ARGC, ARGV, "-auto", &autoFlag );

        if (!dataInfo->outName)
           {
            maxTries = 1;
            timeOut  = 0;
           }

        if (DI->typeFlags[XIFLAG] || DI->typeFlags[DTFLAG] ||
              DI->typeFlags[AIFLAG] || DI->typeFlags[ADFLAG] ||
                DI->typeFlags[DIFLAG] || DI->typeFlags[ACFLAG])
           {
            FPR( stderr, "NMRPipe Error: data type adjustment\n" );
            FPR( stderr, "is not allowed during Transpose.\n" );
            return( 1 );
           }
       }

/***/
/* Initialize, Header:
/*   Suppress hypercomplex TP for real-only data.
/*   Show warnings given data type and transpose type.
/*   Indicate that this function is not a vector-wise one.
/*   Adjust transpose state, swap size, specnum.
/*   Initialize static dimensions for dispersal.
/*   Note special quad state handling for magnitude-mode data.
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
        if (dataInfo->quadState == 1) hyperFlag = 0;

        if (autoFlag)
           {
            if (ISQUAD( CUR_XDIM ) && ISQUAD( CUR_YDIM ))
               hyperFlag = 1;
            else
               hyperFlag = 0;
           }

        if (flagLoc( ARGC, ARGV, "-hyper" ))   hyperFlag = 1;
        if (flagLoc( ARGC, ARGV, "-nohyper" )) hyperFlag = 0;

        if (hyperFlag)
           {
            if (!ISQUAD( CUR_XDIM ) || !ISQUAD( CUR_YDIM ))
               {
                FPR( stderr, "\n" );

                FPR( stderr,
                     "NMRPipe Transpose Warning, Function %d:\n",
                     (int)getParm( HDR, FDPIPECOUNT, 0 ));

                FPR( stderr,
                     "%s%s%s%s%s\n",
                     " Hypercomplex Transpose  of  Ordinary Data.\n",
                     " Consider creating hypercomplex data first;\n",
                     " use the -ac and -ad flags.\n",
                     " Consider using an ordinary transpose;\n",
                     " use the -nohyper flag to select this.\n" );
               }
           }
        else
           {
            if (ISQUAD( CUR_XDIM ) && ISQUAD( CUR_YDIM ))
               {
                FPR( stderr, "\n" );

                FPR( stderr,
                     "NMRPipe Transpose Warning, Function %d:\n",
                     (int)getParm( HDR, FDPIPECOUNT, 0 ));

                FPR( stderr,
                     "%s%s%s%s%s\n",
                     " Ordinary Transpose of Hypercomplex Data.\n",
                     " Consider deleting  imaginary data first;\n",
                     " use the -di flag  to delete imaginaries.\n",
                     " Consider using  hypercomplex  transpose;\n",
                     " use the -hyper flag to select this.\n" );
               }
           }

        dataInfo->scaleFlag = 0;
        dataInfo->vecFlag   = 0;
        sliceCount          = 0;
 
        if (dataInfo->transposed)
           (void) setParm( HDR, FDTRANSPOSED, (float)0.0, 0 );
        else
           (void) setParm( HDR, FDTRANSPOSED, (float)1.0, 0 );

        xID = getParm( HDR, FDDIMORDER1, 0 );
        yID = getParm( HDR, FDDIMORDER2, 0 );

        (void) setParm( HDR, FDDIMORDER1, (float)yID, 0 );
        (void) setParm( HDR, FDDIMORDER2, (float)xID, 0 );

        if (hyperFlag)
           {
            (void) setParm( HDR, NDSIZE, (float)dataInfo->specnum/2, CUR_XDIM );
            (void) setParm( HDR, NDSIZE, (float)dataInfo->size*2,    CUR_YDIM );

            dataInfo->outSize = dataInfo->specnum/2;
            jump              = dataInfo->specnum;
           }
        else
           { 
            (void) setParm( HDR, NDSIZE, (float)dataInfo->specnum, CUR_XDIM );
            (void) setParm( HDR, NDSIZE, (float)dataInfo->size,    CUR_YDIM );

            dataInfo->outSize = dataInfo->specnum;
            jump              = dataInfo->specnum*dataInfo->quadState;
           }

        if (FD_MAGNITUDE == (int)getParm( HDR, FD2DPHASE, 0 ))
           {
            xID = getParm( HDR, NDQUADFLAG, CUR_XDIM );
            yID = getParm( HDR, NDQUADFLAG, CUR_YDIM );

            (void) setParm( HDR, NDQUADFLAG, (float)yID, CUR_XDIM );
            (void) setParm( HDR, NDQUADFLAG, (float)xID, CUR_YDIM );
           }

        i = (dataInfo->sliceCount*dataInfo->size)/dataInfo->outSize;

        if (!dataInfo->varCountFlag)
           {
            (void) setParm( HDR, FDSLICECOUNT, (float)i, 0 );
           }

        length = dataInfo->size*dataInfo->specnum*dataInfo->quadState;
       }

/***/
/* Shutdown:
/*   Deallocate the workspace.
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
        (void) deAlloc( "tp", work, sizeof(float)*length );

        work   = (float *) NULL;
        length = 0;
       }

/***/
/* Error:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* uXZY2ZYX: 3D Transpose Function.
/***/

int uXYZ2ZYX( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    float *srcPtr1, *destPtr;
    int   i, xID, zID, error;

    static int length = 0, pntJump, sliceJump;

    static int xOffset, yOffset, zOffset, xOutSize, yOutSize, zOutSize,
               xQuadState, yQuadState, zQuadState;

    static int maxTries, timeOut;

    static float *work = (float *) NULL;

    error = 0;

/***/
/* Process:
/*   Allocate transpose buffer if this is the first slice.
/*   Disperse the 1D vectors into the 3D matrix:
/*     Alternating odd/even planes go to real/imag (left/right) sides of output.
/*     Real/Imag input points are alternated within a column.
/*   Write results if this slice completes the 3D matrix.
/***/
      
    if (dataInfo->sliceCode > 0)
       {

/* Allocate memory. */

        if (dataInfo->sliceCode == 1)
           {
            if (!(work = fltAlloc( "tp", length )))
               {
                FPR( stderr, "XYZ2ZYX Error allocating memory.\n" );
                return( 1 );
               }
           }

/* Move real input vector into a Z-column. */

        srcPtr1 = dataInfo->rdata;

        if (xQuadState == 1 && yQuadState == 2 && zQuadState == 2)
           destPtr = work + (2*yOffset)*sliceJump + xOffset;
        else
           destPtr = work + yOffset*sliceJump + xOffset;

        if (zQuadState == 2 && (zOffset % 2))
           {
            destPtr += dataInfo->zSize/zQuadState;
           }

        for( i = 0; i < dataInfo->size; i++ )
           {
            *destPtr = *srcPtr1++;
            destPtr += pntJump;
           }

/* Move imaginary input vector into a Z-column. */
       
        if (dataInfo->quadState != 1)
           { 
            srcPtr1 = dataInfo->idata;

            if (xQuadState == 1 && yQuadState == 2 && zQuadState == 2)
               destPtr = work + (2*yOffset + 1)*sliceJump + xOffset;
            else if (xQuadState == 1 && yQuadState == 2 && zQuadState == 1)
               destPtr = work + yOffset*sliceJump + xOffset + dataInfo->zSize;
            else
               destPtr = work + yOffset*sliceJump + xOffset + pntJump/2;

            if (zQuadState == 2 && (zOffset % 2))
               {
                destPtr += dataInfo->zSize/zQuadState;
               }

            for( i = 0; i < dataInfo->size; i++ )
               {
                *destPtr = *srcPtr1++;
                destPtr += pntJump;
               }
           }

/* Test for completion of a 2D plane: */
 
        if (++yOffset >= dataInfo->specnum)
           {
            if (zQuadState == 2)
               {
                if (zOffset % 2) xOffset++;
               }
            else
               {
                xOffset++;
               }

            yOffset = 0;
            zOffset++;
           }

        if (xOffset >= dataInfo->zSize/zQuadState) xOffset = 0;

/* Test for completion of a 3D cube: */
 
        if (zOffset >= dataInfo->zSize)
           {
            error = dataWriteB( dataInfo->outUnit,
                                work, sizeof(float)*length,
                                maxTries, timeOut );

            if (error)
               {
                FPR( stderr, "ZTP Error writing cube.\n" );
                return( error );
               }

            xOffset = 0;
            yOffset = 0;
            zOffset = 0;
           }
       }
 
/***/
/* Command line: 
/*  Test for blocking I/O paramaters, set only for pipe I/O.
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        work = (float *) NULL;

        if (flagLoc( ARGC, ARGV, "-help" ))
           {
            FPR( stderr, "ZTP: 3D Matrix Transpose. [XYZ2ZYX]\n" );
            return( 1 );
           }

        maxTries  = -1;
        timeOut   = getDefTimeout();

        (void) intArgD( ARGC, ARGV, "-mt", &maxTries );
        (void) intArgD( ARGC, ARGV, "-to", &timeOut );

        if (!dataInfo->outName)
           {
            maxTries = 1;
            timeOut  = 0;
           }

        if (DI->typeFlags[XIFLAG] || DI->typeFlags[DTFLAG] ||
              DI->typeFlags[AIFLAG] || DI->typeFlags[ADFLAG] ||
                DI->typeFlags[DIFLAG] || DI->typeFlags[ACFLAG])
           {
            FPR( stderr, "NMRPipe Error: data type adjustment\n" );
            FPR( stderr, "is not allowed during Transpose.\n" );
            return( 1 );
           }
       }

/***/
/* Initialize, Header:
/*   Swap sizes and quad states.
/*   Initialize static dimensions for dispersal.
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
        dataInfo->scaleFlag = 0;
        dataInfo->vecFlag   = 0;
        dataInfo->verbInc   = dataInfo->specnum;

        xQuadState = getQuad( HDR, NDQUADFLAG, CUR_XDIM );
        yQuadState = getQuad( HDR, NDQUADFLAG, CUR_YDIM );
        zQuadState = getQuad( HDR, NDQUADFLAG, CUR_ZDIM );

        dataInfo->outSize = dataInfo->zSize/zQuadState;
        xOutSize          = dataInfo->zSize/zQuadState;
        yOutSize          = dataInfo->specnum;
        zOutSize          = dataInfo->size*xQuadState;

        if (xQuadState == 1 && yQuadState == 2 && zQuadState == 2) yOutSize*=2;
        if (xQuadState == 2 && yQuadState == 2 && zQuadState == 1) yOutSize/=2;

        (void) setParm( HDR, NDSIZE, (float)xOutSize, CUR_XDIM );
        (void) setParm( HDR, NDSIZE, (float)yOutSize, CUR_YDIM );
        (void) setParm( HDR, NDSIZE, (float)zOutSize, CUR_ZDIM );

        xID = getParm( HDR, FDDIMORDER1, 0 );
        zID = getParm( HDR, FDDIMORDER3, 0 );

        (void) setParm( HDR, FDDIMORDER1, (float)zID, 0 );
        (void) setParm( HDR, FDDIMORDER3, (float)xID, 0 );

        if (yQuadState == 1 && zQuadState == 1)
           {
            (void) setParm( HDR, FDQUADFLAG, (float)1.0, 0 );
            dataInfo->outQuadState = 1;
           }
        else
           {
            (void) setParm( HDR, FDQUADFLAG, (float)0.0, 0 );
            dataInfo->outQuadState = 2;
           }

        if (xQuadState == 1 && yQuadState == 2 && zQuadState == 1)
           {
            sliceJump = dataInfo->zSize*dataInfo->quadState;
            pntJump   = dataInfo->zSize*dataInfo->specnum*dataInfo->quadState;
           }
        else
           {
            sliceJump = dataInfo->zSize;
            pntJump   = dataInfo->zSize*dataInfo->specnum*dataInfo->quadState;
           }
 
        xOffset   = 0;
        yOffset   = 0;
        zOffset   = 0;

        length = 
          dataInfo->size*dataInfo->specnum*dataInfo->zSize*dataInfo->quadState;
       }

/***/
/* Shutdown:
/*   Deallocate the workspace.
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
        (void) deAlloc( "tp", work, sizeof(float)*length );

        work   = (float *) NULL;
        length = 0;
       }

/***/
/* Error:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* uFDM: FDM spectral reconstruction.
/***/

int uFDM( dataInfo )

   struct ProcDataInfo *dataInfo;
{
    static struct FDMInfo fInfo;
    static int    xSize = 0,   ySize = 0;
    static float  xSW   = 0.0, ySW   = 0.0;

    int    i, yLoc, procFlag, outSize, tSize, error;

    error = 0;

/***/
/* Process, Write:
/***/

    if (dataInfo->sliceCode > 0)
       {
        switch( fInfo.dimCount )
           {
            case 1:
               (void) fdmPut1D( &fInfo.fInfo1D, DI->rdata, DI->idata );
               procFlag = 1;
               break;

            case 2:
               yLoc     = DI->locList[YLOC];
               procFlag = yLoc == DI->specnum;

               (void) fdmPut2D( &fInfo.fInfo2D, DI->rdata, DI->idata, yLoc );
               break;

            default: 
               return( 1 );
               break;
           }

        if (!procFlag) return( 0 );

        if (error = fdmNDC( &fInfo )) return( error );

        switch( fInfo.dimCount )
           {
            case 1:
               (void) fdmGet1D( &fInfo.fInfo1D, DI->rdata, DI->idata );
               error = putSlice( dataInfo );
               break;

            case 2:
               for( i = 1; i <= fInfo.fInfo2D.NSp2; i++ )
                  {
                   (void) fdmGet2D( &fInfo.fInfo2D, DI->rdata, DI->idata, i );
                   error = putSlice( dataInfo );
                  }
 
               break;

            default:
               return( 1 );
               break;
           }
       }

/***/
/* Command Line:
/***/

    else if (dataInfo->sliceCode == FN_PARAMS)
       {
        (void) fdmNDCNull( &fInfo );

        error = fdmNDCArgs( ARGC, ARGV, &fInfo );
       }

/***/
/* Initialization and Header:
/***/

    else if (dataInfo->sliceCode == FN_INIT)
       {
        xSW   = getParm( HDR, NDSW, CUR_XDIM );
        ySW   = getParm( HDR, NDSW, CUR_YDIM );
        xSize = getParm( HDR, NDAPOD, CUR_XDIM );
        ySize = getParm( HDR, NDAPOD, CUR_YDIM );

        if (error = fdmNDCInit( xSize, ySize, xSW, ySW, &fInfo ))
           {
            return( error );
           }

        if (fInfo.dimCount == 1)
           {
            outSize = fInfo.fInfo1D.NSp;
            tSize   = outSize/2;

            dataInfo->outSize      = outSize;
            dataInfo->timeSize     = tSize, 
            dataInfo->outQuadState = 2;

            (void) setParm( HDR, NDFTFLAG,   (float)1.0,     CUR_XDIM );
            (void) setParm( HDR, NDFTSIZE,   (float)outSize, CUR_XDIM );
            (void) setParm( HDR, NDSIZE,     (float)outSize, CUR_XDIM );
            (void) setParm( HDR, NDAPOD,     (float)tSize,   CUR_XDIM );
            (void) setParm( HDR, NDTDSIZE,   (float)tSize,   CUR_XDIM );
            (void) setParm( HDR, FDQUADFLAG, (float)0.0,     CUR_XDIM );
            (void) setParm( HDR, NDQUADFLAG, (float)0.0,     CUR_XDIM );
           }
        else if (fInfo.dimCount == 2)
           {
            outSize = fInfo.fInfo2D.NSp1;
            tSize   = outSize/2;

            dataInfo->outSize      = outSize;
            dataInfo->timeSize     = tSize, 
            dataInfo->outQuadState = 2;

            (void) setParm( HDR, NDFTFLAG,   (float)1.0,     CUR_XDIM );
            (void) setParm( HDR, NDFTSIZE,   (float)outSize, CUR_XDIM );
            (void) setParm( HDR, NDSIZE,     (float)outSize, CUR_XDIM );
            (void) setParm( HDR, NDAPOD,     (float)tSize,   CUR_XDIM );
            (void) setParm( HDR, NDTDSIZE,   (float)tSize,   CUR_XDIM );
            (void) setParm( HDR, FDQUADFLAG, (float)0.0,     CUR_XDIM );
            (void) setParm( HDR, NDQUADFLAG, (float)0.0,     CUR_XDIM );

            outSize = fInfo.fInfo2D.NSp2;
            tSize   = outSize/2;

            (void) setParm( HDR, NDFTFLAG,   (float)1.0,     CUR_YDIM );
            (void) setParm( HDR, NDFTSIZE,   (float)outSize, CUR_YDIM );
            (void) setParm( HDR, NDSIZE,     (float)outSize, CUR_YDIM );
            (void) setParm( HDR, NDAPOD,     (float)tSize,   CUR_YDIM );
            (void) setParm( HDR, NDTDSIZE,   (float)tSize,   CUR_YDIM );
            (void) setParm( HDR, NDQUADFLAG, (float)1.0,     CUR_YDIM );
           }

        dataInfo->ftFlag = 1;

        if (DI->outSize > DI->maxSize) DI->maxSize = DI->outSize;

        i = (dataInfo->sliceCount*dataInfo->size)/dataInfo->outSize;

        if (!dataInfo->varCountFlag)
           {
            (void) setParm( HDR, FDSLICECOUNT, (float)i, 0 );
           }
       }

/***/
/* Shutdown:
/***/

    else if (dataInfo->sliceCode == FN_SHUTDOWN)
       {
        (void) fdmNDCFree( &fInfo );
       }

/***/
/* Unknown mode:
/***/

    else
       {
        error = 1;
       }

    return( error );
}

/***/
/* rSet: set elements in rvec to rval.
/***/

static int rSet( rvec, length, rval )

   float *rvec, rval;
   int   length;
{
    while( length-- ) *rvec++ = rval;
    return( 0 );
}

/***/
/* rMove: move elements from vec1 to vec2.
/***/

static int rMove( vec1, vec2, length )

   float *vec1, *vec2;
   int length;
{
    while( length-- ) *vec2++ = *vec1++;
    return( 0 );
}

/***/
/* vNeg: negate data in vec;
/**/

int vNeg( vec, length )

   float *vec;
   int   length;
{
    while( length-- )
       {
        *vec = -(*vec);
        vec++;
       }

    return( 0 );
}

/***/
/* vAlt: sign-alternate data in vec (+ -).
/***/

int vAlt( vec, length )

   float *vec;
   int   length;
{
    length /= 2;
    vec++;

    while( length-- )
       {
        *vec = -(*vec);
        vec += 2;
       }

    return( 0 );
}

int slowFFT( rdata, idata, bR, bI, size )

   float *rdata, *idata, *bR, *bI;
   int   size;
{
   float vR, vI, f;
   int   k, kk, n;

   (void) rSet( bR, size, 0.0 );
   (void) rSet( bI, size, 0.0 );

   kk = size/2;

   for( k = 0; k < size; k++ )
      {
       if (kk < 0) kk = size - 1;

       for( n = 0; n < size; n++ )
          {
           f = -2.0*PI*k*n/size;

           vR = cos( f );
           vI = sin( f );

           bR[kk] += rdata[n]*vR - idata[n]*vI;
           bI[kk] += rdata[n]*vI + idata[n]*vR;
          } 

       kk--;
      }

   (void) rMove( bR, rdata, size );
   (void) rMove( bI, idata, size );

   return( 0 );
} 

int slowIFT( rdata, idata, bR, bI, size )

   float *rdata, *idata, *bR, *bI;
   int   size;
{
   float vR, vI, f;
   int   k, kk, n;

   (void) rSet( bR, size, 0.0 );
   (void) rSet( bI, size, 0.0 );

   kk = 0;

   for( k = 0; k < size; k++ )
      {
       for( n = 0; n < size; n++ )
          {
           f = -2.0*PI*k*n/size;

           vR = cos( f )/size;
           vI = sin( f )/size;

           bR[kk] += rdata[n]*vR - idata[n]*vI;
           bI[kk] += rdata[n]*vI + idata[n]*vR;
          } 

       kk++;
      }

   (void) rMove( bR, rdata, size );
   (void) rMove( bI, idata, size );

   return( 0 );
}

/***/
/* nextPower2: returns power of 2 nearest to N.
/***/

static int nextPower2( n )

   int n;
{
    int m;

    m = 1;

    while( n > m ) m *= 2;

    return( m );
}

/***/
/* minMax2: return highest and lowest values in vec without first
/*          initializing maxVal and minVal.
/***/

int minMax2( vec, length, minVal, maxVal )

   float *vec, *minVal, *maxVal;
   int    length;

{
    while( length-- )
       {
        *minVal = (*vec < *minVal) ? *vec : *minVal;
        *maxVal = (*vec > *maxVal) ? *vec : *maxVal;
        vec++;
       }

    return( 0 );
}

/***/
/* Bottom.
/***/