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

/***/
/* FDM: modules to call FDM methods.
/*
/* References:
/*   1. H. Hu et. al., J. Magn. Reson. 134, 76-87 (1998).
/*   2. J. Chen et. al., J. Chem. Phys. 112, 4429-4437 (2000).
/*   3. V. A. Mandelshtam, Prog. NMR Spect. 38, 159-196 (2001).
/*   4. QZ: http://www.netlib.org/toms/535
/***/

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

#include "memory.h"
#include "cmndargs.h"
#include "specunit.h"
#include "fdatap.h"
#include "prec.h"
#include "fdm.h"

#define FCOMPLEX_ALLOC( N ) \
 (FComplex *)voidAlloc( "fdm", sizeof(FComplex)*(N) )

#define FDOUBLE_ALLOC( N ) \
 (FDouble *)voidAlloc( "fdm", sizeof(FDouble)*(N) )

#define FCOMPLEX_DEALLOC( PTR, N ) \
 deAlloc( "fdm", (PTR), sizeof(FComplex)*(N) )

#define FDOUBLE_DEALLOC( PTR, N ) \
 deAlloc( "fdm", (PTR), sizeof(FDouble)*(N) )

#define ICLIP( M, LO, HI ) if (M < LO) M = LO; if (M > HI) M = HI
#define RSWAP( RA, RB )  rtemp = RA; RA = RB; RB = rtemp

#define FPR (void)fprintf

#ifndef LINT_CHECK
#ifndef USE_FDM
void fdm1dmw_() {}
void fdm2dmw_() {}
#endif
#endif

/***/
/* fdmNDCNull: null-initialize FDM structure.
/***/

int fdmNDCNull( fInfo )

   struct FDMInfo *fInfo;
{
    fInfo->dimCount = 1;

    (void) fdm1DCNull( &fInfo->fInfo1D );
    (void) fdm2DCNull( &fInfo->fInfo2D );

    return( 0 );
}

/***/
/* fdmNDCArgs: extract command-line arguments for FDM.
/***/

int fdmNDCArgs( argc, argv, fdata, fInfo )

   struct FDMInfo *fInfo;

   int  argc;
   char **argv;
   float fdata[FDATASIZE];
{
    int error;

    error = fdm1DCArgs( argc, argv, fdata, &fInfo->fInfo1D );
    error = fdm2DCArgs( argc, argv, fdata, &fInfo->fInfo2D );

    (void) intArgD( argc, argv, "-ndim", &fInfo->dimCount );

    if (flagLoc( argc, argv, "-help" ))
       {
        FPR( stderr, " -ndim   [1]        FDM Dimension Count.\n" );
        return( 1 );
       }

    return( error );
}

/***/
/* fdmNDCInit: initialize FDM Engine.
/***/

int fdmNDCInit( xSize, ySize, xSW, ySW, fInfo )

   struct FDMInfo *fInfo;
   int    xSize, ySize;
   float  xSW, ySW;
{
    int error;

    error = 0;

    switch( fInfo->dimCount )
       {
        case 1:
           error = fdm1DCInit( xSize, xSW, &fInfo->fInfo1D );
           break;

        case 2:
           error = fdm2DCInit( xSize, ySize, xSW, ySW, &fInfo->fInfo2D );
           break;

        default:
           FPR( stderr, 
                "FDM Error: Bad Dimension Count %d.\n", 
                fInfo->dimCount );

           return( 1 );
           break;
       }

    return( error );
}

/***/
/* fdmNDC: apply FDM to the current data.
/***/

int fdmNDC( fInfo )

   struct FDMInfo *fInfo;
{
    int error;

    error = 0;

    switch( fInfo->dimCount )
       {
        case 1:
           error = fdm1DC( &fInfo->fInfo1D );
           break;

        case 2:
           error = fdm2DC( &fInfo->fInfo2D );
           break;

        default:
           FPR( stderr,
                "FDM Error: Bad Dimension Count %d.\n",
                fInfo->dimCount );

           return( 1 );
           break;
       }

    return( error );
}

/***/
/* fdm1DCInit: initialize 1D FDM Engine.
/***/

int fdm1DCInit( size, sw, fInfo )

   struct FDM1DInfo *fInfo;
   float  sw;
   int    size;
{
    int   Nb0Max, maxSize, error;
    float rtemp, range;

    error = 0;

    if (sw == 0.0) sw = 1.0;

    ICLIP( fInfo->wmin, -sw/2.0, sw/2.0 );
    ICLIP( fInfo->wmax, -sw/2.0, sw/2.0 );

    maxSize = size - 1;

    if (!(maxSize % 2)) maxSize--;

    if (!fInfo->NSig) fInfo->NSig = maxSize;
    if (!fInfo->NSp)  fInfo->NSp  = 2*size;

    if (!(fInfo->NSig % 2)) fInfo->NSig--;

    ICLIP( fInfo->NSig, 1, maxSize );

    if (fInfo->wmin > fInfo->wmax) {RSWAP( fInfo->wmin, fInfo->wmax );}

    range  = (fInfo->wmax - fInfo->wmin)/sw;
    Nb0Max = 1.0e-5 + range*((fInfo->NSig+1)/2)*fInfo->rho;
 
    fInfo->Nb0 = Nb0Max;

/***/
/*  Setup parameters for multi-window calculation if desired
/***/
    
    if (fInfo->mwFlag)
      {
        fInfo->Nb0  = 100;

        ICLIP( fInfo->Nb0, 2, Nb0Max );

        fInfo->Nbc  = 0;	    
	fInfo->dW   = fInfo->Nb0*(sw/(fInfo->rho*((fInfo->NSig-1)/2 + 1)));
	fInfo->NWin = (fInfo->wmax - fInfo->wmin)/fInfo->dW;
	fInfo->NWin = fInfo->NWin*2 - 1;

	if (fInfo->NWin == 0) fInfo->NWin = 1;

	fInfo->dW   = (fInfo->wmax - fInfo->wmin)/((fInfo->NWin+1)/2);
	fInfo->Nb0  = (fInfo->dW/sw)*((fInfo->NSig-1)/2 + 1)*fInfo->rho;
	fInfo->rho  = fInfo->Nb0*sw/(fInfo->dW*((fInfo->NSig-1)/2 + 1));
      }
    else
      {
	fInfo->NWin = 1;
	fInfo->dW   = (fInfo->wmax - fInfo->wmin);
      }

/***/
/*  End multi-window adjustment
/***/

    ICLIP( fInfo->Nb0, 2, (fInfo->NSig-1)/2 + 1 );

    fInfo->Nb    = fInfo->Nb0 + fInfo->Nbc;
    fInfo->delta = 2.0*PI/sw;
   
    if (fInfo->gamma < 0.0) fInfo->gamma = sw/fInfo->NSp;

    error = fdm1DCAlloc( fInfo );

    return( error );
}

