/***/
/* yz.M: extract the first YZ plane of a 3D or 4D hypercomplex FID.
/*
/* Requires command-line definition of strings "inName" and "outName",
/* and variable "fidLoc" e.g.:
/*
/*   M yz.M -str inName fid/test%03.0f.fid outName yz.fid -var fidLoc 1 -quit
/***/

/***/
/* Parameters taken from "fdatap.h"; these will be pre-defined
/* in future implementations of "M":
/***/

    FDATASIZE = 512;

    float fdata[FDATASIZE];
    char  thisName[256];

    FDDIMCOUNT   =   9;
    FDDIMORDER1  =  24;
    FDDIMORDER2  =  25;
    FDDIMORDER3  =  26;
    FDDIMORDER4  =  27;
    FDFILECOUNT  = 442;
    NULL_DIM     = 0;

/***/
/* Open the first 4D plane.
/* Read the header from the first 4D plane.
/* Extract the size parameters.
/***/

    (void) sprintf( thisName, inName, 1, 1 );

    inFD  = dOpen( thisName, "r" );
    outFD = dOpen( outName,  "w" );

    (void) dReadHdr( inFD, fdata, 0 );
    (void) dClose( inFD );

    dimCount = getParm( fdata, FDDIMCOUNT, NULL_DIM );

    xSize = getParm( fdata, NDSIZE, CUR_XDIM );
    ySize = getParm( fdata, NDSIZE, CUR_YDIM );
    zSize = getParm( fdata, NDSIZE, CUR_ZDIM );
    aSize = getParm( fdata, NDSIZE, CUR_ADIM );

    if (dimCount < 3)
       {
        (void) printf( "Error: the dimension count is too small.\n" );
        exit( 1 );
       };

    if (dimCount < 4)
       {
        aSize = 1;
       };

/***/
/* Adjust header of input to reflect an output with Y and Z axis exchanged:
/***/

    (void) setParm( fdata, NDSIZE, ySize/2, CUR_XDIM );
    (void) setParm( fdata, NDSIZE, zSize,   CUR_YDIM );
    (void) setParm( fdata, NDSIZE, xSize*2, CUR_ZDIM );

    xID = getParm( fdata, FDDIMORDER1, NULL_DIM );
    yID = getParm( fdata, FDDIMORDER2, NULL_DIM );
    zID = getParm( fdata, FDDIMORDER3, NULL_DIM );

    (void) setParm( fdata, FDDIMORDER1, yID, NULL_DIM );
    (void) setParm( fdata, FDDIMORDER2, zID, NULL_DIM );
    (void) setParm( fdata, FDDIMORDER3, xID, NULL_DIM );

    (void) setParm( fdata, FDFILECOUNT,  2*xSize*aSize, NULL_DIM );

    (void) dWrite( outFD, fdata, wordLen*FDATASIZE );

/***/
/* Read a complex 1D point from every vector in every XY plane in cube A=1:
/***/

    float buff[ySize*zSize], pR[1], pI[1];

    ia = 1;

    for( iz = 1; iz <= zSize; iz++ )
       {
        if (dimCount == 3)
           {
            (void) sprintf( thisName, inName, iz );
           }
        else
           {
            (void) sprintf( thisName, inName, ia, iz );
           }

        inFD = dOpen( thisName, "r" );
        iR   = ySize*(iz - 1);
        iI   = iR + ySize/2;

        for( iy = 1; iy <= ySize/2; iy++ )
           {
            offset = wordLen*4*xSize*(iy - 1) + wordLen*(fidLoc - 1);

            (void) dSeek( inFD, wordLen*FDATASIZE + offset );
            (void) dRead( inFD, pR, wordLen );

            offset += wordLen*2*xSize;

            (void) dSeek( inFD, wordLen*FDATASIZE + offset );
            (void) dRead( inFD, pI, wordLen );

            buff[iR++] = pR[0];
            buff[iI++] = pI[0];
           }

        (void) dClose( inFD ); 
       }

    (void) dWrite( outFD, buff, wordLen*ySize*zSize );

    (void) dClose( outFD );
