Index: trunk/psLib/src/fft/psVectorFFT.c
===================================================================
--- trunk/psLib/src/fft/psVectorFFT.c	(revision 1625)
+++ trunk/psLib/src/fft/psVectorFFT.c	(revision 1864)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-08-25 21:10:08 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 18:31:49 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,4 +23,6 @@
 #include "psImageExtraction.h"
 
+#include "psDataManipErrors.h"
+
 #define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
 
@@ -42,5 +44,10 @@
 
     if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {
-        psError(__func__, "Input image must be a 32-bit float or complex image (type=%d)", type);
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_IMAGE_TYPE_UNSUPPORTED,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -48,5 +55,10 @@
 
     if ((type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE)) {
-        psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_REVERSE_NOT_COMPLEX,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -55,5 +67,10 @@
 
     if ((type != PS_TYPE_F32) && (direction == PS_FFT_FORWARD)) {
-        psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_FORWARD_NOT_REAL,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -94,5 +111,7 @@
     /* check if a plan exists now */
     if (plan == NULL) {
-        psError(__func__, "Failed to create FFTW plan.");
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_FFTW_PLAN_NULL);
         psFree(out);
         return NULL;
@@ -142,6 +161,22 @@
             outVec[i] = crealf(inVec[i]);
         }
-    } else {
-        psError(__func__, "Can not extract real component from given vector type (%d).", type);
+    } else if (type == PS_TYPE_C64) {
+        psF64* outVec;
+        psC64* inVec = in->data.C64;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_F64);
+        out->n = numElements;
+        outVec = out->data.F64;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = creal(inVec[i]);
+        }
+    } else {
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorReal",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_TYPE_UNSUPPORTED,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -186,6 +221,22 @@
             outVec[i] = cimagf(inVec[i]);
         }
-    } else {
-        psError(__func__, "Can not extract imaginary component from given vector type (%d).", type);
+    } else if (type == PS_TYPE_C64) {
+        psF64* outVec;
+        psC64* inVec = in->data.C64;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_F64);
+        out->n = numElements;
+        outVec = out->data.F64;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = cimag(inVec[i]);
+        }
+    } else {
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorImaginary",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_TYPE_UNSUPPORTED,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -213,11 +264,12 @@
 
     if (imag->type.type != type) {
-        psError(__func__, "The inputs to psVectorComplex must be the same type.");
-        psFree(out);
-        return NULL;
-    }
-
-    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
-        psError(__func__, "The inputs to psVectorComplex can not be complex.");
+        char* typeStrReal;
+        char* typeStrImag;
+        PS_TYPE_NAME(typeStrReal,type);
+        PS_TYPE_NAME(typeStrImag,imag->type.type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorComplex",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_REAL_IMAG_TYPE_MISMATCH,
+                   typeStrReal,typeStrImag);
         psFree(out);
         return NULL;
@@ -236,6 +288,23 @@
             outVec[i] = realVec[i] + I * imagVec[i];
         }
-    } else {
-        psError(__func__, "Can not merge real and imaginary portions for given vector type (%d).", type);
+    } else if (type == PS_TYPE_F64) {
+        psC64* outVec;
+        psF64* realVec = real->data.F64;
+        psF64* imagVec = imag->data.F64;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_C64);
+        out->n = numElements;
+        outVec = out->data.C64;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = realVec[i] + I * imagVec[i];
+        }
+    } else {
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorComplex",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_NONREAL_NOTSUPPORTED,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -281,6 +350,22 @@
             outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]);
         }
-    } else {
-        psError(__func__, "Can not compute complex conjugate for given vector type (%d).", type);
+    } else if (type == PS_TYPE_C64) {
+        psC64* outVec;
+        psC64* inVec = in->data.C64;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_C64);
+        out->n = numElements;
+        outVec = out->data.C64;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]);
+        }
+    } else {
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorConjugate",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_NONCOMPLEX_NOTSUPPORTED,
+                   typeStr);
         psFree(out);
         return NULL;
@@ -309,11 +394,4 @@
     outNumElements = inHalfNumElements + 1;
 
-    /* if not a complex number, this is not implemented */
-    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
-        psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
-        psFree(out);
-        return NULL;
-    }
-
     if (type == PS_TYPE_C32) {
         psF32* outVec;
@@ -340,11 +418,40 @@
         inAbs1 = cabsf(inVec[inHalfNumElements]);
         outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
-    } else {
-        psError(__func__, "Can not power spectrum for given vector type (%d).", type);
-        psFree(out);
-        return NULL;
-    }
-
-    return out;
-
-}
+    } else if (type == PS_TYPE_C64) {
+        psF64* outVec;
+        psC64* inVec = in->data.C64;
+        psF64 inAbs1;
+        psF64 inAbs2;
+
+        out = psVectorRecycle(out, outNumElements, PS_TYPE_F64);
+        out->n = outNumElements;
+        outVec = out->data.F64;
+
+        // from ADD: P_0 = |C_0|^2/N^2
+        inAbs1 = cabs(inVec[0]);
+        outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
+
+        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
+        for (unsigned int i = 1; i < inHalfNumElements; i++) {
+            inAbs1 = cabs(inVec[i]);
+            inAbs2 = cabs(inVec[inNumElements - i]);
+            outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
+        }
+
+        // from ADD: P_N/2 = |C_N/2|^2/N^2
+        inAbs1 = cabs(inVec[inHalfNumElements]);
+        outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
+    } else {
+        char* typeStr;
+        PS_TYPE_NAME(typeStr,type);
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorPowerSpectrum",
+                   PS_ERR_BAD_PARAMETER_TYPE, true,
+                   PS_ERRORTEXT_psVectorFFT_NONCOMPLEX_NOTSUPPORTED,
+                   typeStr);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}