/***/
/* fdm2DCInit: initialize 2D FDM Engine.
/***/

int fdm2DCInit( xSize, ySize, xSW, ySW, fInfo )

   struct FDM2DInfo *fInfo;
   int    xSize, ySize;
   float  xSW, ySW;
{
    int    NSig1H, NSig2H, NWin1H, NWin2H, Nb1Max, Nb2Max;
    int    maxSizeX, maxSizeY, error;
    double rangeX, rangeY;
    float  rtemp;

    error = 0;

    if (xSW == 0.0) xSW = 1.0;
    if (ySW == 0.0) ySW = 1.0;

    ICLIP( fInfo->wmin1, -xSW/2.0, xSW/2.0 );
    ICLIP( fInfo->wmax1, -xSW/2.0, xSW/2.0 );

    ICLIP( fInfo->wmin2, -ySW/2.0, ySW/2.0 );
    ICLIP( fInfo->wmax2, -ySW/2.0, ySW/2.0 );

    maxSizeX = xSize - 1;
    maxSizeY = ySize - 1;

    if (!(maxSizeX % 2)) maxSizeX--;
    if (!(maxSizeY % 2)) maxSizeY--;

    if (!fInfo->NSig1) fInfo->NSig1 = maxSizeX;
    if (!fInfo->NSig2) fInfo->NSig2 = maxSizeY;

    ICLIP( fInfo->NSig1, 1, maxSizeX );
    ICLIP( fInfo->NSig2, 1, maxSizeY );

    if (!fInfo->NSp1)  fInfo->NSp1  = 2*xSize;
    if (!fInfo->NSp2)  fInfo->NSp2  = 2*ySize;

    if (!(fInfo->NSig1 % 2)) fInfo->NSig1--;
    if (!(fInfo->NSig2 % 2)) fInfo->NSig2--;

    if (fInfo->wmin1 > fInfo->wmax1) {RSWAP( fInfo->wmin1, fInfo->wmax1 );}
    if (fInfo->wmin2 > fInfo->wmax2) {RSWAP( fInfo->wmin2, fInfo->wmax2 );}

    NSig1H = (fInfo->NSig1 + 1)/2;
    NSig2H = (fInfo->NSig2 + 1)/2;

    rangeX = (fInfo->wmax1 - fInfo->wmin1)/xSW;
    rangeY = (fInfo->wmax2 - fInfo->wmin2)/ySW;

    Nb1Max = 1.0e-5 + rangeX*fInfo->rho1*(float)NSig1H;
    Nb2Max = 1.0e-5 + rangeY*fInfo->rho2*(float)NSig2H;

    fInfo->Nb1 = Nb1Max;
    fInfo->Nb2 = Nb2Max;

/***/
/*  Setup parameters for multi-window calculation if desired
/***/
    
    if (fInfo->mwxFlag)
       {
        fInfo->Nb1   = 10;

        ICLIP( fInfo->Nb1, 1, Nb1Max );

	fInfo->dW1   = fInfo->Nb1*(xSW/(fInfo->rho1*NSig1H));
	fInfo->NWin1 = (fInfo->wmax1 - fInfo->wmin1)/fInfo->dW1;
	fInfo->NWin1 = fInfo->NWin1*2 - 1;

	if (fInfo->NWin1 == 0) fInfo->NWin1 = 1;

        NWin1H = (fInfo->NWin1 + 1)/2;

	fInfo->dW1   = (fInfo->wmax1 - fInfo->wmin1)/NWin1H;
	fInfo->Nb1   = (fInfo->dW1/xSW)*NSig1H*fInfo->rho1;
	fInfo->rho1  = fInfo->Nb1*xSW/(fInfo->dW1*NSig1H);
       }
    else
       {
	fInfo->NWin1  = 1;
	fInfo->dW1    = fInfo->wmax1 - fInfo->wmin1;
       }

    if (fInfo->singFlag)
       {
        fInfo->Nb1   = 1;
	fInfo->dW1   = fInfo->Nb1*(xSW/(fInfo->rho1*NSig1H));
	fInfo->NWin1 = (fInfo->wmax1 - fInfo->wmin1)/fInfo->dW1;
	fInfo->NWin1 = fInfo->NWin1*2 - 1;

	if (fInfo->NWin1 == 0) fInfo->NWin1 = 1;

	fInfo->dW1   = (fInfo->wmax1 - fInfo->wmin1)/NWin1H;
	fInfo->Nb1   = (fInfo->dW1/xSW)*NSig1H*fInfo->rho1;
	fInfo->rho1  = fInfo->Nb1*xSW/(fInfo->dW1*NSig1H);
       }

    if (fInfo->mwyFlag)
       {
        fInfo->Nb2   = 10;

        ICLIP( fInfo->Nb2, 1, Nb2Max );

	fInfo->dW2   = fInfo->Nb2*(ySW/(fInfo->rho2*NSig2H));
	fInfo->NWin2 = (fInfo->wmax2 - fInfo->wmin2)/fInfo->dW2;
	fInfo->NWin2 = fInfo->NWin2*2 - 1;

	if (fInfo->NWin2 == 0) fInfo->NWin2 = 1;

        NWin2H = (fInfo->NWin2 + 1)/2;

	fInfo->dW2   = (fInfo->wmax2 - fInfo->wmin2)/NWin2H;
	fInfo->Nb2   = (fInfo->dW2/ySW)*NSig2H*fInfo->rho2;
	fInfo->rho2  = fInfo->Nb2*ySW/(fInfo->dW2*NSig2H);
      }
    else
       {
	fInfo->NWin2 = 1;
	fInfo->dW2   = fInfo->wmax2 - fInfo->wmin2;
       }

/***/
/*  End multi-window adjustment
/***/

    fInfo->Nb     = fInfo->Nb1*fInfo->Nb2;

    fInfo->delta1 = 2.0*PI/xSW;
    fInfo->delta2 = 2.0*PI/ySW;
   
    if (fInfo->gamma1 < 0.0) 
       {
        fInfo->gamma1 = xSW/fInfo->NSp1;
       }

    if (fInfo->gamma2 < 0.0) 
       {
        fInfo->gamma2 = ySW/fInfo->NSp2;
       }

    error = fdm2DCAlloc( fInfo );

    return( error );
}
    
