int main(int argc, char* argv[]){

	const int NP=256;				// Number of points in curve
	const int N1=31, N2=N1*N1;		// Size of graph
	float pixels[N2];				// Accumulation buffer

	Complex<double> phase(1,0), freq;
	

	printf("\nUnit circle\n");
	mem::zero(pixels, N2);
	freq.fromPhase(M_2PI/NP);
	for(int i=0; i<NP; ++i){
		int ix = posToInd(phase[0], N1);
		int iy = posToInd(phase[1], N1);
		phase *= freq;
		pixels[iy*N1 + ix] += 0.125;
	}
	print2D(pixels, N1, N1);


	printf("\nHalf circle\n");
	mem::zero(pixels, N2);
	phase(0.5, 0);
	freq.fromPhase(M_2PI/NP);
	for(int i=0; i<NP; ++i){
		int ix = posToInd(phase[0], N1);
		int iy = posToInd(phase[1], N1);
		phase *= freq;
		pixels[iy*N1 + ix] += 0.125;
	}
	print2D(pixels, N1, N1);	

	return 0;
}

Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Eo() const
{
    const uniformDimensionedVectorField& g =
        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");

    return
        mag(g)*sqr(phase().d())
       *(otherPhase().rho() - phase().rho())
       /fluid().sigma();
}

main()
{
  complex a, b, c, n[2];;
  double real1, imag1, real2, imag2, d, e, f, g, mag, pha;
/*
  a.real = 1.0; a.imag = 1.0;
  b.real = 2.0; b.imag = 2.0;
*/
  printf("\nEnter real part of first number = ");
  scanf("%lf", &real1);
  printf("\nEnter imaginary part of first number = ");
  scanf("%lf", &imag1);
  printf("\nEnter real part of second number = ");
  scanf("%lf", &real2);
  printf("\nEnter imaginary part of second number = ");
  scanf("%lf", &imag2);
  a.real = real1;
  a.imag = imag1;
  b.real = real2;
  b.imag = imag2;
  n[1].real = real1;
  n[1].imag = imag1;
  n[2].real = real2;
  n[2].imag = imag2;

  printf("a = "); print_complex(a); printf("\n");
  printf("b = "); print_complex(b); printf("\n");
  c = add(a,b);
  printf("a + b = "); print_complex(c); printf("\n");
  c = subtract(a,b);
  printf("a - b = "); print_complex(c); printf("\n");
  c = multiply(a,b);
  printf("a * b = "); print_complex(c); printf("\n");
  c = divide(a,b);
  printf("a / b = "); print_complex(c); printf("\n");
  d = realpart(a);
  e = imaginarypart(a);
  printf("real part of a = %lf\n", d);
  printf("imaginary part of a = %lf\n", e);
  f = realpart(b);
  g = imaginarypart(b);
  printf("real part of b = %lf\n", f);
  printf("imaginary part of b = %lf\n", g);
  mag = magnitude(a);
  printf("magnitude of a = %lf\n", mag);
  pha =phase(a);
  printf("phase in radians of a = %lf\n", pha);
  mag = magnitude(b);
  printf("magnitude of b = %lf\n", mag);
  pha = phase(b);
  printf("phase in radians of b = %lf\n", pha);
} 

	double Execute(/*float* data=NULL, int smp*/)
	{
		float x, y, y1, y2, y3, mag, maxL=0.0f, maxR=0.0f, dbin, fc, fcR, fcL, phL, phR;
		/*if (data) memcpy(_inBuf, data, sizeof(float) * chunk_size * 2);*/
		fftwf_execute(_plan);
		int i, maxLpos, maxRpos;
		for (i=0;i<(chunk_size/2+1);++i)
		{
			x = _outBuf[i*2][0];
			y = _outBuf[i*2][1];
			mag = sqrtf(x*x+y*y);
			if (mag > maxL)
			{
				maxL = mag;
				phL = phase(x,y);
				maxLpos = i;
			}
			x = _outBuf[i*2+1][0];
			y = _outBuf[i*2+1][1];
			mag = sqrtf(x*x+y*y);
			if (mag > maxR)
			{
				maxR = mag;
				phR = phase(x,y);
				maxRpos = i;
			}
		}
		x = _outBuf[(maxLpos-1)*2][0];
		y = _outBuf[(maxLpos-1)*2][1];
		y1 = sqrtf(x*x+y*y);
		y2 = maxL;
		x = _outBuf[(maxLpos+1)*2][0];
		y = _outBuf[(maxLpos+1)*2][1];
		y3 = sqrtf(x*x+y*y);
		dbin = ((y3-y1)/(y1+y2+y3));
		fcL = (float(maxLpos)+dbin)*(SAMPLERATE/chunk_size);
		x = _outBuf[(maxRpos-1)*2+1][0];
		y = _outBuf[(maxRpos-1)*2+1][1];
		y1 = sqrtf(x*x+y*y);
		y2 = maxR;
		x = _outBuf[(maxRpos+1)*2+1][0];
		y = _outBuf[(maxRpos+1)*2+1][1];
		y3 = sqrtf(x*x+y*y);
		dbin = ((y3-y1)/(y1+y2+y3));
		fcR = (float(maxRpos)+dbin)*(SAMPLERATE/chunk_size);
		fc = (fcL+fcR)*0.5f;

		//printf("_speed %f fc %f magL %f magR %f phaseL - phaseR %f\n", 1000.0/fc, fc, maxL, maxR, phL- phR);
		if (maxL < 1.0f)
			return 0.0;
		return 2000.0/fc;
	}

