/// <summary>
		/// Calculate the square root of a complex number
		/// </summary>
		/// <param name="c"></param>
		/// <returns></returns>
		static public ComplexF	Sqrt( ComplexF c ) {
			double	x	= c.Re;
			double	y	= c.Im;

			double	modulus	= Math.Sqrt( x*x + y*y );
			int		sign	= ( y < 0 ) ? -1 : 1;

			c.Re		= (float)( _halfOfRoot2 * Math.Sqrt( modulus + x ) );
			c.Im	= (float)( _halfOfRoot2 * sign * Math.Sqrt( modulus - x ) );

			return	c;
		}

        //---------------------------------------------------------------------------------------------------
        /// <summary>
        /// Calculate the power of a complex number
        /// </summary>
        /// <param name="c"></param>
        /// <param name="exponent"></param>
        /// <returns></returns>
        public static ComplexF Pow( ComplexF c, double exponent )
        {
            double	x	= c.Re;
            double	y	= c.Im;

            double	modulus		= Math.Pow( x*x + y*y, exponent * 0.5 );
            double	argument	= Math.Atan2( y, x ) * exponent;

            c.Re		= (float)( modulus * System.Math.Cos( argument ) );
            c.Im = (float)( modulus * System.Math.Sin( argument ) );

            return	c;
        }

		static private ComplexF		SumRecursion( ComplexF[] data, int start, int end ) {
			Debug.Assert( 0 <= start, "start = " + start );
			Debug.Assert( start < end, "start = " + start + " and end = " + end );
			Debug.Assert( end <= data.Length, "end = " + end + " and data.Length = " + data.Length );
			if( ( end - start ) <= 1000 ) {
				ComplexF sum = ComplexF.Zero;
				for( int i = start; i < end; i ++ ) {
					sum += data[ i ];
				
				}
				return	sum;
			}
			else {
				int middle = ( start + end ) >> 1;
				return	SumRecursion( data, start, middle ) + SumRecursion( data, middle, end );
			}
		}

    void InitWaveGenerator()
    {
        // Wind restricted to one direction, reduces calculations
        Vector2 wind = new Vector2 (windx, 0.0f);

        // Initialize wave generator
        for (int y=0; y<height; y++) {
            for (int x=0; x<width; x++) {
                float yc = y < height / 2f ? y : -height + y;
                float xc = x < width / 2f ? x : -width + x;
                Vector2 vec_k = new Vector2 (2.0f * Mathf.PI * xc / size.x, 2.0f * Mathf.PI * yc / size.z);
                h0 [width * y + x] = new ComplexF (GaussianRnd (), GaussianRnd ()) * 0.707f * (float)System.Math.Sqrt (P_spectrum (vec_k, wind));
            }
        }

        for (int y=0; y<n_height; y++) {
            for (int x=0; x<n_width; x++) {
                float yc = y < n_height / 2f ? y : -n_height + y;
                float xc = x < n_width / 2f ? x : -n_width + x;
                Vector2 vec_k = new Vector2 (2.0f * Mathf.PI * xc / (size.x / normal_scale), 2.0f * Mathf.PI * yc / (size.z / normal_scale));
                n0 [n_width * y + x] = new ComplexF (GaussianRnd (), GaussianRnd ()) * 0.707f * (float)System.Math.Sqrt (P_spectrum (vec_k, wind));
            }
        }
    }

		/// <summary>
		/// Invert each element in the array
		/// </summary>
		/// <param name="array"></param>
		static public void Invert( ComplexF[] array ) {
			for( int i = 0; i < array.Length; i ++ ) {
				array[i] = ((ComplexF) 1 ) / array[i];
			}
		}

		/// <summary>
		/// Copy an array
		/// </summary>
		/// <param name="dest"></param>
		/// <param name="source"></param>
		static public void Copy( ComplexF[] dest, ComplexF[] source ) {
			Debug.Assert( dest != null );
			Debug.Assert( source != null );
			Debug.Assert( dest.Length == source.Length );
			for( int i = 0; i < dest.Length; i ++ ) {
				dest[i] = source[i];
			}
		}

		/// <summary>
		/// Divide each element in target array with corresponding element in rhs array
		/// </summary>
		/// <param name="target"></param>
		/// <param name="rhs"></param>
		static public void Divide( ComplexF[] target, ComplexF[] rhs ) {
			ComplexArray.Divide( target, rhs, target );
		}

		/// <summary>
		/// Multiply each element in target array with corresponding element in rhs array
		/// </summary>
		/// <param name="target"></param>
		/// <param name="rhs"></param>
		static public void Multiply( ComplexF[] target, ComplexF[] rhs ) {
			ComplexArray.Multiply( target, rhs, target );
		}

		/// <summary>
		/// Multiply each element in the array by a specific value
		/// </summary>
		/// <param name="array"></param>
		/// <param name="scale"></param>
		static public void Scale( ComplexF[] array, ComplexF scale ) {
			Debug.Assert( array != null );

			int length = array.Length;
			for( int i = 0; i < length; i ++ ) {
				array[i] *= scale;
			}
		}

		static private void		LockBufferCF( int length, ref ComplexF[] buffer ) {
			
			if( length != _bufferCF.Length ) {
				_bufferCF	= new ComplexF[ length ];
			}
			buffer =	_bufferCF;
		}

		//-------------------------------------------------------------------------------------



		static private void ReorderArray( ComplexF[] data ) {

			int length = data.Length;

			int[] reversedBits = GetReversedBits( Log2( length ) );
			for( int i = 0; i < length; i ++ ) {
				int swap = reversedBits[ i ];
				if( swap > i ) {
					ComplexF temp = data[ i ];
					data[ i ] = data[ swap ];
					data[ swap ] = temp;
				}
			}
		}

		/// <summary>
		/// Swap two complex numbers
		/// </summary>
		/// <param name="a"></param>
		/// <param name="b"></param>
		static public void Swap( ref ComplexF a, ref ComplexF b ) {
			ComplexF temp = a;
			a = b;
			b = temp;
		}

		/// <summary>
		/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
		/// </summary>
		/// <param name="data"></param>
		/// <param name="direction"></param>
		public static void	FFT( ComplexF[] data, FourierDirection direction ) {
			if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			Fourier.FFT( data, data.Length, direction );
		}

		/// <summary>
		/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
		/// </summary>
		/// <param name="data"></param>
		/// <param name="length"></param>
		/// <param name="direction"></param>
		public static void	FFT_Quick( ComplexF[] data, int length, FourierDirection direction ) {
			/*if( data == null ) {
				throw new ArgumentNullException( "data" );
			}
			if( data.Length < length ) {
				throw new ArgumentOutOfRangeException( "length", length, "must be at least as large as 'data.Length' parameter" );
			}
			if( Fourier.IsPowerOf2( length ) == false ) {
				throw new ArgumentOutOfRangeException( "length", length, "must be a power of 2" );
			}

			Fourier.SyncLookupTableLength( length );*/

			int ln = Fourier.Log2( length );
			
			// reorder array
			Fourier.ReorderArray( data );
			
			// successive doubling
			int N = 1;
			int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;
			
			for( int level = 1; level <= ln; level ++ ) {
				int M = N;
				N <<= 1;

				float[] uRLookup = _uRLookupF[ level, signIndex ];
				float[] uILookup = _uILookupF[ level, signIndex ];

				for( int j = 0; j < M; j ++ ) {
					float uR = uRLookup[j];
					float uI = uILookup[j];

					for( int even = j; even < length; even += N ) {
						int odd	 = even + M;
						
						float	r = data[ odd ].Re;
						float	i = data[ odd ].Im;

						float	odduR = r * uR - i * uI;
						float	odduI = r * uI + i * uR;

						r = data[ even ].Re;
						i = data[ even ].Im;
						
						data[ even ].Re	= r + odduR;
						data[ even ].Im	= i + odduI;
						
						data[ odd ].Re	= r - odduR;
						data[ odd ].Im	= i - odduI;
					}
				}
			}

		}

		static private void Swap( ref ComplexF a, ref ComplexF b ) {
			ComplexF temp = a;
			a = b;
			b = temp;
		}

		/// <summary>
		///  Determine whether the elements in the two arrays are the same
		/// </summary>
		/// <param name="array1"></param>
		/// <param name="array2"></param>
		/// <param name="tolerance"></param>
		/// <returns></returns>
		static public bool		IsEqual( ComplexF[] array1, ComplexF[] array2, float tolerance ) {
			if ( array1.Length != array2.Length ) {
				return false;
			}
			for( int i = 0; i < array1.Length; i ++ ) {
				if( ComplexF.IsEqual( array1[i], array2[i], tolerance ) == false ) {
					return false;
				}
			}
			return true;
		}

		/// <summary>
		/// Add a specific value to each element in the array
		/// </summary>
		/// <param name="array"></param>
		/// <param name="offset"></param>
		static public void Offset( ComplexF[] array, ComplexF offset ) {
			int length = array.Length;
			for( int i = 0; i < length; i ++ ) {
				array[i] += offset;
			}
		}

		static private void		UnlockBufferCF( ref ComplexF[] buffer ) {
			buffer = null;
		}

		private static void	LinearFFT_Quick( ComplexF[] data, int start, int inc, int length, FourierDirection direction ) {
			
			// copy to buffer
			ComplexF[]	buffer = null;
			LockBufferCF( length, ref buffer );
			int j = start;
			for( int i = 0; i < length; i ++ ) {
				buffer[ i ] = data[ j ];
				j += inc;
			}

			FFT( buffer, length, direction );

			// copy from buffer
			j = start;
			for( int i = 0; i < length; i ++ ) {
				data[ j ] = buffer[ i ];
				j += inc;
			}
			UnlockBufferCF( ref buffer );
		}

		//======================================================================================
		 

		/// <summary>
		/// Compute a 1D fast Fourier transform of a dataset of complex numbers.
		/// </summary>
		/// <param name="data"></param>
		/// <param name="length"></param>
		/// <param name="direction"></param>
		public static void	FFT( ComplexF[] data, int length, FourierDirection direction ) {

			SyncLookupTableLength( length );

			int ln = Log2( length );
			
			// reorder array
			ReorderArray( data );
			
			// successive doubling
			int N = 1;
			int signIndex = ( direction == FourierDirection.Forward ) ? 0 : 1;
			
			for( int level = 1; level <= ln; level ++ ) {
				int M = N;
				N <<= 1;

				float[] uRLookup = _uRLookupF[ level, signIndex ];
				float[] uILookup = _uILookupF[ level, signIndex ];

				for( int j = 0; j < M; j ++ ) {
					float uR = uRLookup[j];
					float uI = uILookup[j];

					for( int even = j; even < length; even += N ) {
						int odd	 = even + M;
						
						float	r = data[ odd ].Re;
						float	i = data[ odd ].Im;

						float	odduR = r * uR - i * uI;
						float	odduI = r * uI + i * uR;

						r = data[ even ].Re;
						i = data[ even ].Im;
						
						data[ even ].Re	= r + odduR;
						data[ even ].Im	= i + odduI;
						
						data[ odd ].Re	= r - odduR;
						data[ odd ].Im	= i - odduI;
					}
				}
			}

		}