/***/
/* fdm1DCArgs: extract command-line arguments for 1D FDM.
/***/

int fdm1DCArgs( argc, argv, fdata, fInfo )

   struct FDM1DInfo *fInfo;

   int   argc;
   char  **argv;
   float fdata[FDATASIZE];
{
    float rtemp, swx;
    int   xSize;

    if (flagLoc( argc, argv, "-help" ))
       {
        FPR( stderr, "1D FDM:\n" );
        FPR( stderr, " -tsize [TDSIZE]   Number of Input Points to Use.\n"  );
        FPR( stderr, " -fsize [2*TDSIZE] Desired Final Points.\n" );
        FPR( stderr, " -nc    [0]        Number of Coarse Basis Points.\n"  );
        FPR( stderr, " -x1    [-SW/2]    Start of Desired Freq Window.\n" );
        FPR( stderr, " -xn    [+SW/2]    End of Desired Freq Window.\n" );
        FPR( stderr, " -rho   [1.0]      Basis Function Density.\n" );
        FPR( stderr, " -gamma [Auto]     Smoothing Factor.\n" );
	FPR( stderr, " -mw               Multi-Window Mode ON.\n" );
        return( 1 );
       }

    (void) intArgD( argc, argv, "-tsize", &fInfo->NSig );
    (void) intArgD( argc, argv, "-fsize", &fInfo->NSp );
    (void) intArgD( argc, argv, "-nc",    &fInfo->Nbc );

    (void) dblArgD( argc, argv, "-rho",   &fInfo->rho );
    (void) dblArgD( argc, argv, "-gamma", &fInfo->gamma );

    if (flagLoc( argc, argv, "-mw" ))   fInfo->mwFlag = 1;
    if (flagLoc( argc, argv, "-nomw" )) fInfo->mwFlag = 0;

    if (fInfo->Nbc) fInfo->msFlag = 1;

    xSize  = getParm( fdata, NDSIZE, CUR_XDIM );
    swx    = getParm( fdata, NDSW,   CUR_XDIM );

    if (xSize <= 1) return( 1 );
            
    if (isInterleaved( fdata, CUR_XDIM )) xSize /= 2;

    fInfo->sx1 = 1.0;
    fInfo->sxn = xSize;

    (void) specFltArgD( argc, argv, "-x1", &fInfo->sx1, fdata, CUR_XDIM );
    (void) specFltArgD( argc, argv, "-xn", &fInfo->sxn, fdata, CUR_XDIM );

    if (fInfo->sx1 > fInfo->sxn) {RSWAP( fInfo->sx1, fInfo->sxn );}

    fInfo->wmin = -swx/2.0 + swx*(fInfo->sx1 - 1)/(xSize - 1);
    fInfo->wmax = -swx/2.0 + swx*(fInfo->sxn - 1)/(xSize - 1);

    return( 0 );
}

/***/
/* fdm2DCArgs: extract command-line arguments for 2D FDM.
/***/

