/*
 Copyright (C) 2000-2004

 Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane

 This file is part of xmds.
 
 This program is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License
 as published by the Free Software Foundation; either version 2
 of the License, or (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

/*
  $Id: xmdsfield.cc,v 1.32 2005/04/27 05:44:19 joehope Exp $
*/

/*! @file xmdsfield.cc
  @brief Field element parsing classes and methods

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsFieldGeometry
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

// ******************************************************************************
void xmdsFieldGeometry::clear() {
  if(debugFlag) {
    printf("xmdsFieldGeometry::clear\n");
  }

  myDimensionsList.clear();
};

// ******************************************************************************
unsigned long xmdsFieldGeometry::nDims() const {
  if(debugFlag) {
    printf("xmdsFieldGeometry::nDims\n");
  }

  return myDimensionsList.size();
};

// ******************************************************************************
const dimensionStruct* xmdsFieldGeometry::dimension(
						    const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsFieldGeometry::dimension\n");
  }

  if(index>=myDimensionsList.size()) {
    throw xmdsException("Internal range error in xmdsFieldGeometry::dimension");
  }

  list<dimensionStruct>::const_iterator pdimensionStruct = myDimensionsList.begin();

  for(unsigned long i=0; i<index; i++) {
    pdimensionStruct++;
  }

  return &*pdimensionStruct;
};

// ******************************************************************************
void xmdsFieldGeometry::setDimension(
				     const unsigned long& index,
				     const dimensionStruct& newDimension) {
  if(debugFlag) {
    printf("xmdsFieldGeometry::setDimension\n");
  }

  if(index>myDimensionsList.size()) {
    throw xmdsException("Internal range error in xmdsFieldGeometry::setDimension");
  }

  if(index==myDimensionsList.size()) {
    myDimensionsList.push_back(newDimension);
  }

  list<dimensionStruct>::iterator pdimensionStruct = myDimensionsList.begin();

  for(unsigned long i=0; i<index; i++) {
    pdimensionStruct++;
  }

  *pdimensionStruct=newDimension;
};

// ******************************************************************************
void xmdsFieldGeometry::addDimension(
				     const dimensionStruct& newDimension) {
  if(debugFlag) {
    printf("xmdsFieldGeometry::addDimension\n");
  }

  myDimensionsList.push_back(newDimension);
};

// ******************************************************************************
bool xmdsFieldGeometry::getDimNumber(
				     const XMLString& dimName,
				     unsigned long& dimNumber) const {
  if(debugFlag) {
    printf("xmdsFieldGeometry::getDimNumber\n");
  }

  unsigned long i = 0;
  for(list<dimensionStruct>::const_iterator pDim = myDimensionsList.begin(); pDim != myDimensionsList.end(); pDim++) {
    if(pDim->name == dimName) {
      dimNumber=i;
      return 1;
    }
    i++;
  }

  return 0;
};

// ******************************************************************************
unsigned long xmdsFieldGeometry::fullSpace() const {
  if(debugFlag) {
    printf("xmdsFieldGeometry::fullSpace\n");
  }

  unsigned long two2n=1;
  for(unsigned long i=0; i<myDimensionsList.size(); i++) {
    two2n *= 2;
  }

  return two2n-1;
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsField Public
// ******************************************************************************
// ******************************************************************************

long nxmdsFields=0;  //!< Number of xmds field objects

// ******************************************************************************
xmdsField::xmdsField(
		     const xmdsSimulation *const yourSimulation,
		     const bool& yourVerboseMode) :
  xmdsElement(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsFields++;
    printf("xmdsField::xmdsField\n");
    printf("nxmdsFields=%li\n",nxmdsFields);
  }
};

// ******************************************************************************
xmdsField::~xmdsField() {
  if(debugFlag) {
    nxmdsFields--;
    printf("xmdsField::~xmdsField\n");
    printf("nxmdsFields=%li\n",nxmdsFields);
  }

  // destroy xmdsVectors, since they are not managed by the xmdsElement class

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    delete (*ppxmdsVector);
  }
};

// ******************************************************************************
void xmdsField::processElement(
			       const Element *const yourElement) {
  if(debugFlag) {
    printf("xmdsField::processElement\n");
  }

  if(verbose()) {
    printf("Processing field element ...\n");
  }

  myNeedsFFTWPlans=0;

  list<XMLString> anXMLStringList;
  list<XMLString>::const_iterator pXMLString;

  // ************************************
  // find name

  getAssignmentStrings(yourElement,"name",0,1,anXMLStringList);

  if(anXMLStringList.size()==1) {
    myName = *anXMLStringList.begin();
    if(verbose()) {
      printf("field name is '%s'\n",myName.c_str());
    }
  }
  else {
    myName = "main";
    printf("field name defaulting to '%s'\n",myName.c_str());
  }

  list<unsigned long> aLatticeList;
  list<unsigned long>::const_iterator pULong;
  list<domainStruct> aDomainList;
  list<domainStruct>::const_iterator pDomain;

  // ************************************
  // find transverse dimensions

  getAssignmentStrings(yourElement,"dimensions",0,0,anXMLStringList);

  if(anXMLStringList.size()>63) {
    throw xmdsException(yourElement,"xmds currently limited to 63 transverse dimensions");
  }

  if(anXMLStringList.size()>0) {

    unsigned long i;
	
    // ************************************
    // find transverse lattice

    getAssignmentULongs(yourElement,"lattice",1,anXMLStringList.size(),aLatticeList);
    i=0;
    for(pULong = aLatticeList.begin(); pULong != aLatticeList.end(); pULong++) {

      if(*pULong < 2) {
	sprintf(errorMessage(),	"lattice for dimension #%li must be > 1",i+1);
	throw xmdsException(yourElement,errorMessage());
      }

      if(verbose()) {
	printf("transverse lattice #%li  has %li points\n",i+1,*pULong);
      }

      i++;
    }

    // ************************************
    // find transverse domains

    getAssignmentDomains(yourElement,"domains",1,anXMLStringList.size(),aDomainList);

    if(verbose()) {
      i=0;
      for(pDomain = aDomainList.begin(); pDomain != aDomainList.end(); pDomain++) {
	printf("transverse domain #%li:\n",i+1);
	printf(" 	begin = %f\n",pDomain->begin);
	printf(" 	end   = %f\n",pDomain->end);
	i++;
      }
    }
  }
  else {
    if(verbose()) {
      printf("No transverse dimensions found.\n");
    }
  }

  // build field geometry dimensions

  myGeometry.clear();

  pXMLString = anXMLStringList.begin();
  pULong = aLatticeList.begin();
  pDomain = aDomainList.begin();

  for(unsigned long i=0; i< anXMLStringList.size(); i++) {

    dimensionStruct nextDimension;

    for(unsigned long j=0;j<i;j++) {
      if(*pXMLString==myGeometry.dimension(j)->name) {
	sprintf(errorMessage(),"duplicate dimension '%s'",pXMLString->c_str());
	throw xmdsException(yourElement,errorMessage());
      }
    }
    if(verbose()) {
      printf("adding transverse dimension '%s'\n",pXMLString->c_str());
    }

    nextDimension.name = *pXMLString;
    nextDimension.lattice = *pULong;
    nextDimension.domain = *pDomain;

    myGeometry.addDimension(nextDimension);

    pXMLString++;
    pULong++;
    pDomain++;
  }

  // ************************************
  // find samples

  getAssignmentULongs(yourElement,"samples",1,0,mySamplesList);

  // ************************************
  // find and process vector elements
  // ************************************

  const NodeList* candidateElements;

  candidateElements = yourElement->getElementsByTagName("vector",0);

  if(candidateElements->length()==0) {
    throw xmdsException(yourElement,"at least one vector expected");
  }

  for(unsigned long i=0; i<candidateElements->length(); i++) {

    const Element* nextElement = dynamic_cast<const Element*>(candidateElements->item(i));

    xmdsElement* newxmdsVectorElement = createxmdsVectorElement();

    newxmdsVectorElement->processElement(nextElement);
  }
};

// ******************************************************************************
void xmdsField::processVectors(
			       const list<XMLString>& vectorNamesList,
			       const unsigned long& space) const {
  if(debugFlag) {
    printf("xmdsField::processVectors\n");
  }

  for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {

    const xmdsVector* nextVector;

    if(!getVector(*pXMLString,nextVector)) {
      sprintf(errorMessage(),"Vector '%s' unknown",pXMLString->c_str());
      throw xmdsException(errorMessage());
    }

    for(list<XMLString>::const_iterator pXMLString2 = vectorNamesList.begin(); pXMLString2 != pXMLString; pXMLString2++) {
      if(*pXMLString == *pXMLString2) {
	sprintf(errorMessage(),"Duplicate vector '%s'",pXMLString->c_str());
	throw xmdsException(errorMessage());
      }
    }

    if(nextVector->initialSpace() != space){

      if(nextVector->vectorType() != COMPLEX) {

	const char* spaceDescription;

	if(space == 0) {
	  spaceDescription = "x-space";
	}
	else if(space == myGeometry.fullSpace()) {
	  spaceDescription = "k-space";
	}
	else {
	  spaceDescription = "mixed fourier space";
	}

	sprintf(errorMessage(),"Cannot convert vector '%s' to %s since it is not of type 'complex'",
		pXMLString->c_str(),spaceDescription);
	throw xmdsException(errorMessage());
      }

      myNeedsFFTWPlans = 1;

      nextVector->setNeedsFFTWRoutines();
    }

//  Always write FFTW Plans for main vector 
//  (otherwise a rare exception for deterministic MPI which uses FFTW to split field across processes)
//  Could make this MPI specifc, or else define the appropriate

    if(pXMLString->c_str()=="main")
      myNeedsFFTWPlans = 1;


    if(verbose()) {
      printf("adding vector '%s' to accessible vectors list\n",pXMLString->c_str());
    }
  }
};

// ******************************************************************************
void xmdsField::outputSampleCount() const {
  if(debugFlag) {
    printf("xmdsField::outputSampleCount\n");
  }

  if(mySamplesList.size() != simulation()->output()->nMomentGroups()) {
    throw xmdsException("number of samples in field element does not match number of output moment groups");
  }

  simulation()->output()->addSamples(mySamplesList);
};

// ******************************************************************************
bool xmdsField::needsFFTWPlans() const {
  if(debugFlag) {
    printf("xmdsField::needsFFTWPlans\n");
  }

  return myNeedsFFTWPlans;
};

// ******************************************************************************
void xmdsField::writeInitialisationCalls(
					 FILE *const outfile,
					 const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::writeInitialisationCalls\n");
  }

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writeInitialisationCall(outfile,indent);
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writeSampleCalls(
				 FILE *const outfile,
				 const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::writeSampleCalls\n");
  }

  unsigned long i=0;
  for(list<unsigned long>::const_iterator pULong = mySamplesList.begin(); pULong != mySamplesList.end(); pULong++) {
    if(*pULong>0) {
      fprintf(outfile,"%s_mg%li_sample();\n",indent,i);
    }
    i++;
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::assignActiveVectorPointers(
					   FILE *const outfile,
					   const char *const tempVectorName) const {
  if(debugFlag) {
    printf("xmdsField::assignActiveVectorPointers\n");
  }

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    fprintf(outfile,"_active_%s_%s = %s_%s_%s;\n",
	    myName.c_str(),(*ppxmdsVector)->name()->c_str(),tempVectorName,myName.c_str(),(*ppxmdsVector)->name()->c_str());
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
const XMLString* xmdsField::name() const {
  if(debugFlag) {
    printf("xmdsField::name\n");
  }
	
  return &myName;
};

// ******************************************************************************
const xmdsFieldGeometry* xmdsField::geometry() const{
  if(debugFlag) {
    printf("xmdsField::geometry\n");
  }

  return &myGeometry;
};

// ******************************************************************************
bool xmdsField::getVector(
			  const XMLString& vectorName,
			  const xmdsVector*& theVector) const {
  if(debugFlag) {
    printf("xmdsField::getVector\n");
  }

  theVector=0;

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    if(*(*ppxmdsVector)->name()==vectorName) {
      theVector = *ppxmdsVector;
      return 1;
    }
  }
  return 0;
};

// ******************************************************************************
void xmdsField::vectorNames(
			    list<XMLString>& vectorNamesList) const {
  if(debugFlag) {
    printf("xmdsField::vectorNames\n");
  }

  vectorNamesList.clear();

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    vectorNamesList.push_back(*(*ppxmdsVector)->name());
  }
};

// ******************************************************************************
void xmdsField::vectors2space(
			      FILE *const outfile,
			      const unsigned long& space,
			      const list<XMLString>& vectorNamesList,
			      const char *const indent) const {
  if(debugFlag) {
    printf("xmdsField::vectors2space\n");
  }

  if(myGeometry.nDims()==0) {
    return;
  }

  for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {

    const xmdsVector* nextVector;

    if(!getVector(*pXMLString,nextVector)) {
      throw xmdsException("Internal error in xmdsField::vectors2space: cannot find vector");
    }

    if(nextVector->needsFFTWRoutines()) {
      fprintf(outfile,"%s_%s_%s_go_space(%li);\n",indent,myName.c_str(),pXMLString->c_str(),space);
      fprintf(outfile,"\n");
    }
  }
};

// ******************************************************************************
void xmdsField::openLoops(
			  FILE *const outfile,
			  unsigned long space,
						  const list<XMLString>& vectorNamesList) const {
	if(debugFlag) {
		printf("xmdsField::openLoops\n");
	}
	
	for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
		fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",myName.c_str(),pXMLString->c_str());
	}
	fprintf(outfile,"\n");
	
	bool swapped=false;
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
		if((space&1)&((space>>1)&1)) {
			swapped=true;
		}
	}
	
	if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)&swapped) {
		fprintf(outfile,
				"double k%s = local_y_start_after_transpose*_%s_dk1;\n"
				"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n"
				"	 k%s -= _%s_lattice1*_%s_dk1;\n"
				"\n"    
				"for(long _i0=0; _i0<local_ny_after_transpose; _i0++) {\n\n",
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str(),
				myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());
		
		space >>= 1;
		
		fprintf(outfile,
				"	"
				"double k%s = 0;\n"
				"\n"
				"	"
				"for(long _i1=0; _i1<_%s_lattice0; _i1++) {\n"
				"\n",
				myGeometry.dimension(0)->name.c_str(),
				myName.c_str());
		
		space >>= 1;
		
		for(unsigned long i=2; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
    }
	else if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)) {
		if(space&1) {
			fprintf(outfile,"double k%s = local_x_start*_%s_dk0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str());
			fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",
					myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());	
			fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",
					myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
		}
		else
			fprintf(outfile,"double %s = _%s_xmin0+local_x_start*_%s_dx0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
		fprintf(outfile,"\n");
		
		fprintf(outfile,"for(long _i0=0; _i0<local_nx; _i0++) {\n");
		fprintf(outfile,"\n");
		
		space >>= 1;
		
		for(unsigned long i=1; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
    }
	else {
		for(unsigned long i=0; i<myGeometry.nDims(); i++) {
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			if(space&1) {
				fprintf(outfile,"double k%s = 0;\n",myGeometry.dimension(i)->name.c_str());
			}
			else {
				fprintf(outfile,"double %s = _%s_xmin%li;\n",myGeometry.dimension(i)->name.c_str(),myName.c_str(),i);
			}
			fprintf(outfile,"\n");
			
			for(unsigned long j=0; j<i; j++) {
				fprintf(outfile,"	");
			}
			fprintf(outfile,"for(long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,myName.c_str(),i,i);
			fprintf(outfile,"\n");
			
			space >>= 1;
		}
	}
						  };

// ******************************************************************************
void xmdsField::closeLoops(
			   FILE *const outfile,
			   unsigned long space,
			   const list<XMLString>& vectorNamesList) const {
  if(debugFlag) {
    printf("xmdsField::closeLoops\n");
  }

  if(myGeometry.nDims()>0) {
    fprintf(outfile,"\n");
    for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
      for(unsigned long j=0; j<myGeometry.nDims(); j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
	      myName.c_str(),pXMLString->c_str(),myName.c_str(),pXMLString->c_str());
    }
  }

  bool *isKSpace = new bool[myGeometry.nDims()];
  unsigned long tempSpace = space;
  for(long unsigned int i=0; i<myGeometry.nDims(); i++) {
	 if(tempSpace&1) {
		isKSpace[i] = true;
		}
	 else {
		isKSpace[i] = false;
		}
	tempSpace>>=1;
  }

  bool swapped=false;
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    if((space&1)&((space>>1)&1)) {
      swapped=true;
    }
  }

  if((!simulation()->parameters()->stochastic)&(simulation()->parameters()->usempi)&swapped) {
      for(unsigned long i=myGeometry.nDims(); i>2; i--) {
		  fprintf(outfile,"\n");
		  if(isKSpace[i-1]) {
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"k%s += _%s_dk%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1);
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li)\n",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);	
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);
		  }
		  else {
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"%s += _%s_dx%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1);
		  }
		  
		  for(unsigned long j=0; j<i; j++) {
			  fprintf(outfile,"	");
		  }
		  fprintf(outfile,"}\n");
	  }
	  
      fprintf(outfile,"\n");
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"k%s += _%s_dk0;\n",myGeometry.dimension(0)->name.c_str(),myName.c_str());
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",
			  myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());	
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",
			  myGeometry.dimension(0)->name.c_str(),myName.c_str(),myName.c_str());
	  
      for(unsigned long j=0; j<2; j++) {
		  fprintf(outfile,"	");
      }
      fprintf(outfile,"}\n");
	  
      fprintf(outfile,"\n");
	  
      fprintf(outfile,"	");
      fprintf(outfile,"k%s += _%s_dk1;\n",myGeometry.dimension(1)->name.c_str(),myName.c_str());
	  
      fprintf(outfile,"	");
      fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",
			  myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());	
	  
      fprintf(outfile,"	");
      fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",
			  myGeometry.dimension(1)->name.c_str(),myName.c_str(),myName.c_str());
	  
      fprintf(outfile,"	");
      fprintf(outfile,"}\n");
  }
  else
    {
      for(unsigned long i=myGeometry.nDims(); i>0; i--) {
		  
		  fprintf(outfile,"\n");
		  
		  if(isKSpace[i-1]) {
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"k%s += _%s_dk%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1);
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li)\n",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);	
			  
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",
					  myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1,myName.c_str(),i-1);
		  }
		  else {
			  for(unsigned long j=0; j<i; j++) {
				  fprintf(outfile,"	");
			  }
			  fprintf(outfile,"%s += _%s_dx%li;\n",myGeometry.dimension(i-1)->name.c_str(),myName.c_str(),i-1);
		  }
		  
		  for(unsigned long j=0; j<i; j++) {
			  fprintf(outfile,"	");
		  }
		  fprintf(outfile,"}\n");
      }
    }
	  
  delete isKSpace;

};

// ******************************************************************************
void xmdsField::writePlanCreationCalls(
				       FILE *const outfile,
				       const bool& useFFTWMeasure,
				       const bool& useWisdom) const {

  if(debugFlag) {
    printf("xmdsField::writePlanCreationCalls\n");
  }

  if(!myNeedsFFTWPlans) {
    return;
  }

  const char *useMeasureStr = "";
  const char *useWisdomStr = "";
  //const char *wisdomFname = simulation()->parameters()->simulationName.c_str();

  if(useFFTWMeasure) {
    useMeasureStr = "|FFTW_MEASURE";
  }

  if(useWisdom) {
    useWisdomStr = "|FFTW_USE_WISDOM";
  }

  for(unsigned long i=0;i<myGeometry.nDims();i++) {
    fprintf(outfile,"_fftw_lattice[%li]=_%s_lattice%li;\n",i,myName.c_str(),i);
  }
  if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)){
    if (useWisdom) {
      fprintf(outfile, 
          "{\n"
	      "\nFILE *fWisdom;\n"
	      "ifstream fIn;\n"
	      "ofstream fOut;\n"
	      "// this is a dodgy hack, but it should work\n"
	      "string findHostString, findHomeString, hostStuff, homeStuff, rmString, rankString;\n"
	      "ostringstream outString;\n"
	      "outString << rank;\n"
	      "rankString = outString.str();\n"
	      "hostStuff = \"host_rank_\" + rankString + \".stuff\";\n"
	      "findHostString = \"uname -n > \" + hostStuff;\n"
	      "system(findHostString.c_str());\n"
	      "homeStuff = \"home_rank_\" + rankString + \".stuff\";\n"
	      "findHomeString = \"echo $HOME > \" + homeStuff;\n"
	      "system(findHomeString.c_str());\n"
	      "string hostName, homeDir;\n"
	      "fIn.open(hostStuff.c_str());\n"
	      "if (fIn.fail()) {\n"
	      "    // do something\n"
	      "}\n"
	      "fIn >> hostName;\n"
	      "fIn.close();\n"
	      "fIn.open(homeStuff.c_str());\n"
	      "if (fIn.fail()) {\n"
	      "    // do something\n"
	      "}\n"
	      "fIn >> homeDir;\n"
	      "fIn.close();\n"
	      "rmString = \"rm \" + hostStuff + \" \" + homeStuff;\n"
	      "system(rmString.c_str());\n\n"
	      "string testFname, wisdomFname;\n"
	      "testFname = homeDir + \"/.xmds/wisdom/test-\" + rankString;\n"
	      "fOut.open(testFname.c_str());\n"
	      "if (fOut.fail()) {\n"
	      "    cout << \"Warning: ~/.xmds/wisdom directory doesn't seem to exist.\\n\";\n"
	      "    cout << \"Using current directory instead\\n\";\n"
	      "    wisdomFname = hostName + \".wisdom\";\n"
	      "}\n"
	      "else {\n"
	      "    fOut.close();\n"
	      "    rmString = \"rm \" + testFname;\n"
	      "    system(rmString.c_str());\n"
	      "    wisdomFname = homeDir + \"/.xmds/wisdom/\" + hostName + \".wisdom\";\n"
	      "}\n"
	      "printf(\"Performing fftw calculations\\n\");\n"
	      "if ((fWisdom = fopen(wisdomFname.c_str(), \"r\")) != NULL) {\n"
	      "    cout << \"Standing upon the shoulders of giants... (Importing wisdom)\\n\";\n"
	      "    fftw_import_wisdom_from_file(fWisdom);\n"
	      "    fclose(fWisdom);\n"
	      "}\n"
	      ); 
    }
    fprintf(outfile,"printf(\"Making forward plan\\n\");\n");
    fprintf(outfile,"_%s_forward_plan = fftwnd_mpi_create_plan(MPI_COMM_WORLD,_%s_ndims,_fftw_lattice,FFTW_FORWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
    // Using FFTW_TRANSPOSED_ORDER, so backward plan has different array sizes
    fprintf(outfile,"_fftw_lattice[0]=_%s_lattice1;\n",myName.c_str());
    fprintf(outfile,"_fftw_lattice[1]=_%s_lattice0;\n",myName.c_str());
    for(unsigned long i=2;i<myGeometry.nDims();i++) {
      fprintf(outfile,"_fftw_lattice[%li]=_%s_lattice%li;\n",i,myName.c_str(),i);
    }
    fprintf(outfile,"printf(\"Making backward plan\\n\");\n");
    fprintf(outfile,"_%s_backward_plan = fftwnd_mpi_create_plan(MPI_COMM_WORLD,_%s_ndims,_fftw_lattice,FFTW_BACKWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
    if (useWisdom) {
      fprintf(outfile,"printf(\"Keeping accumulated wisdom\\n\");\n");
      fprintf(outfile,"fWisdom = fopen(wisdomFname.c_str(), \"w\");\n");
      fprintf(outfile,"fftw_export_wisdom_to_file(fWisdom);\n");
      fprintf(outfile,"fclose(fWisdom);\n");
      fprintf(outfile,"}\n");
      fprintf(outfile,"printf(\"Finished fftw calculations\\n\");\n");
    }
  }
  else {
    if (useWisdom) {
      fprintf(outfile, 
          "{\n"
	      "FILE *fWisdom;\n"
	      "ifstream fIn;\n"
	      "ofstream fOut;\n"
	      "// this is a dodgy hack, but it should work\n"
	      "string findHostString, findHomeString, hostStuff, homeStuff, rmString, rankString;\n");
      if (!simulation()->parameters()->usempi) {
	fprintf(outfile, "int rank = 0;\n");
      }
      fprintf(outfile, 
	      "ostringstream outString;\n"
	      "outString << rank;\n"
	      "rankString = outString.str();\n"
	      "hostStuff = \"host_rank_\" + rankString + \".stuff\";\n"
	      "findHostString = \"uname -n > \" + hostStuff;\n"
	      "system(findHostString.c_str());\n"
	      "homeStuff = \"home_rank_\" + rankString + \".stuff\";\n"
	      "findHomeString = \"echo $HOME > \" + homeStuff;\n"
	      "system(findHomeString.c_str());\n"
	      "string hostName, homeDir;\n"
	      "fIn.open(hostStuff.c_str());\n"
	      "if (fIn.fail()) {\n"
	      "    // do something\n"
	      "}\n"
	      "fIn >> hostName;\n"
	      "fIn.close();\n"
	      "fIn.open(homeStuff.c_str());\n"
	      "if (fIn.fail()) {\n"
	      "    // do something\n"
	      "}\n"
	      "fIn >> homeDir;\n"
	      "fIn.close();\n"
	      "rmString = \"rm \" + hostStuff + \" \" + homeStuff;\n"
	      "system(rmString.c_str());\n\n"
	      "string testFname, wisdomFname;\n"
	      "testFname = homeDir + \"/.xmds/wisdom/test-\" + rankString;\n"
	      "fOut.open(testFname.c_str());\n"
	      "if (fOut.fail()) {\n"
	      "    cout << \"Warning: ~/.xmds/wisdom directory doesn't seem to exist.\\n\";\n"
	      "    cout << \"Using current directory instead\\n\";\n"
	      "    wisdomFname = hostName + \".wisdom\";\n"
	      "}\n"
	      "else {\n"
	      "    fOut.close();\n"
	      "    wisdomFname = homeDir + \"/.xmds/wisdom/\" + hostName + \".wisdom\";\n"
	      "    rmString = \"rm \" + testFname;\n"
	      "    system(rmString.c_str());\n"
	      "}\n"
	      "printf(\"Performing fftw calculations\\n\");\n"
	      "if ( (fWisdom = fopen(wisdomFname.c_str(), \"r\")) != NULL) {\n"
	      "    cout << \"Standing upon the shoulders of giants... (Importing wisdom)\\n\";\n"
	      "    fftw_import_wisdom_from_file(fWisdom);\n"
	      "    fclose(fWisdom);\n"
	      "}\n");
    }
    fprintf(outfile,"printf(\"Making forward plan\\n\");\n");
    fprintf(outfile,"_%s_forward_plan = fftwnd_create_plan(_%s_ndims,_fftw_lattice,FFTW_FORWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
    fprintf(outfile,"printf(\"Making backward plan\\n\");\n");
    fprintf(outfile,"_%s_backward_plan = fftwnd_create_plan(_%s_ndims,_fftw_lattice,FFTW_BACKWARD,FFTW_IN_PLACE%s%s);\n",myName.c_str(),myName.c_str(),useMeasureStr,useWisdomStr);
    if (useWisdom) {
      fprintf(outfile,
	      "printf(\"Keeping accumulated wisdom\\n\");\n"
	      "fWisdom = fopen(wisdomFname.c_str(),\"w\");\n"
	      "fftw_export_wisdom_to_file(fWisdom);\n"
	      "fclose(fWisdom);\n"
		  "}\n"
	      "printf(\"Finished fftw calculations\\n\");\n");
    }
  }
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writePlanDeletionCalls(
				       FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsField::writePlanDeletionCalls\n");
  }

  if(!myNeedsFFTWPlans) {
    return;
  }

  if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)){
    fprintf(outfile,"fftwnd_mpi_destroy_plan(_%s_forward_plan);\n",myName.c_str());
    fprintf(outfile,"fftwnd_mpi_destroy_plan(_%s_backward_plan);\n",myName.c_str());
  }
  else {
    fprintf(outfile,"fftwnd_destroy_plan(_%s_forward_plan);\n",myName.c_str());
    fprintf(outfile,"fftwnd_destroy_plan(_%s_backward_plan);\n",myName.c_str());
  }
  fprintf(outfile,"\n");
};

// ******************************************************************************
xmdsVector* xmdsField::createxmdsVector() {
  if(debugFlag) {
    printf("xmdsField::createxmdsVector\n");
  }

  xmdsVector* newVector = new xmdsVector(this);
  myVectorsList.push_back(newVector);
  return newVector;
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsField protected
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsField::writeDefines(
			     FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsField::writeDefines\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// field %s defines\n",myName.c_str());
  fprintf(outfile,"\n");

  fprintf(outfile,"#define _%s_ndims %li\n",myName.c_str(),myGeometry.nDims());

  for(unsigned long i=0; i<myGeometry.nDims(); i++) {
    const dimensionStruct* dimI = myGeometry.dimension(i);
    fprintf(outfile,"#define _%s_lattice%li %li\n",myName.c_str(),i,dimI->lattice);
    fprintf(outfile,"#define _%s_xmin%li %e\n",myName.c_str(),i,dimI->domain.begin);
    fprintf(outfile,"#define _%s_dx%li ((%e - %e)/(double)_%s_lattice%li)\n",myName.c_str(),i,dimI->domain.end,dimI->domain.begin,myName.c_str(),i);
    fprintf(outfile,"#define _%s_dk%li (2*M_PI/(%e - %e))\n",myName.c_str(),i,dimI->domain.end,dimI->domain.begin);
    fprintf(outfile,"#define d%s _%s_dx%li\n",dimI->name.c_str(),myName.c_str(),i);
    fprintf(outfile,"#define dk%s _%s_dk%li\n",dimI->name.c_str(),myName.c_str(),i);
  }

  fprintf(outfile,"\n");

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writeDefines(outfile);
  }
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writeGlobals(
			     FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsField::writeGlobals\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// field %s globals\n",myName.c_str());
  fprintf(outfile,"\n");

  if(myGeometry.nDims()>0) {
    fprintf(outfile,"const unsigned long _%s_size=_%s_lattice0",myName.c_str(),myName.c_str());
    for(unsigned long i=1; i<myGeometry.nDims(); i++) {
      fprintf(outfile,"*_%s_lattice%li",myName.c_str(),i);
    }
    fprintf(outfile,";\n");
  }
  else
    fprintf(outfile,"const unsigned long _%s_size=1;\n",myName.c_str());
  fprintf(outfile,"\n");

  if(myNeedsFFTWPlans) {
    if((!(simulation()->parameters()->stochastic))&(simulation()->parameters()->usempi)) {
	fprintf(outfile,"fftwnd_mpi_plan _%s_forward_plan;\n",myName.c_str());
	fprintf(outfile,"fftwnd_mpi_plan _%s_backward_plan;\n",myName.c_str());
      }
    else {
	fprintf(outfile,"fftwnd_plan _%s_forward_plan;\n",myName.c_str());
	fprintf(outfile,"fftwnd_plan _%s_backward_plan;\n",myName.c_str());
      }
  }

  fprintf(outfile,"\n");

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writeGlobals(outfile);
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writePrototypes(
				FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsField::writePrototypes\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// field %s prototypes\n",myName.c_str());
  fprintf(outfile,"\n");

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writePrototypes(outfile);
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::writeRoutines(
			      FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsField::writeRoutines\n");
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"//		field %s routines\n",myName.c_str());
  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"\n");

  for(list<const xmdsVector*>::const_iterator ppxmdsVector = myVectorsList.begin(); ppxmdsVector != myVectorsList.end(); ppxmdsVector++) {
    (*ppxmdsVector)->writeRoutines(outfile);
  }

  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsField::setName(
			const XMLString& yourName) {
  if(debugFlag) {
    printf("xmdsField::setName\n");
  }

  myName=yourName;
};

// ******************************************************************************
void xmdsField::setGeometry(
			    const xmdsFieldGeometry& yourGeometry) {
  if(debugFlag) {
    printf("xmdsField::setGeometry\n");
  }

  myGeometry = yourGeometry;
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsField private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
xmdsVectorElement* xmdsField::createxmdsVectorElement() {
  if( debugFlag) {
    printf("xmdsField::createxmdsVectorElement\n");
  }

  xmdsVectorElement* newVectorElement = new xmdsVectorElement(simulation(),verbose(),this);
  myVectorsList.push_back(newVectorElement);
  return newVectorElement;
};