void Wf_return::setVals(Array2 <log_value<doublevar> > & v ) {
  is_complex=0;
  int nfunc=v.GetDim(0);
  int nst=v.GetDim(1);
  Resize(nfunc,nst);
  for(int f=0; f< nfunc; f++) {
    amp(f,0)=v(f,0).logval;
    phase(f,0)=v(f,0).sign<0?pi:0.0;
    for(int s=1; s< nst; s++) {
      amp(f,s)=v(f,s).val();
      phase(f,s)=0.0;
    }
  }
}

template<int N> template<class dtype> void LogPhaseProb<N>::set_abcd( const datatypes::ABCD<dtype>& abcd )
{
  if ( !abcd.missing() ) {
    ftype c, s;
    for ( int p = 0; p < q.size(); p++ ) {
      c = cos( phase(p) );
      s = sin( phase(p) );
      q[p] = abcd.a()*c + abcd.b()*s
	+ abcd.c()*(c*c-s*s) + abcd.d()*(2.0*c*s);
    }
  } else {
    for ( int p = 0; p < q.size(); p++ ) q[p] = 0.0;
  }
}

NEMO_PLUGIN_DLL_PUBLIC
void
cpu_update_neurons(
		unsigned start, unsigned end,
		unsigned cycle,
		float* paramBase, size_t paramStride,
		float* stateBase, size_t stateHistoryStride, size_t stateVarStride,
		unsigned fbits,
		unsigned fstim[],
		RNG rng[],
		float /*currentEPSP*/[],
		float /*currentIPSP*/[],
		float /*currentExternal*/[],
		uint64_t recentFiring[],
		unsigned fired[],
		void* rcm_ptr)
{
	const nemo::runtime::RCM& rcm = *static_cast<nemo::runtime::RCM*>(rcm_ptr);

	const float* frequency = paramBase + 0 * paramStride;
 	const float* g_Cmean = paramBase + 1 * paramStride;

	const float* phase0 = phase(stateBase, stateHistoryStride, cycle);// current
	float* phase1 = phase(stateBase, stateHistoryStride, cycle+1);    // next

	std::vector<float> weight;
	std::vector<float> sourcePhase;

	for(unsigned n=start; n < end; n++) {
		float h = 0.05;	
		const float f = frequency[n];
		float targetPhase = phase0[n];

		const unsigned indegree = loadIncoming(rcm, n, cycle, stateBase, stateHistoryStride, weight, sourcePhase);
		float Cmean = g_Cmean[n] > 0 ?  g_Cmean[n] : 1;
		

		float k0 = f + (sumN(weight, sourcePhase, indegree, targetPhase          )/Cmean);
		float k1 = f + (sumN(weight, sourcePhase, indegree, targetPhase+k0*0.5f*h)/Cmean);
		float k2 = f + (sumN(weight, sourcePhase, indegree, targetPhase+k1*0.5f*h)/Cmean);
		float k3 = f + (sumN(weight, sourcePhase, indegree, targetPhase+k2*h     )/Cmean);

		//! \todo ensure neuron is valid
		//! \todo use precomputed factor and multiply
		targetPhase += h*(k0 + 2*k1 + 2*k2 + k3)/6.0f;
		phase1[n] = fmodf(targetPhase, 2.0f*M_PI) + (targetPhase < 0.0f ? 2.0f*M_PI: 0.0f);
	}
}