int fdm2DCArgs( argc, argv, fdata, fInfo )

   struct FDM2DInfo *fInfo;

   int   argc;
   char  **argv;
   float fdata[FDATASIZE];
{
    float rtemp, swx, swy;
    int   xSize, ySize;
 
    if (flagLoc( argc, argv, "-help" ))
       {
        FPR( stderr, "2D FDM:\n" );
        FPR( stderr, " -xtsize [TDSIZE]   Number of Input Points, X.\n" );
        FPR( stderr, " -ytsize [TDSIZE]   Number of Input Points, Y.\n" );
        FPR( stderr, " -xfsize [2*TDSIZE] Desired Final Points, X.\n" );
        FPR( stderr, " -yfsize [2*TDSIZE] Desired Final Points, Y.\n" );
        FPR( stderr, " -x1     [1]        Start of Desired Freq Window, X.\n" );
        FPR( stderr, " -xn     [SIZE]     End of Desired Freq Window, X.\n" );
        FPR( stderr, " -y1     [1]        Start of Desired Freq Window, Y.\n" );
        FPR( stderr, " -yn     [SIZE]     End of Desired Freq Window, Y.\n" );
        FPR( stderr, " -xgamma [Auto]     Smoothing Factor, X.\n" );
        FPR( stderr, " -ygamma [Auto]     Smoothing Factor, Y.\n" );
        FPR( stderr, " -xrho   [1.1]      Basis Function Density, X.\n" );
        FPR( stderr, " -yrho   [1.1]      Basis Function Density, Y.\n" );
        FPR( stderr, " -q2     [1.0e-4]   Regularization Parameter.\n" );
	FPR( stderr, " -mwx               Multi-Window Mode ON, X.\n" );
	FPR( stderr, " -mwy               Multi-Window Mode ON, Y.\n" );
	FPR( stderr, " -sing              Single-Basis Mode ON.\n" );
        FPR( stderr, " -ct                Constant-Time Mode ON.\n" );
        FPR( stderr, " -mode   [1]        1=Gauss 2=Lorentz, 3=Pseudo.\n" );
        return( 1 );
       }

    (void) intArgD( argc, argv, "-xtsize", &fInfo->NSig1 );
    (void) intArgD( argc, argv, "-ytsize", &fInfo->NSig2 );
    (void) intArgD( argc, argv, "-xfsize", &fInfo->NSp1 );
    (void) intArgD( argc, argv, "-yfsize", &fInfo->NSp2 );
    (void) intArgD( argc, argv, "-mode",   &fInfo->iModel );

    (void) dblArgD( argc, argv, "-xrho",   &fInfo->rho1 );
    (void) dblArgD( argc, argv, "-yrho",   &fInfo->rho2 );
    (void) dblArgD( argc, argv, "-xgamma", &fInfo->gamma1 );
    (void) dblArgD( argc, argv, "-ygamma", &fInfo->gamma2 );

    (void) dblArgD( argc, argv, "-q2",     &fInfo->q2 );

    if (flagLoc( argc, argv, "-mwx" ))   fInfo->mwxFlag = 1;
    if (flagLoc( argc, argv, "-nomwx" )) fInfo->mwxFlag = 0;

    if (flagLoc( argc, argv, "-mwy" ))   fInfo->mwyFlag = 1;
    if (flagLoc( argc, argv, "-nomwy" )) fInfo->mwyFlag = 0;

    if (flagLoc( argc, argv, "-sing" ))   fInfo->singFlag = 1;
    if (flagLoc( argc, argv, "-nosing" )) fInfo->singFlag = 0;

    if (flagLoc( argc, argv, "-ct" ))   fInfo->ctFlag = 1;
    if (flagLoc( argc, argv, "-noct" )) fInfo->ctFlag = 0;

    xSize  = getParm( fdata, NDSIZE, CUR_XDIM );
    ySize  = getParm( fdata, NDSIZE, CUR_YDIM );
    swx    = getParm( fdata, NDSW,   CUR_XDIM );
    swy    = getParm( fdata, NDSW,   CUR_YDIM );
            
    if (xSize <= 1 || ySize <= 1) return( 1 );

    if (isInterleaved( fdata, CUR_XDIM )) xSize /= 2;
    if (isInterleaved( fdata, CUR_YDIM )) ySize /= 2;

    fInfo->sx1 = 1.0;
    fInfo->sxn = xSize;
    fInfo->sy1 = 1.0;
    fInfo->syn = ySize;

    (void) specFltArgD( argc, argv, "-x1", &fInfo->sx1, fdata, CUR_XDIM );
    (void) specFltArgD( argc, argv, "-xn", &fInfo->sxn, fdata, CUR_XDIM );
    (void) specFltArgD( argc, argv, "-y1", &fInfo->sy1, fdata, CUR_YDIM );
    (void) specFltArgD( argc, argv, "-yn", &fInfo->syn, fdata, CUR_YDIM );

    if (fInfo->sx1 > fInfo->sxn) {RSWAP( fInfo->sx1, fInfo->sxn );}
    if (fInfo->sy1 > fInfo->syn) {RSWAP( fInfo->sy1, fInfo->syn );}

    fInfo->wmin1 = -swx/2.0 + swx*(fInfo->sx1 - 1)/(xSize - 1);
    fInfo->wmax1 = -swx/2.0 + swx*(fInfo->sxn - 1)/(xSize - 1);
    fInfo->wmin2 = -swy/2.0 + swy*(fInfo->sy1 - 1)/(ySize - 1);
    fInfo->wmax2 = -swy/2.0 + swy*(fInfo->syn - 1)/(ySize - 1);
    
    return( 0 );
}

/***/
/* fdmPut1D: move 1D time-domain input vector into FDM workspace.
/***/

int fdmPut1D( fInfo, rdata, idata )

   float *rdata, *idata;
   struct FDM1DInfo *fInfo;
{
    FComplex *zPtr;
    int      i;

    zPtr = fInfo->td;

    for( i = 0; i <= fInfo->NSig; i++ )
       {
        zPtr->rdata = rdata[i];
        zPtr->idata = idata[i];
        zPtr++;
       }

    return( 0 );
}

/***/
/* fdmPut2D: move 1D time-domain input vector into 2D FDM workspace.
/***/

int fdmPut2D( fInfo, rdata, idata, yLoc )

   float *rdata, *idata;
   struct FDM2DInfo *fInfo;
   int    yLoc;
{
    FComplex *zPtr;
    int      i;

    if (yLoc > fInfo->NSig2 + 1) return( 0 );

    zPtr = fInfo->td + (fInfo->NSig1 + 1)*(yLoc - 1);

    for( i = 0; i <= fInfo->NSig1; i++ )
       {
        zPtr->rdata = rdata[i];
        zPtr->idata = idata[i];
        zPtr++;
       }

    return( 0 );
}

/***/
/* fdmGet1D: get 1D freq-domain output vector from FDM workspace.
/***/

int fdmGet1D( fInfo, rdata, idata )

   float *rdata, *idata;
   struct FDM1DInfo *fInfo;
{
    FComplex *zPtr;
    int      i;

    zPtr = fInfo->riSpec;

    for( i = 0; i < fInfo->NSp; i++ )
       {
        rdata[i] = -zPtr->idata;
        idata[i] = -zPtr->rdata;

        zPtr++;
       }

    return( 0 );
}

/***/
/* fdmGet2D: get 1D freq-domain output vector from 2D FDM workspace.
/***/

int fdmGet2D( fInfo, rdata, idata, yLoc )

   float *rdata, *idata;
   struct FDM2DInfo *fInfo;
   int    yLoc;
{
    FComplex *zPtr;
    int      i;

    zPtr = fInfo->riSpec + fInfo->NSp1*(yLoc - 1);

    for( i = 0; i < fInfo->NSp1; i++ )
       {
        rdata[i] = -zPtr->idata;
        idata[i] = -zPtr->rdata;

        zPtr++;
       }

    return( 0 );
}

/***/
/* fdm1DC: invoke 1D FDM on the stored input data. 
/***/

int fdm1DC( fInfo )

   struct FDM1DInfo *fInfo;
{
    int    NSig, NSp, Nb0, Nbc, Nb, NWin, msFlag;
    double delta, wmin, wmax, dW, gamma;
    int    error;

    error  = 0;

    NSig   = fInfo->NSig;
    NSp    = fInfo->NSp;
    delta  = fInfo->delta;
    Nb0    = fInfo->Nb0;
    Nbc    = fInfo->Nbc;
    Nb     = fInfo->Nb;
    NWin   = fInfo->NWin;
    wmin   = fInfo->wmin;
    wmax   = fInfo->wmax;
    dW     = fInfo->dW;
    msFlag = fInfo->msFlag;
    gamma  = fInfo->gamma;

    (void) fdm1dmw_( fInfo->td, 
                     &NSig, 
		     &NSp,
                     &delta, 
		     &Nb0, 
		     &Nbc, 
		     &Nb, 
		     &NWin,
		     &wmin, 
		     &wmax,
		     &dW,
		     &msFlag,
		     &gamma,
		     fInfo->riSpec,
		     fInfo->wk, 
		     fInfo->dk,
		     fInfo->U, 
		     fInfo->f, 
		     fInfo->g, 
		     fInfo->diag, 
		     fInfo->z, 
		     fInfo->coefW, 
		     fInfo->zz, 
		     fInfo->U0, 
		     fInfo->beta,
		     fInfo->wSpec,
		     &error );

   return( error );
}

/***/
/* fdm2DC: invoke 2D FDM on the stored input data.
/***/

int fdm2DC( fInfo )

   struct FDM2DInfo *fInfo;
{
    double delta1, delta2, wmin1, wmax1, wmin2, wmax2, dW1, dW2, gamma1, gamma2, q2, s;
    int    NSig1, NSig2, Nb1, Nb2, Nb, NSp1, NSp2, NWin1, NWin2, iModel, riMode, ctFlag;

    int    i, error;

    error  = 0;

    NSig1  = fInfo->NSig1;
    NSig2  = fInfo->NSig2;
    NSp1   = fInfo->NSp1;
    NSp2   = fInfo->NSp2;
    Nb1    = fInfo->Nb1;
    Nb2    = fInfo->Nb2;
    Nb     = fInfo->Nb;
    NWin1  = fInfo->NWin1;
    NWin2  = fInfo->NWin2;
    delta1 = fInfo->delta1;
    delta2 = fInfo->delta2;
    wmin1  = fInfo->wmin1;
    wmax1  = fInfo->wmax1;
    wmin2  = fInfo->wmin2;
    wmax2  = fInfo->wmax2;
    dW1    = fInfo->dW1;
    dW2    = fInfo->dW2;
    gamma1 = fInfo->gamma1;
    gamma2 = fInfo->gamma2;
    ctFlag = fInfo->ctFlag;
    q2     = fInfo->q2;
    iModel = fInfo->iModel;
    riMode = fInfo->riMode;
    s      = fInfo->scale;

    (void) fdm2dmw_( fInfo->td, 
		     &NSig1, 
		     &NSig2,
		     &NSp1,
		     &NSp2,
		     &delta1, 
		     &delta2, 
		     &Nb1, 
		     &Nb2, 
		     &Nb, 
		     &NWin1,
		     &NWin2,
		     &wmin1, 
		     &wmax1, 
		     &wmin2, 
		     &wmax2,
		     &dW1,
		     &dW2,
		     &gamma1,
		     &gamma2,
		     &q2,
		     &ctFlag,
		     &iModel,
		     &riMode,
		     &s,
		     fInfo->riSpec,
		     fInfo->wk, 
		     fInfo->dk,
		     fInfo->U, 
		     fInfo->zz, 
		     fInfo->bk, 
		     fInfo->coefW, 
		     fInfo->f, 
		     fInfo->g, 
		     fInfo->gz, 
		     fInfo->gsave, 
		     fInfo->gzct, 
		     fInfo->z2_Mt, 
		     fInfo->cm1, 
		     fInfo->cm2, 
		     fInfo->beta,
		     fInfo->wSpec,
		     fInfo->v,
		     &error );

    return( error );
}

/***/
/* fdm1DCNull: null-initialize 1D FDM structure.
/***/

int fdm1DCNull( fInfo )

   struct FDM1DInfo *fInfo;
{
    fInfo->NSig     = 0; 
    fInfo->NSp      = 0;     
    fInfo->delta    = 0.0;   
    fInfo->gamma    = -1.0;   
    fInfo->Nb0      = 0;    
    fInfo->Nbc      = 0;    
    fInfo->Nb       = 0; 
    fInfo->NWin     = 0;
    fInfo->wmin     = -LARGENUM;  
    fInfo->wmax     =  LARGENUM;
    fInfo->dW       = 0.0;
    fInfo->rho      = 1.1;  
    fInfo->msFlag   = 0; 
    fInfo->mwFlag   = 1;