Foam::tmp<Foam::fvScalarMatrix>
Foam::diameterModels::IATEsources::phaseChange::R
(
    const volScalarField& alphai,
    volScalarField& kappai
) const
{
    if (!iDmdtPtr_)
    {
        iDmdtPtr_ = &alphai.mesh().lookupObject<volScalarField>
        (
            IOobject::groupName("iDmdt", pairName_)
        );
    }

    const volScalarField& iDmdt = *iDmdtPtr_;

    return -fvm::SuSp
    (
        (1.0/3.0)
       *iDmdt()
       /(alphai()*phase().rho()()),
        kappai
    );
}

void WaveTable::process(em_audio_stream_t buf, em_integer buflen)
{

    curFreq=frequency*transpose;
    float ph=phase()*mWaveTableSize/2;
    for(int i=0;i<buflen;i++)
    {



        float tablepos=(((curTablePos+ph)>=mWaveTableSize)?(curTablePos+ph-mWaveTableSize):(curTablePos+ph));
        tablepos=((tablepos<0)?(mWaveTableSize+tablepos):tablepos);

        int intPos=tablepos;
        float delta=tablepos-intPos;
        buf[i]+=volume()*(mWaveTable[intPos]+delta*(mWaveTable[(intPos==mWaveTableSize-1)?0:(intPos+1)]-mWaveTable[intPos])
                +delta*delta*(mWaveTable[(intPos==0)?(mWaveTableSize-1):(intPos-1)]+mWaveTable[(intPos==mWaveTableSize-1)?0:(intPos+1)]-2*mWaveTable[intPos]));
        //buf[i]*=mWindow[i];
        //buffer[i]*=ampl;
        curTablePos+=curFreq/BASE_FREQ;

        if(curTablePos>=mWaveTableSize)
        {
            curTablePos=0;

        }

    }
    SoundGenerator::process(buf,buflen);
}

    // Transform to a single k point:
    std::complex<double> XTable::kval(double kx, double ky) const 
    {
        check_array();
        // Don't evaluate if k not in fundamental period 
        kx*=_dx; ky*=_dx;
#ifdef FFT_DEBUG
        if (std::abs(kx) > M_PI || std::abs(ky) > M_PI) 
            throw FFTOutofRange("XTable::kval() args out of range");
#endif
        std::complex<double> I(0.,1.);
        std::complex<double> dxphase=std::exp(-I*kx);
        std::complex<double> dyphase=std::exp(-I*ky);
        std::complex<double> phase(1.,0.);
        std::complex<double> z;
        std::complex<double> sum=0.;

        const double* zptr=_array.get();
        std::complex<double> yphase=std::exp(I*(ky*_N/2));
        for (int iy=0; iy< _N; iy++) {
            phase = yphase;
            phase *= std::exp(I*(kx*_N/2));
            for (int ix=0; ix< _N ; ix++) {
                sum += phase* (*(zptr++));
                phase *= dxphase;
            }
            yphase *= dyphase;
        }
        sum *= _dx*_dx;
        return sum;
    }

/**
 * Build a tree using a given string
 * @param txt the text to build it from
 * @return the finished tree's root
 */
node *build_tree( char *txt )
{
    // init globals
    e = 0;
    root=NULL;
    f=NULL;
    current=NULL;
    memset( &last, 0, sizeof(pos) );
    memset( &old_beta, 0, sizeof(pos) );
    old_j = 0;
    str = txt;
    slen = strlen(txt);
    // actually build the tree
    root = node_create( 0, 0 );
    if ( root != NULL )
    {
        f = node_create_leaf( 0 );
        if ( f != NULL )
        {
            int i;
            node_add_child( root, f );
            for ( i=1; i<=slen; i++ )
                phase(i);
            set_e( root );
        }
    }
    return root;
}