    fInfo->td       = (FComplex *)NULL;  
    fInfo->wk       = (FComplex *)NULL;    
    fInfo->dk       = (FComplex *)NULL;    
    fInfo->U        = (FComplex *)NULL;     
    fInfo->f        = (FComplex *)NULL;     
    fInfo->g        = (FComplex *)NULL;     
    fInfo->diag     = (FComplex *)NULL;  
    fInfo->z        = (FComplex *)NULL;     
    fInfo->beta     = (FDouble *)NULL; 
    fInfo->coefW    = (FComplex *)NULL; 
    fInfo->zz       = (FComplex *)NULL;    
    fInfo->U0       = (FComplex *)NULL;    
    fInfo->riSpec   = (FComplex *)NULL;    
    fInfo->wSpec    = (FComplex *)NULL;    

    return( 0 );
}

/***/
/* fdm2DCNull: null-initialize 2D FDM structure.
/***/

int fdm2DCNull( fInfo )

   struct FDM2DInfo *fInfo;
{
    fInfo->NSig1    = 0;   
    fInfo->NSig2    = 0;   
    fInfo->NSp1     = 0;     
    fInfo->NSp2     = 0;     
    fInfo->Nb1      = 0;     
    fInfo->Nb2      = 0;     
    fInfo->Nb       = 0;     
    fInfo->NWin1    = 0;
    fInfo->NWin2    = 0;
    fInfo->sx1      = 1.0;
    fInfo->sxn      = LARGENUM;
    fInfo->sy1      = 1.0;
    fInfo->syn      = LARGENUM;
    fInfo->wmin1    = -LARGENUM;  
    fInfo->wmax1    =  LARGENUM;  
    fInfo->wmin2    = -LARGENUM;  
    fInfo->wmax2    =  LARGENUM;
    fInfo->dW1      = 0.0;
    fInfo->dW2      = 0.0;
    fInfo->delta1   = 0.0;   
    fInfo->delta2   = 0.0;   
    fInfo->gamma1   = -1.0;   
    fInfo->gamma2   = -1.0;   
    fInfo->rho1     = 1.1;  
    fInfo->rho2     = 1.1;  
    fInfo->q2       = 1.0e-4;
    fInfo->mwxFlag  = 1;
    fInfo->mwyFlag  = 0;
    fInfo->singFlag = 0;
    fInfo->ctFlag   = 0;

    fInfo->iModel   = 1;
    fInfo->riMode   = 1;
    fInfo->scale    = 1.0;

    fInfo->td       = (FComplex *)NULL;  
    fInfo->wk       = (FComplex *)NULL;    
    fInfo->dk       = (FComplex *)NULL;    
    fInfo->U        = (FComplex *)NULL;     
    fInfo->zz       = (FComplex *)NULL;    
    fInfo->bk       = (FComplex *)NULL;     
    fInfo->coefW    = (FComplex *)NULL; 
    fInfo->f        = (FComplex *)NULL;     
    fInfo->g        = (FComplex *)NULL;     
    fInfo->gz       = (FComplex *)NULL;     
    fInfo->gsave    = (FComplex *)NULL;     
    fInfo->gzct     = (FComplex *)NULL;     
    fInfo->z2_Mt    = (FComplex *)NULL;  
    fInfo->cm1      = (FComplex *)NULL;    
    fInfo->cm2      = (FComplex *)NULL;    
    fInfo->beta     = (FDouble *)NULL; 
    fInfo->v        = (FComplex *)NULL;    
    fInfo->wSpec    = (FComplex *)NULL; 
    fInfo->riSpec   = (FComplex *)NULL;

    return( 0 );
}

/***/
/* fdm1DCAlloc: Allocate workspace for 1D FDM.
/***/

int fdm1DCAlloc( fInfo )

   struct FDM1DInfo *fInfo;
{
    int Nb, NSig, NSp, error;
 
    error = 0; 
    Nb    = fInfo->Nb + 1;
    NSig  = fInfo->NSig + 1;
    NSp   = fInfo->NSp + 1;
 
    MTEST( "FDM", fInfo->td     = FCOMPLEX_ALLOC( NSig+1 ));
    MTEST( "FDM", fInfo->wk     = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->dk     = FCOMPLEX_ALLOC( Nb )); 
    MTEST( "FDM", fInfo->U      = FCOMPLEX_ALLOC( Nb*Nb*2 ));
    MTEST( "FDM", fInfo->f      = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->g      = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->diag   = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->z      = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->coefW  = FCOMPLEX_ALLOC( NSig ));
    MTEST( "FDM", fInfo->zz     = FCOMPLEX_ALLOC( Nb*Nb ));
    MTEST( "FDM", fInfo->U0     = FCOMPLEX_ALLOC( Nb*Nb ));
    MTEST( "FDM", fInfo->beta   = FDOUBLE_ALLOC( Nb ));
    MTEST( "FDM", fInfo->riSpec = FCOMPLEX_ALLOC( NSp ));
    MTEST( "FDM", fInfo->wSpec  = FCOMPLEX_ALLOC( NSp ));

shutdown:
 
    return( error );
}

/***/
/* fdm2DCAlloc: Allocate workspace for 2D FDM.
/***/

int fdm2DCAlloc( fInfo )

   struct FDM2DInfo *fInfo;
{
    int Nb, Nb1, Nb2, NSig1, NSig2, NSp1, NSp2, NSigMax, NSpMax, NBMax, error;
 
    error = 0; 

    Nb      = fInfo->Nb + 1;
    Nb1     = fInfo->Nb1 + 1;
    Nb2     = fInfo->Nb2 + 1;
    NSig1   = fInfo->NSig1 + 1;
    NSig2   = fInfo->NSig2 + 1;
    NSp1    = fInfo->NSp1 + 1;
    NSp2    = fInfo->NSp2 + 1;
    NSigMax = NSig1 > NSig2 ? NSig1 : NSig2;
    NSpMax  = NSp1  > NSp2  ? NSp1  : NSp2;
    NBMax   = Nb1   > Nb2   ? Nb1   : Nb2;
 
    MTEST( "FDM", fInfo->td     = FCOMPLEX_ALLOC( (NSig1+1)*(NSig2+1) ));
    MTEST( "FDM", fInfo->wk     = FCOMPLEX_ALLOC( Nb*2 ));
    MTEST( "FDM", fInfo->dk     = FCOMPLEX_ALLOC( Nb*Nb )); 
    MTEST( "FDM", fInfo->U      = FCOMPLEX_ALLOC( Nb*Nb*3 ));
    MTEST( "FDM", fInfo->zz     = FCOMPLEX_ALLOC( Nb*Nb*2 ));
    MTEST( "FDM", fInfo->bk     = FCOMPLEX_ALLOC( Nb*2 ));
    MTEST( "FDM", fInfo->coefW  = FCOMPLEX_ALLOC( NSigMax+1 ));
    MTEST( "FDM", fInfo->f      = FCOMPLEX_ALLOC( (NSigMax+1)*NBMax*3*2 ));
    MTEST( "FDM", fInfo->g      = FCOMPLEX_ALLOC( Nb*3*3 ));
    MTEST( "FDM", fInfo->gz     = FCOMPLEX_ALLOC( Nb*3*3*3 ));
    MTEST( "FDM", fInfo->gsave  = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->gzct   = FCOMPLEX_ALLOC( Nb ));
    MTEST( "FDM", fInfo->z2_Mt  = FCOMPLEX_ALLOC( NBMax ));
    MTEST( "FDM", fInfo->cm1    = FCOMPLEX_ALLOC( Nb*Nb ));
    MTEST( "FDM", fInfo->cm2    = FCOMPLEX_ALLOC( Nb*Nb ));
    MTEST( "FDM", fInfo->beta   = FDOUBLE_ALLOC( Nb ));

    MTEST( "FDM", fInfo->wSpec  = FCOMPLEX_ALLOC( NSp1*NSp2 ));
    MTEST( "FDM", fInfo->riSpec = FCOMPLEX_ALLOC( NSp1*NSp2 ));
    MTEST( "FDM", fInfo->v      = FCOMPLEX_ALLOC( Nb*NSpMax ));

shutdown:
 
    return( error );
}

/***/
/* fdmNDCFree: De-Allocate workspace for FDM.
/***/

int fdmNDCFree( fInfo )

   struct FDMInfo *fInfo;
{            
    if (fInfo->dimCount == 1)
       (void) fdm1DCFree( &fInfo->fInfo1D );
    else if (fInfo->dimCount == 2)
       (void) fdm2DCFree( &fInfo->fInfo2D );

    (void) fdmNDCNull( fInfo );

    return( 0 );
}

/***/
/* fdm1DCFree: De-Allocate workspace for 1D FDM.
/***/

int fdm1DCFree( fInfo )

   struct FDM1DInfo *fInfo;
{
    int Nb, NSig, NSp;

    Nb   = fInfo->Nb + 1;
    NSig = fInfo->NSig + 1;
    NSp  = fInfo->NSp + 1;

    FCOMPLEX_DEALLOC( fInfo->wk,     Nb );
    FCOMPLEX_DEALLOC( fInfo->dk,     Nb );
    FCOMPLEX_DEALLOC( fInfo->td,     NSig+1 );
    FCOMPLEX_DEALLOC( fInfo->U,      Nb*Nb*2 );
    FCOMPLEX_DEALLOC( fInfo->f,      Nb );
    FCOMPLEX_DEALLOC( fInfo->g,      Nb );
    FCOMPLEX_DEALLOC( fInfo->diag,   Nb );
    FCOMPLEX_DEALLOC( fInfo->z,      Nb );
    FCOMPLEX_DEALLOC( fInfo->coefW,  NSig );
    FCOMPLEX_DEALLOC( fInfo->zz,     Nb*Nb );
    FCOMPLEX_DEALLOC( fInfo->U0,     Nb*Nb );
    FDOUBLE_DEALLOC( fInfo->beta,    Nb );
    FCOMPLEX_DEALLOC( fInfo->riSpec, NSp );
    FCOMPLEX_DEALLOC( fInfo->wSpec,  NSp );

    (void) fdm1DCNull( fInfo );

    return( 0 );
}

/***/
/* fdm2DCFree: De-Allocate workspace for 2D FDM.
/***/

int fdm2DCFree( fInfo )

   struct FDM2DInfo *fInfo;
{
    int Nb, Nb1, Nb2, NSig1, NSig2, NSp1, NSp2, NSigMax, NSpMax, NBMax, error;
 
    error = 0; 

    Nb      = fInfo->Nb + 1;
    Nb1     = fInfo->Nb1 + 1;
    Nb2     = fInfo->Nb2 + 1;
    NSig1   = fInfo->NSig1 + 1;
    NSig2   = fInfo->NSig2 + 1;
    NSp1    = fInfo->NSp1 + 1;
    NSp2    = fInfo->NSp2 + 1;
    NSigMax = NSig1 > NSig2 ? NSig1 : NSig2;
    NSpMax  = NSp1  > NSp2  ? NSp1  : NSp2;
    NBMax   = Nb1   > Nb2   ? Nb1   : Nb2;