void CheMPS2::Heff::addDiagonal2f3spin1(const int ikappa, double * memHeffDiag, const Sobject * denS, TensorOperator * Dtensor) const{

   int N2 = denS->gN2(ikappa);
   
   if (N2==1){

      int theindex = denS->gIndex();
      int ptr = denS->gKappa2index(ikappa);
      
      int NL = denS->gNL(ikappa);
      int TwoSL = denS->gTwoSL(ikappa);
      int IL = denS->gIL(ikappa);
      int NR = denS->gNR(ikappa);
      int TwoSR = denS->gTwoSR(ikappa);
      int IR = denS->gIR(ikappa);
      int TwoJ = denS->gTwoJ(ikappa);
      int N1 = denS->gN1(ikappa);
      
      int dimL     = denBK->gCurrentDim(theindex,  NL,TwoSL,IL);
      int dimR     = denBK->gCurrentDim(theindex+2,NR,TwoSR,IR);
      
      int fase = phase(TwoSR + TwoSL + 2*TwoJ + ((N1==1)?1:0) + 1);
      const double alpha = fase * (TwoJ+1) * sqrt(3.0*(TwoSR+1)) * Wigner::wigner6j(TwoJ,TwoJ,2,1,1,((N1==1)?1:0)) * Wigner::wigner6j(TwoJ,TwoJ,2,TwoSR,TwoSR,TwoSL);
      
      double * Dblock = Dtensor->gStorage(NR,TwoSR,IR,NR,TwoSR,IR);
      for (int cntR=0; cntR<dimR; cntR++){
         for (int cntL=0; cntL<dimL; cntL++){
            memHeffDiag[ptr + cntL + dimL*cntR] += alpha * Dblock[(dimR+1)*cntR];
         }
      }

   }
   
}

double Animation::calculateTimeToEffectChange(bool forwards, double localTime, double timeToNextIteration) const
{
    const double start = startTimeInternal() + specifiedTiming().startDelay;
    const double end = start + activeDurationInternal();

    switch (phase()) {
    case PhaseBefore:
        ASSERT(start >= localTime);
        return forwards
            ? start - localTime
            : std::numeric_limits<double>::infinity();
    case PhaseActive:
        return 0;
    case PhaseAfter:
        ASSERT(localTime >= end);
        // If this Animation is still in effect then it will need to update
        // when its parent goes out of effect. We have no way of knowing when
        // that will be, however, so the parent will need to supply it.
        return forwards
            ? std::numeric_limits<double>::infinity()
            : localTime - end;
    default:
        ASSERT_NOT_REACHED();
        return std::numeric_limits<double>::infinity();
    }
}

/*
 * \return in-degree of the neuron
 */
unsigned
loadIncoming(const nemo::runtime::RCM& rcm,
		unsigned target,
		int cycle,
		float* phaseBase,
		size_t phaseStride,
		std::vector<float>& weight,
		std::vector<float>& sourcePhase)
{
	unsigned indegree = rcm.indegree(target);
	weight.resize(indegree);
	sourcePhase.resize(indegree);

	if(indegree) {
		unsigned si = 0U;
		const std::vector<size_t>& warps = rcm.warps(target);
		for(std::vector<size_t>::const_iterator wi = warps.begin();
				wi != warps.end(); ++wi) {

			const nemo::RSynapse* rsynapse_p = rcm.data(*wi);
			const float* weight_p = rcm.weight(*wi);

			for(unsigned ri=0; ri < rcm.WIDTH && si < indegree; ri++, si++) {
				weight[si] = weight_p[ri];
				sourcePhase[si] = phase(phaseBase, phaseStride, cycle-int(rsynapse_p[ri].delay-1))[rsynapse_p[ri].source];
			}
		}
	}
	return indegree;
}