    FCOMPLEX_DEALLOC( fInfo->td,     (NSig1+1)*(NSig2+1) );
    FCOMPLEX_DEALLOC( fInfo->wk,     Nb*2 );
    FCOMPLEX_DEALLOC( fInfo->dk,     Nb*Nb ); 
    FCOMPLEX_DEALLOC( fInfo->U,      Nb*Nb*3 );
    FCOMPLEX_DEALLOC( fInfo->zz,     Nb*Nb*2 );
    FCOMPLEX_DEALLOC( fInfo->bk,     Nb*2 );
    FCOMPLEX_DEALLOC( fInfo->coefW,  NSigMax+1 );
    FCOMPLEX_DEALLOC( fInfo->f,      (NSigMax+1)*NBMax*3*2 );
    FCOMPLEX_DEALLOC( fInfo->g,      Nb*3*3 );
    FCOMPLEX_DEALLOC( fInfo->gz,     Nb*3*3*3 );
    FCOMPLEX_DEALLOC( fInfo->gsave,  Nb );
    FCOMPLEX_DEALLOC( fInfo->gzct,   Nb );
    FCOMPLEX_DEALLOC( fInfo->z2_Mt,  NBMax );
    FCOMPLEX_DEALLOC( fInfo->cm1,    Nb*Nb );
    FCOMPLEX_DEALLOC( fInfo->cm2,    Nb*Nb );
    FDOUBLE_DEALLOC( fInfo->beta,    Nb );

    FCOMPLEX_DEALLOC( fInfo->wSpec,  NSp1*NSp2 );
    FCOMPLEX_DEALLOC( fInfo->riSpec, NSp1*NSp2 );
    FCOMPLEX_DEALLOC( fInfo->v,      Nb*NSpMax );

shutdown:

    return( error );
}

 
int showFDM1DInfo( fInfo )

   struct FDM1DInfo *fInfo;
{
    FPR( stderr, "NSig: %6d\n", fInfo->NSig );
    FPR( stderr, "NSp:  %6d\n", fInfo->NSp );
    FPR( stderr, "Nb0:  %6d\n", fInfo->Nb0 );
    FPR( stderr, "Nbc:  %6d\n", fInfo->Nbc );
    FPR( stderr, "Nb:   %6d\n", fInfo->Nb );
    FPR( stderr, "NWin: %6d\n", fInfo->NWin );

    FPR( stderr, "sx1:   % e\n", fInfo->sx1 );      
    FPR( stderr, "sxn:   % e\n", fInfo->sxn );      

    FPR( stderr, "wmin:  % le\n", fInfo->wmin );     
    FPR( stderr, "wmax:  % le\n", fInfo->wmax );     
    FPR( stderr, "dW:    % le\n", fInfo->dW );       
    FPR( stderr, "gamma: % le\n", fInfo->gamma );    
    FPR( stderr, "delta: % le\n", fInfo->delta );    
    FPR( stderr, "rho:   % le\n", fInfo->rho );      

    FPR( stderr, "mwFlag: %s\n", fInfo->mwFlag ? "True" : "False" );
    FPR( stderr, "msFlag: %s\n", fInfo->msFlag ? "True" : "False" );

    return( 0 );
}

int showFDM2DInfo( fInfo )

   struct FDM2DInfo *fInfo;
{
    FPR( stderr, "NSig1: %6d\n", fInfo->NSig1 );
    FPR( stderr, "NSig2: %6d\n", fInfo->NSig2 );
    FPR( stderr, "NSp1:  %6d\n", fInfo->NSp1 );
    FPR( stderr, "NSp2:  %6d\n", fInfo->NSp2 );
    FPR( stderr, "Nb1:  %6d\n", fInfo->Nb1 );
    FPR( stderr, "Nb2:  %6d\n", fInfo->Nb2 );
    FPR( stderr, "Nb:   %6d\n", fInfo->Nb );
    FPR( stderr, "NWin1: %6d\n", fInfo->NWin1 );
    FPR( stderr, "NWin2: %6d\n", fInfo->NWin2 );

    FPR( stderr, "sx1:   % e\n", fInfo->sx1 );      
    FPR( stderr, "sxn:   % e\n", fInfo->sxn );      
    FPR( stderr, "sy1:   % e\n", fInfo->sy1 );      
    FPR( stderr, "syn:   % e\n", fInfo->syn );      

    FPR( stderr, "wmin1:  % le\n", fInfo->wmin1 );
    FPR( stderr, "wmin2:  % le\n", fInfo->wmin2 );
    FPR( stderr, "wmax1:  % le\n", fInfo->wmax1 );
    FPR( stderr, "wmax2:  % le\n", fInfo->wmax2 );
    FPR( stderr, "dW1:    % le\n", fInfo->dW1 );
    FPR( stderr, "dW2:    % le\n", fInfo->dW2 );
    FPR( stderr, "delta1: % le\n", fInfo->delta1 );
    FPR( stderr, "delta2: % le\n", fInfo->delta2 );
    FPR( stderr, "gamma1: % le\n", fInfo->gamma1 );
    FPR( stderr, "gamma2: % le\n", fInfo->gamma2 );
    FPR( stderr, "rho1:   % le\n", fInfo->rho1 );
    FPR( stderr, "rho2:   % le\n", fInfo->rho2 );
    FPR( stderr, "q2:     % le\n", fInfo->q2 );

    FPR( stderr, "mwxFlag:  %s\n", fInfo->mwxFlag  ? "True" : "False" );
    FPR( stderr, "mwyFlag:  %s\n", fInfo->mwyFlag  ? "True" : "False" );
    FPR( stderr, "singFlag: %s\n", fInfo->singFlag ? "True" : "False" );
    FPR( stderr, "ctFlag:   %s\n", fInfo->ctFlag   ? "True" : "False" );

    FPR( stderr, "iModel: %d\n", fInfo->iModel );
    FPR( stderr, "riMode: %d\n", fInfo->riMode );

    FPR( stderr, "scale:  % le\n", fInfo->scale );

    return( 0 );
}

int showFDMInfo( fInfo )

   struct FDMInfo *fInfo;
{
    FPR( stderr, "dimCount: %d\n", fInfo->dimCount );

    if (fInfo->dimCount == 1)
       (void) showFDM1DInfo( &fInfo->fInfo1D );
    else if (fInfo->dimCount == 2)
       (void) showFDM2DInfo( &fInfo->fInfo2D );

    return( 0 );
}