vector<complex<double> > fDSQPSKModulator(vector<int> bitsIn,
	vector<int> goldseq, double phi)
{
	// imaginary unit
	complex<double> j(0.0,1.0);

	int inLength = bitsIn.size();
	vector<complex<double> > tmp(inLength/2 + inLength%2, 0);

	int outLength = tmp.size() * goldseq.size();
	vector<complex<double> > out(outLength, 0);

	// create new bit vector
	// move odd bits to imaginary part, even bits to real part
	for (int i = 0; i < inLength; i++)
	{
		if (!bitsIn[i]) bitsIn[i] = -1;
		complex<double> bit(bitsIn[i],0.0);
		if (i%2) tmp[i/2] += j*bit;
		else tmp[i/2] += bit;
	}

	complex<double> phase(cos(phi*2*3.14159/360),sin(phi*2*3.14159/360));

	// spread with pn-code and multiply with exp(j*phi)
	for (int i = 0; i < outLength; i++)
	{
		complex<double> pnVal(goldseq[i%goldseq.size()], 0.0);
		out[i] = phase * tmp[i/goldseq.size()] * pnVal;
	}
	return out;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// generate
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RawPcmData::Ptr AdditiveSynthesizer::generate( size_t length )
{
   RawPcmData* result = new RawPcmData( getSamplingInfo(), length );

   const std::vector< HarmonicInfo >& harmonicsInfo = getHarmonicsInfo();
   assert( harmonicsInfo.size() > 0 );

   std::vector< double > phase( harmonicsInfo.size() );
   for ( size_t iHarmonic = 0; iHarmonic < harmonicsInfo.size(); ++iHarmonic )
   {
      phase[ iHarmonic ] = harmonicsInfo[ iHarmonic ].getPhase();
   }

   for ( size_t iSample = 0; iSample < length; ++iSample )
   {
      double val = 0;
      for ( size_t iHarmonic = 0; iHarmonic < harmonicsInfo.size(); ++iHarmonic )
      {
         phase[ iHarmonic ] += harmonicsInfo[ iHarmonic ].getPhaseStep();
         val += harmonicsInfo[ iHarmonic ].getAmplitude() * sin( phase[ iHarmonic ] );
      }
      (*result)[ iSample ] = val * getCurrentSampleAmplitude();
      nextSample();
   }

   setPhase( fmod( phase[ 0 ], 2 * M_PI ) );

   return RawPcmData::Ptr( result );
}

cv::Mat NFringeStructuredLight::WrapPhase(vector<cv::Mat> fringeImages)
{
  Utils::AssertOrThrowIfFalse(fringeImages.size() == m_numberOfFringes, 
	"Invalid number of fringes passed into phase wrapper");

  // Should be the same size as our fringe images 
  // and floating point precision for decimal phase values
  cv::Mat phase(fringeImages[0].size(), CV_32F, 0.0f);

  for(int row = 0; row < phase.rows; ++row)
  {
	for(int col = 0; col < phase.cols; ++col)
	{
	  float sine = 0.0f;
	  float cosine = 0.0f;
	  for(int fringe = 0; fringe < m_numberOfFringes; ++fringe)
	  {
		sine += ( float( fringeImages[fringe].at<uchar>(row, col) ) / 255.0 ) * sin(2.0 * M_PI * float(fringe) / float(m_numberOfFringes));
		cosine += ( float( fringeImages[fringe].at<uchar>(row, col) ) / 255.0 ) * cos(2.0 * M_PI * float(fringe) / float(m_numberOfFringes));
	  }
	  phase.at<float>(row, col) = atan2(sine, cosine);
	}
  }

  return phase;
}

void SplashScreen::think()
{
    if (phase() == Phase::Done)
    {
        if (m_fileBrowser.m_view)
            m_fileBrowser.m_view.reset();
        return;
    }

    ModalWindow::think();
    if (m_fileBrowser.m_view)
        m_fileBrowser.m_view->think();

    if (m_openButt.m_view)
        m_openButt.m_view->think();

    if (m_newProjBind.m_deferPath.size())
    {
        Log.report(logvisor::Info, _S("Making project '%s'"), m_newProjBind.m_deferPath.c_str());
        m_vm.projectManager().newProject(m_newProjBind.m_deferPath);
        m_newProjBind.m_deferPath.clear();
    }
    else if (m_openProjBind.m_deferPath.size())
    {
        Log.report(logvisor::Info, _S("Opening project '%s'"), m_openProjBind.m_deferPath.c_str());
        m_vm.projectManager().openProject(m_openProjBind.m_deferPath);
        m_openProjBind.m_deferPath.clear();
    }
    else if (m_extractProjBind.m_deferPath.size())
    {
        Log.report(logvisor::Info, _S("Extracting game '%s'"), m_extractProjBind.m_deferPath.c_str());
        m_vm.projectManager().extractGame(m_extractProjBind.m_deferPath);
        m_extractProjBind.m_deferPath.clear();
    }
}

void ppu::reg_write(const unsigned int addr,const unsigned char val, const unsigned int cycle) {
    //cout<<"!! "<<phases[phase(cycle)]<<" Write address: 0x"<<hex<<addr<<" Value: 0x"<<int(val)<<dec<<" Cycle: "<<cycle<<endl;
    ppu_phase curphase = phase(cycle);
    if(curphase == RENDER) ++s.writes_during_render;
    switch(addr) {
        case 0x2000:
        {
            int bg_table = control0.bg_table;
            int spr_table = control0.sprite_table;
            control0.reg = val;
            vram_ptr_reset.fields.page_y = (control0.base_addr & 2)/2;
            vram_ptr_reset.fields.page_x = (control0.base_addr & 1);

            int xscroll = (control0.base_addr%2) * 32 + vram_ptr_reset.fields.tile_x;
            int yscroll = (control0.base_addr/2) * 30 + vram_ptr_reset.fields.tile_y;
            if(curphase != RENDER && curphase != POST_RENDER) { //keeps track of where the screen will start
                set_scroll(xscroll*8+fine_x, yscroll*8+vram_ptr_reset.fields.fine_y);
                cart.ppu_change(cycle, BG_CHANGE_ADDR, control0.bg_table);
                cart.ppu_change(cycle, SPR_CHANGE_ADDR, control0.sprite_table);
                //if(control0.bg_table != bg_table) cout<<"Changed the bg table (to "<<control0.bg_table<<"), cycle: "<<cycle<<endl;
                //if(control0.sprite_table != spr_table) cout<<"Changed the sprite table (to "<<control0.sprite_table<<"), cycle: "<<cycle<<endl;
            } else { //keeps track of when (during rendering) the reset shifts, but only worry about the X reset
                set_scroll(xscroll*8+fine_x, yscroll*8+vram_ptr_reset.fields.fine_y, cycle%CLK_PER_LINE, cycle/CLK_PER_LINE);
                if(control0.bg_table != bg_table) {
                    //cout<<"Changed the bg table mid-frame! (to "<<control0.bg_table<<"), cycle: "<<cycle<<endl;
                    cart.ppu_change(cycle, BG_CHANGE_ADDR, 1 - control0.bg_table);
                }
                if(control0.sprite_table != spr_table) {
                    //cout<<"Changed the sprite table mid-frame! (to "<<control0.sprite_table<<"), cycle: "<<cycle<<endl;
                    cart.ppu_change(cycle, SPR_CHANGE_ADDR, 1 - control0.sprite_table);
                }
            }
        }
            break;
        case 0x2001:
            /*
            cout<<"Wrote "<<hex<<val<<dec<<" to 2001 (Color disabled: "<<control1.bits.color_disabled
                <<" show left bg: "<<control1.bits.show_left_bg<<" show left sprites: "<<control1.bits.show_left_sprites
                <<" show bg: "<<control1.bits.show_bg<<" show sprites: "<<control1.bits.show_sprites<<endl;
            */
            control1.reg = val;
            break;
        case 0x2003:
            spraddr = val;
            break;
        case 0x2004:
            if(spraddr/4 == 0) recalc_s0 = true;
            recalc_overflow = true;
            sprites[spraddr/4].val[spraddr%4] = val;
            spraddr++;
            break;
        case 0x2005:
            if(!scroll_latch) {
                vram_ptr_reset.fields.tile_x = (val & 0xf8)>>(3) ;
                fine_x = (val & 0x7);
                //cout<<"Set X reset to ("<<(control0.base_addr&1)*256+vram_ptr_reset.fields.tile_x*8+fine_x<<", "
                //    <<(control0.base_addr&2)*120+vram_ptr_reset.fields.tile_y*8+vram_ptr_reset.fields.fine_y<<") at cycle "
                //    <<cycle<<" ("<<cycle/CLK_PER_LINE<<", "<<cycle%CLK_PER_LINE<<") ("<<phases[curphase]<<")"<<endl;
            }
            else {

// Transform to a single k point:
DComplex
xTable::kval(double kx, double ky) const {
  // check this:  don't evaluate if x not in fundamental period:
  kx*=dx; ky*=dx;
  if (kx > 2*PI || ky > 2*PI) throw FFTOutofRange();
  DComplex I(0.,1.);
  DComplex dxphase=exp(-I*kx);
  DComplex dyphase=exp(-I*ky);
  DComplex phase(1.,0.);
  DComplex z;
  DComplex sum=0.;

  double *zptr=array;
  DComplex yphase=exp(I*(ky*N/2));
  for (int i=0; i< N; i++) {
    phase = yphase;
    phase *= exp(I*(kx*N/2));
    for (int j=0; j< N ; j++) {
      sum += phase* (*(zptr++));
      phase *= dxphase;
    }
    yphase *= dyphase;
  }
  sum *= dx*dx*scaleby;
  return sum;
}