/*
 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: xmdsmomentgroup.cc,v 1.43 2005/09/29 05:35:50 joehope Exp $
*/

/*! @file xmdsmomentgroup.cc
  @brief Moment group handling classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//				xmdsMomentGroup public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsMomentGroups=0;  //!< Number of xmds moment group objects

// ******************************************************************************
xmdsMomentGroup::xmdsMomentGroup(
				 const xmdsSimulation *const yourSimulation,
				 const bool& yourVerboseMode,
				 const unsigned long& yourGroupNumber) :
  xmdsField(yourSimulation,yourVerboseMode),
  myGroupNumber(yourGroupNumber) {
  if(debugFlag) {
    nxmdsMomentGroups++;
    printf("xmdsMomentGroup::xmdsMomentGroup\n");
    printf("nxmdsMomentGroups=%li\n",nxmdsMomentGroups);
  }
  nSamples=0;
};

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

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

  if(verbose()) {
    printf("Processing moment group %li ...\n",myGroupNumber);
  }

  const unsigned long nDims = simulation()->field()->geometry()->nDims();

  unsigned long i;
  const NodeList* candidateElements;

  myPost2MainDimList.clear();
  myMain2PostDimList.clear();

  myRequiresIntegrations=0;

  myPost2MainDimList.push_back(0);

  list<bool> rawSpaceList;
  rawSpaceList.push_back(0);

  dimensionStruct nextDimension;
  tempGeometry.addDimension(nextDimension);

  // ************************************
  // find sampling element

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

  if(candidateElements->length()==0) {
    throw xmdsException(yourElement,"<sampling> element required");
  }

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <sampling> elements defined");
  }

  const Element* samplingElement = dynamic_cast<const Element*> (candidateElements->item(0));

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

  // ************************************
  // get sampling field geometery

  if(nDims>0) {

    // ************************************
    // find sampling space

    list<bool> samplingSpaceList;
    getAssignmentBools(samplingElement,"fourier_space",0,simulation()->field()->geometry()->nDims(),samplingSpaceList);

    if(samplingSpaceList.size() == 0) {
      printf("Sampling space for moment group #%li defaulting to x-space.\n",myGroupNumber+1);
      mySamplingSpace = 0;
    }
    else {
      mySamplingSpace = spaceList2ULong(samplingSpaceList);

      if(verbose()) {
	for(i=0;i<nDims;i++) {
	  if(samplingSpace(i)) {
	    printf("sampling will be performed with dimension #%li in fourier space\n",i+1);
	  }
	  else {
	    printf("sampling will be performed with dimension #%li in normal space\n",i+1);
	  }
	}
      }
    }

    // ************************************
    // find sampling lattice

    getAssignmentULongs(samplingElement,"lattice",1,nDims,mySamplingLatticeList);

    list<unsigned long>::const_iterator pULong = mySamplingLatticeList.begin();
    for(i=0;i<nDims;i++) {

      const unsigned long fieldLatticeI = simulation()->field()->geometry()->dimension(i)->lattice;
      const unsigned long latticeI=*pULong;

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

      if(latticeI>1) {

	nextDimension.name = simulation()->field()->geometry()->dimension(i)->name;
	nextDimension.lattice = latticeI;
	nextDimension.domain = simulation()->field()->geometry()->dimension(i)->domain;

	tempGeometry.addDimension(nextDimension);

	myPost2MainDimList.push_back(i);
	rawSpaceList.push_back(samplingSpace(i));

	if(fieldLatticeI - latticeI*(fieldLatticeI/latticeI)) {
	  sprintf(errorMessage(),"moments group sampling lattice dimension[%li] does not divide evenly into same field lattice",i+1);
	  throw xmdsException(samplingElement,errorMessage());
	}

	if(verbose()) {
	  printf("sampling lattice dimension #%li divides same field lattice by %li\n",i+1,fieldLatticeI/latticeI);
	}
      }
      else if(latticeI==1) {

	// cross section
	if(verbose()) {
	  printf("transverse dimension #%li will be cross-sectioned at ",i+1);
	  if(samplingSpace(i)) {
	    printf("k = 0\n");
	  }
	  else {
	    printf("x = domain centre\n");
	  }
	}
      }

      else {
	// integrate
	myRequiresIntegrations = 1;

	if(verbose()) {
	  printf("transverse dimension #%li will be integrated in ",i+1);
	  if(samplingSpace(i)) {
	    printf(" fourier space\n");
	  }
	  else {
	    printf(" normal space\n");
	  }
	}
      }

      myMain2PostDimList.push_back(myPost2MainDimList.size()-1);
      pULong++;
    }
  }
  else
	{// No transverse dimensions, set sampling space to something safe
	mySamplingSpace=0;
	}

  myRawSpace=spaceList2ULong(rawSpaceList);

  // ************************************
  // find vectors

  getAssignmentStrings(samplingElement,"vectors",0,0,myVectorNamesList);

  if(myVectorNamesList.size()==0) {
    // no vectors specified, therefore assume using only main vector
    myVectorNamesList.push_back("main");
  }

  simulation()->field()->processVectors(myVectorNamesList,mySamplingSpace);

  // ************************************
  // find sampling moments

  getAssignmentStrings(samplingElement,"moments",1,0,mySamplingMomentsList);

  if(mySamplingMomentsList.size()==0) {
    throw xmdsException(samplingElement,"No sampling moments defined!");
  }
  if(verbose()) {
    for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
      printf("adding sampling moment '%s'\n",pXMLString->c_str());
    }
  }

  // ************************************
  // find sampling code

  mySamplingCode=*samplingElement->textContent(0);

  // check for non-white space charaters in code:

  if(mySamplingCode.isAllWhiteSpace()) {
    throw xmdsException(samplingElement,"No sampling code defined!");
  }

  if(verbose()) {
    printf("moment group sampling code loaded\n");
  }

  // ************************************
  // find post_propagation element

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

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <post_propagation> elements defined");
  }

  else if(candidateElements->length()==1) {

    // there is a post_propagation element

    const Element* postPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

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

    // ************************************
    // find post space

    list<bool> aSpaceList;

    getAssignmentBools(postPropagationElement,"fourier_space",0,myPost2MainDimList.size(),aSpaceList);

    if(aSpaceList.size() == 0) {
      printf("In moment group #%li, using sampling space for space of remaining post propagation dimensions.\n",myGroupNumber+1);
      for(i=0;i<myPost2MainDimList.size();i++) {
	aSpaceList.push_back(samplingSpace(post2MainDim(i)));
      }
    }

    myPostSpace = spaceList2ULong(aSpaceList);

    if(verbose()) {
      for(i=0;i<myPost2MainDimList.size();i++) {
	if((myPostSpace >> i)&1) {
	  printf("post propagation will be performed with remaining dimension #%li in fourier space\n",i+1);
	}
	else {
	  printf("post propagation will be performed with remaining dimension #%li in normal space\n",i+1);
	}
      }
    }

    // ************************************
    // find post moments

    getAssignmentStrings(postPropagationElement,"moments",1,0,myPostMomentsList);

    if(myPostMomentsList.size()==0) {
      throw xmdsException(postPropagationElement,"No post propagation moments defined!");
    }
    if(verbose()) {
      for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
	printf("post propagation moment '%s' added\n",pXMLString->c_str());
      }
    }

    // ************************************
    // find post code

    myPostCode=*postPropagationElement->textContent(0);

    // check for non-white space charaters in code:

    if(myPostCode.isAllWhiteSpace()) {
      throw xmdsException(postPropagationElement,"No post propagation code defined!");
    }

    if(verbose()) {
      printf("moment group post propagation code loaded\n");
    }
  }
  else {
    // no post_propagation element
    // need to assign across raw space and moments

    myPostSpace = myRawSpace;
    myPostMomentsList=mySamplingMomentsList;
  }

  char tempName[256];
  sprintf(tempName,"mg%li",myGroupNumber);

  setName((XMLString) tempName);

  // add moment group vectors

  // determine if raw sample vector needs to be complex

  complexRawVector = 0;

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

    const xmdsVector* nextVector;

    simulation()->field()->getVector(*pXMLString,nextVector);

    complexRawVector |= (nextVector->vectorType()==COMPLEX);
  }

  complexRawVector |= (myPostSpace != myRawSpace);

  xmdsVector* nextVector = createxmdsVector();
  nextVector->setName("raw");
  myPostProcessingVectorNamesList.push_back("raw");
  if(complexRawVector) {
    nextVector->setVectorType(COMPLEX);
  }
  else {
    nextVector->setVectorType(DOUBLE);
  }
  nextVector->setInitialSpace(myRawSpace);
  nextVector->setComponents(mySamplingMomentsList);

  nextVector = createxmdsVector();
  nextVector->setName("fullstep");
  myPostProcessingVectorNamesList.push_back("fullstep");
  nextVector->setVectorType(DOUBLE);
  nextVector->setInitialSpace(myPostSpace);
  nextVector->setComponents(myPostMomentsList);

  if(simulation()->parameters()->stochastic) {
    nextVector = createxmdsVector();
    nextVector->setName("fullstep_sd");
    myPostProcessingVectorNamesList.push_back("fullstep_sd");
    nextVector->setVectorType(DOUBLE);
    nextVector->setInitialSpace(myPostSpace);
    nextVector->setComponents(myPostMomentsList);
  }

  if(simulation()->parameters()->errorCheck) {
    nextVector = createxmdsVector();
    nextVector->setName("halfstep");
    myPostProcessingVectorNamesList.push_back("halfstep");
    nextVector->setVectorType(DOUBLE);
    nextVector->setInitialSpace(myPostSpace);
    nextVector->setComponents(myPostMomentsList);
    if(simulation()->parameters()->stochastic) {
      nextVector = createxmdsVector();
      nextVector->setName("halfstep_sd");
      myPostProcessingVectorNamesList.push_back("halfstep_sd");
      nextVector->setVectorType(DOUBLE);
      nextVector->setInitialSpace(myPostSpace);
      nextVector->setComponents(myPostMomentsList);
    }
  }

  processVectors(myPostProcessingVectorNamesList,myPostSpace);

};

// ******************************************************************************
void xmdsMomentGroup::addSamples(
				 const unsigned long& n2add) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::addSamples\n");
  }

  if(verbose())
    printf("Adding %li samples to moment group %li\n",n2add,myGroupNumber+1);

  nSamples += n2add;
};

// ******************************************************************************
void xmdsMomentGroup::finaliseGeometry() {
  if(debugFlag) {
    printf("xmdsMomentGroup::finaliseGeometry\n");
  }

  dimensionStruct firstDimension;

  firstDimension.name = simulation()->parameters()->propDimName;
  firstDimension.lattice = nSamples;
  firstDimension.domain.begin = 0;
  firstDimension.domain.end = 1;

  tempGeometry.setDimension(0,firstDimension);

  setGeometry(tempGeometry);
};

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

  return myRequiresIntegrations;
};

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

  return complexRawVector;
};

// ******************************************************************************
// ******************************************************************************
//				xmdsMomentGroup private
// ******************************************************************************
// ******************************************************************************

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

  if(verbose()) {
    printf("Writing moment group %li defines ...\n",myGroupNumber);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li defines\n",myGroupNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"#define _%s_ndims %li\n",mgFieldName,geometry()->nDims());
  fprintf(outfile,"#define _%s_lattice0 %li\n",mgFieldName,geometry()->dimension(0)->lattice);
  fprintf(outfile,"#define _%s_dx0 ((_%s_x0[_%s_lattice0-1] - _%s_x0[0])/(double)%li)\n",
	  mgFieldName,mgFieldName,mgFieldName,mgFieldName,geometry()->dimension(0)->lattice);
  fprintf(outfile,"#define _%s_dk0 (2*M_PI/(_%s_x0[_%s_lattice0-1] - _%s_x0[0]))\n",mgFieldName,mgFieldName,mgFieldName,mgFieldName);

  for(unsigned long i=1; i<geometry()->nDims(); i++) {

    const dimensionStruct* dimI = geometry()->dimension(i);

    fprintf(outfile,"#define _%s_lattice%li %li\n",mgFieldName,i,dimI->lattice);
    fprintf(outfile,"#define _%s_xmin%li %e\n",mgFieldName,i,dimI->domain.begin);
    fprintf(outfile,"#define _%s_dx%li ((%e - %e)/(double)%li)\n",mgFieldName,i,dimI->domain.end,dimI->domain.begin,dimI->lattice);
    fprintf(outfile,"#define _%s_dk%li (2*M_PI/(%e - %e))\n",mgFieldName,i,dimI->domain.end,dimI->domain.begin);
  }
  fprintf(outfile,"\n");

  fprintf(outfile,"#define _%s_raw_ncomponents %li\n",mgFieldName,(long)mySamplingMomentsList.size());
  fprintf(outfile,"#define _%s_fullstep_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());

  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"#define _%s_halfstep_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());
  }

  if(simulation()->parameters()->stochastic) {

    fprintf(outfile,"#define _%s_fullstep_sd_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());

    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"#define _%s_halfstep_sd_ncomponents %li\n",mgFieldName,(long)myPostMomentsList.size());
    }
  }
  fprintf(outfile,"\n");

};

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

  if(verbose()) {
    printf("Writing moment group %li globals ...\n",myGroupNumber);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li globals\n",myGroupNumber);
  fprintf(outfile,"\n");
  xmdsField::writeGlobals(outfile);
  fprintf(outfile,"unsigned long _%s_sample_pointer;\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"double *_%s_x0 = new double[_%s_lattice0];\n",mgFieldName,mgFieldName);
  fprintf(outfile,"\n");
};

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

  if(verbose()) {
    printf("Writing moment group %li prototypes ...\n",myGroupNumber);
  }

  const char* mgFieldName = name()->c_str();

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// moment group %li prototypes\n",myGroupNumber);
  fprintf(outfile,"\n");
  xmdsField::writePrototypes(outfile);
  fprintf(outfile,"void _%s_sample();\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _%s_process();\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _%s_write_out(\n",mgFieldName);
  fprintf(outfile,"	FILE*);\n");
  fprintf(outfile,"\n");
};

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

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const unsigned long myNDims = geometry()->nDims();
  const char* fieldName = simulation()->field()->name()->c_str();
  const char* mgFieldName = name()->c_str();

  unsigned long i;
  unsigned long j;

  if(verbose()) {
    printf("Writing moment group %s routines ...\n",mgFieldName);
  }

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"//		moment group %li routines\n",myGroupNumber);
  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"\n");

  xmdsField::writeRoutines(outfile);

  // ***************************************
  // ********  sample routine  *************
  // ***************************************

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_sample() {\n",mgFieldName);
  fprintf(outfile,"\n");

  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    if(complexRawVector) {
      fprintf(outfile,"complex %s;\n",pXMLString->c_str());
    }
    else {
      fprintf(outfile,"double %s;\n",pXMLString->c_str());
    }
  }
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,mySamplingSpace,myVectorNamesList,"");

  list<unsigned long>::const_iterator psamplingLatticeI = mySamplingLatticeList.begin();

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

  if(!myRequiresIntegrations&(simulation()->parameters()->stochastic|(!simulation()->parameters()->usempi))) {
    {
      fprintf(outfile,"unsigned long _%s_raw_index_pointer=_%s_sample_pointer",mgFieldName,mgFieldName);
      for(i=1;i<myNDims;i++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
      }
      fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);
      fprintf(outfile,"\n");
    }
  }

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      if(swapped) {
	  int firstlattice = *psamplingLatticeI;
	  psamplingLatticeI++;
	  int secondlattice = *psamplingLatticeI;
	  psamplingLatticeI++;
        
	  if(secondlattice==0) {
	    // integration over this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
            
	    fprintf(outfile,"double k%s = local_y_start_after_transpose*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	    fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);	
	    fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
                                    
	    fprintf(outfile,"\n");

	    fprintf(outfile,"for(long _i1=0; _i1<local_ny_after_transpose; _i1++) {\n");
	    fprintf(outfile,"\n");
	  }
	  else if(secondlattice==1) {
	    // cross-section in this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }

	    fprintf(outfile,"\n");

	    fprintf(outfile,"double k%s = 0;\n",simulation()->field()->geometry()->dimension(1)->name.c_str());
	    fprintf(outfile,"unsigned long _i1 = 0;\n");                    
	    fprintf(outfile,"if(local_y_start_after_transpose==0) {\n");                    
	  }
	  else {
	    // normal sampling
	    fprintf(outfile,"long first_x_pointer, last_x_pointer;\n");
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"double k%s;\n\n",simulation()->field()->geometry()->dimension(1)->name.c_str());
                    
	    fprintf(outfile,"if(local_y_start_after_transpose<(_%s_lattice%li-1)/2+1) {\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   first_x_pointer=local_y_start_after_transpose;\n");
	    fprintf(outfile,"   k%s=local_y_start_after_transpose*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"   }\n");
	    fprintf(outfile,"else if(local_y_start_after_transpose>(_%s_lattice%li-1)/2+_%s_lattice1-_%s_lattice%li) {\n"
		    ,mgFieldName,main2PostDim(1),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   first_x_pointer=local_y_start_after_transpose-_%s_lattice1+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   k%s=(local_y_start_after_transpose-_%s_lattice1)*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"   }\n");
	    fprintf(outfile,"else {\n");
	    fprintf(outfile,"   first_x_pointer=(_%s_lattice%li-1)/2+1;\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   k%s=((_%s_lattice%li-1)/2-_%s_lattice%li+1)*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),mgFieldName,main2PostDim(1),fieldName);
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		fprintf(outfile,"   _%s_%s_index_pointer=(first_x_pointer+_%s_lattice1-_%s_lattice%li-local_y_start_after_transpose)",fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(1));
		fprintf(outfile,"*_%s_lattice0",fieldName);
		for(i=2;i<nDims;i++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	    fprintf(outfile,"   }\n\n");
                    
	    fprintf(outfile,"if(local_y_start_after_transpose+local_ny_after_transpose-1<(_%s_lattice%li-1)/2+1)\n",mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   last_x_pointer=local_y_start_after_transpose+local_ny_after_transpose-1;\n");
	    fprintf(outfile,"else if(local_y_start_after_transpose+local_ny_after_transpose-1>(_%s_lattice%li-1)/2+_%s_lattice1-_%s_lattice%li)\n"
		    ,mgFieldName,main2PostDim(1),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"   last_x_pointer=local_y_start_after_transpose+local_ny_after_transpose-1-_%s_lattice1+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"else \n");
	    fprintf(outfile,"   last_x_pointer=(_%s_lattice%li-1)/2;\n\n",mgFieldName,main2PostDim(1));
    
	    fprintf(outfile,"for(long _i1=first_x_pointer;_i1<last_x_pointer+1;_i1++) {\n");
	    fprintf(outfile,"\n");
	  }
                

	  const char* dimName = simulation()->field()->geometry()->dimension(0)->name.c_str();

	  fprintf(outfile,"	");
	  fprintf(outfile,"double k%s = 0;\n",dimName);
	  fprintf(outfile,"\n");

	  if(firstlattice==0) {
	    // integration over this dimension
	    fprintf(outfile,"	");
	    fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",fieldName);
	    fprintf(outfile,"\n");
	  }
	  else if(firstlattice==1) {
	    // cross-section in this dimension
	    fprintf(outfile,"	");
	    fprintf(outfile,"unsigned long _i0 = 0;\n");
	  }
	  else {
	    // normal sampling
	    fprintf(outfile,"	");
	    fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_lattice%li;_i0++) {\n",mgFieldName,main2PostDim(0));
	    fprintf(outfile,"\n");
	  }
                

	  for(i=2;i<nDims;i++) {
    
            const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();
    
            for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
            if(samplingSpace(i)) {
	      fprintf(outfile,"double k%s = 0;\n",dimName);
	    }
            else {
	      fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	    }
            fprintf(outfile,"\n");
    
            if(*psamplingLatticeI==0) {
	      // integration over this dimension
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	"); 
	      }
	      fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
	      fprintf(outfile,"\n");
	    }
            else if(*psamplingLatticeI==1) {
	      // cross-section in this dimension
	      if(samplingSpace(i)) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	      }
	      else {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"%s += _i%li*_%s_dx%li;\n",dimName,i,fieldName,i);
    
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _i%li*_%s_%s_ncomponents",
			  fieldName,pXMLString->c_str(),i,fieldName,pXMLString->c_str());
		  for(j=i+1;j<nDims;j++) {
		    fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		  }
		  fprintf(outfile,";\n");
		}
	      }
	    }
            else {
	      // normal sampling
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,main2PostDim(i),i);
	      fprintf(outfile,"\n");
	    }
    
            psamplingLatticeI++;
	  }
        }
      else
        {
	  // Make the first one smaller.  Then loop.
	  // Add the k-space possibility to this first loop everywhere.

	  const char* dimName = simulation()->field()->geometry()->dimension(0)->name.c_str();

	  if(*psamplingLatticeI==0) {
	    // integration over this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
                
	    if(samplingSpace(0)) {
		fprintf(outfile,"double k%s = local_x_start*_%s_dk0;\n",dimName,fieldName);
		fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",dimName,fieldName,fieldName);	
		fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",dimName,fieldName,fieldName);
	      }
	    else {
	      fprintf(outfile,"double %s = _%s_xmin0 + local_x_start*_%s_dx0;\n",dimName,fieldName,fieldName);
	    }
                                    
	    fprintf(outfile,"\n");

	    fprintf(outfile,"for(long _i0=0; _i0<local_nx; _i0++) {\n");
	    fprintf(outfile,"\n");
	  }
	  else if(*psamplingLatticeI==1) {
	    // cross-section in this dimension
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      fprintf(outfile,"long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
	    }
	    if(!myRequiresIntegrations) {
		fprintf(outfile,"unsigned long _%s_raw_index_pointer=_%s_sample_pointer",mgFieldName,mgFieldName);
		for(i=1;i<myNDims;i++) {
		  fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		}
		fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);
		fprintf(outfile,"\n");
	      }

	    fprintf(outfile,"\n");

	    if(samplingSpace(0)) {
		fprintf(outfile,"double k%s = 0;\n",dimName);
		fprintf(outfile,"unsigned long _i0 = 0;\n");                    
		fprintf(outfile,"if(local_x_start==0) {\n");                    
	      }
	    else
	      {
		fprintf(outfile,"double %s = _%s_xmin0;\n",dimName,fieldName);
		fprintf(outfile,"\n");

		fprintf(outfile,"long _i0=_%s_lattice0/2-local_x_start;\n",fieldName);
    
		fprintf(outfile,"if((_i0+1>0)&(_i0<local_nx)) {\n");
    
		fprintf(outfile,"  %s += (_i0+local_x_start)*_%s_dx0;\n",dimName,fieldName);
    
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"  _%s_%s_index_pointer += _i0*_%s_%s_ncomponents",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		  for(j=1;j<nDims;j++) {
		    fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		  }
		  fprintf(outfile,";\n\n");
		}
	      }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(0)) {
		fprintf(outfile,"long first_x_pointer, last_x_pointer;\n");
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"long _%s_%s_index_pointer;\n",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"double k%s;\n\n",dimName);
                    
		fprintf(outfile,"if(local_x_start<(_%s_lattice%li-1)/2+1) {\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   first_x_pointer=local_x_start;\n");
		fprintf(outfile,"   k%s=local_x_start*_%s_dk0;\n",dimName,fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"   }\n");
		fprintf(outfile,"else if(local_x_start>(_%s_lattice%li-1)/2+_%s_lattice0-_%s_lattice%li) {\n"
			,mgFieldName,main2PostDim(0),fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   first_x_pointer=local_x_start-_%s_lattice0+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   k%s=(local_x_start-_%s_lattice0)*_%s_dk0;\n",dimName,fieldName,fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  fprintf(outfile,"   _%s_%s_index_pointer=0;",fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"   }\n");
		fprintf(outfile,"else {\n");
		fprintf(outfile,"   first_x_pointer=(_%s_lattice%li-1)/2+1;\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   k%s=((_%s_lattice%li-1)/2-_%s_lattice%li+1)*_%s_dk0;\n",dimName,mgFieldName,main2PostDim(0),mgFieldName,main2PostDim(0),fieldName);
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		    fprintf(outfile,"   _%s_%s_index_pointer=(first_x_pointer+_%s_lattice0-_%s_lattice%li-local_x_start)",fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(0));
		    for(i=1;i<nDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		    }
		    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
		  }
		fprintf(outfile,"   }\n\n");
                    
		fprintf(outfile,"if(local_x_start+local_nx-1<(_%s_lattice%li-1)/2+1)\n",mgFieldName,main2PostDim(0));
		fprintf(outfile,"   last_x_pointer=local_x_start+local_nx-1;\n");
		fprintf(outfile,"else if(local_x_start+local_nx-1>(_%s_lattice%li-1)/2+_%s_lattice0-_%s_lattice%li)\n"
			,mgFieldName,main2PostDim(0),fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"   last_x_pointer=local_x_start+local_nx-1-_%s_lattice0+_%s_lattice%li;\n",fieldName,mgFieldName,main2PostDim(0));
		fprintf(outfile,"else \n");
		fprintf(outfile,"   last_x_pointer=(_%s_lattice%li-1)/2;\n\n",mgFieldName,main2PostDim(0));
 
		if(!myRequiresIntegrations) {
		    fprintf(outfile,"unsigned long _%s_raw_index_pointer=_%s_sample_pointer",mgFieldName,mgFieldName);
		    for(i=1;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents + first_x_pointer",mgFieldName);
		    for(i=2;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);            
                
		    fprintf(outfile,"\n");
		  }
                       
		fprintf(outfile,"for(long _i0=first_x_pointer;_i0<last_x_pointer+1;_i0++) {\n");
		fprintf(outfile,"\n");
	      }
	    else
	      {
		fprintf(outfile,"long first_x_pointer=_%s_lattice1-(_%s_lattice1*(_%s_lattice0-local_x_start))/_%s_lattice0;\n"
			,mgFieldName,mgFieldName,fieldName,fieldName);
		fprintf(outfile,"long last_x_pointer=(_%s_lattice1*(local_x_start+local_nx-1))/_%s_lattice0;\n\n",mgFieldName,fieldName);
 
		if(!myRequiresIntegrations) {
		    fprintf(outfile,"unsigned long _%s_raw_index_pointer=_%s_sample_pointer",mgFieldName,mgFieldName);
		    for(i=1;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents + first_x_pointer",mgFieldName);
		    for(i=2;i<myNDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
		    }
		    fprintf(outfile,"*_%s_raw_ncomponents;\n",mgFieldName);            
                
		    fprintf(outfile,"\n");
		  }
                       
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		    fprintf(outfile,"unsigned long _%s_%s_index_pointer=(first_x_pointer*(_%s_lattice0/_%s_lattice1)-local_x_start)"
                            ,fieldName,pXMLString->c_str(),fieldName,mgFieldName);
		    for(i=1;i<nDims;i++) {
		      fprintf(outfile,"*_%s_lattice%li",fieldName,i);
		    }
		    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
		  }
		fprintf(outfile,"\n");
    
		fprintf(outfile,"double %s = _%s_xmin0 + first_x_pointer*_%s_dx1;\n",dimName,fieldName,mgFieldName);
		fprintf(outfile,"\n");
            
		fprintf(outfile,"for(long _i0=first_x_pointer;_i0<last_x_pointer+1;_i0++) {\n");
		fprintf(outfile,"\n");
	      }
        
	  }

	  psamplingLatticeI++;

	  for(i=1;i<nDims;i++) {
    
            const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();
    
            for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
            if(samplingSpace(i)) {
	      fprintf(outfile,"double k%s = 0;\n",dimName);
	    }
            else {
	      fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	    }
            fprintf(outfile,"\n");
    
            if(*psamplingLatticeI==0) {
	      // integration over this dimension
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
	      fprintf(outfile,"\n");
	    }
            else if(*psamplingLatticeI==1) {
	      // cross-section in this dimension
	      if(samplingSpace(i)) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	      }
	      else {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"%s += _i%li*_%s_dx%li;\n",dimName,i,fieldName,i);
    
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _i%li*_%s_%s_ncomponents",
			  fieldName,pXMLString->c_str(),i,fieldName,pXMLString->c_str());
		  for(j=i+1;j<nDims;j++) {
		    fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		  }
		  fprintf(outfile,";\n");
		}
	      }
	    }
            else {
	      // normal sampling
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,main2PostDim(i),i);
	      fprintf(outfile,"\n");
	    }
    
            psamplingLatticeI++;
	  }
        }
    } 
  else 
    {
      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	fprintf(outfile,"unsigned long _%s_%s_index_pointer=0;\n",fieldName,pXMLString->c_str());
      }
      fprintf(outfile,"\n");
    
      for(i=0;i<nDims;i++) {
    
        const char* dimName = simulation()->field()->geometry()->dimension(i)->name.c_str();

        for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
        if(samplingSpace(i)) {
	  fprintf(outfile,"double k%s = 0;\n",dimName);
	}
        else {
	  fprintf(outfile,"double %s = _%s_xmin%li;\n",dimName,fieldName,i);
	}
        fprintf(outfile,"\n");

        if(*psamplingLatticeI==0) {
	  // integration over this dimension
	  for(j=0;j<i;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,fieldName,i,i);
	  fprintf(outfile,"\n");
	}
        else if(*psamplingLatticeI==1) {
	  // cross-section in this dimension
	  if(samplingSpace(i)) {
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"unsigned long _i%li = 0;\n",i);
	  }
	  else {
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"unsigned long _i%li=_%s_lattice%li/2;\n",i,fieldName,i);
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"%s += _i%li*_%s_dx%li;\n",dimName,i,fieldName,i);

	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _i%li*_%s_%s_ncomponents",
		      fieldName,pXMLString->c_str(),i,fieldName,pXMLString->c_str());
	      for(j=i+1;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,";\n");
	    }
	  }
	}
        else {
	  // normal sampling
	  for(j=0;j<i;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,main2PostDim(i),i);
	  fprintf(outfile,"\n");
	}

        psamplingLatticeI++;
      }
    }
    
  fprintf(outfile,"// ************ Sampling code! ******************\n");
  fprintf(outfile,"%s\n",mySamplingCode.c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  if(myRequiresIntegrations|swapped) {
    for(j=0;j<nDims;j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"unsigned long _%s_raw_index_pointer=(_%s_sample_pointer",mgFieldName,mgFieldName);
    for(i=1;i<myNDims;i++) {
      fprintf(outfile,"*_%s_lattice%li",mgFieldName,i);
    }
    for(i=1;i<myNDims;i++) {
      fprintf(outfile," + _i%li",post2MainDim(i));
      for(j=i+1;j<myNDims;j++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,j);
      }
    }
    fprintf(outfile,")*_%s_raw_ncomponents;\n",mgFieldName);
    fprintf(outfile,"\n");
  }

  i=0;
  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    for(j=0;j<nDims;j++) {
	fprintf(outfile,"	");
    }
    if(myRequiresIntegrations) {
      fprintf(outfile,"_%s_raw[_%s_raw_index_pointer+%li] += %s",mgFieldName,mgFieldName,i,pXMLString->c_str());
      j=0;
      psamplingLatticeI = mySamplingLatticeList.begin();
      while(j<nDims) {
	if(*psamplingLatticeI==0) {
	  if(samplingSpace(j)) {
	    fprintf(outfile,"*_%s_dk%li",fieldName,j);
	  }
	  else {
	    fprintf(outfile,"*_%s_dx%li",fieldName,j);
	  }
	}
	psamplingLatticeI++;
	j++;
      }
      fprintf(outfile,";\n");
    }
    else
      fprintf(outfile,"_%s_raw[_%s_raw_index_pointer+%li] = %s;\n",mgFieldName,mgFieldName,i,pXMLString->c_str());
    i++;
  }

  // close nested loops

  if(!myRequiresIntegrations&!swapped) {
    if(nDims>0) {
      for(i=0;i<nDims;i++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_raw_index_pointer += _%s_raw_ncomponents;\n",mgFieldName,mgFieldName);
      fprintf(outfile,"\n");
    }
  }

  list<unsigned long>::const_reverse_iterator psamplingLatticeIr = mySamplingLatticeList.rbegin();

  if(nDims>0) {
    if(swapped) {
        for(i=nDims;i>2;i--) {

	  const char* dimName = simulation()->field()->geometry()->dimension(i-1)->name.c_str();

	  fprintf(outfile,"\n");
	  if(*psamplingLatticeIr==0) {
	    // integration
	    if(i==nDims) {
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	    }
	    if(samplingSpace(i-1)) {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,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",dimName,fieldName,i-1,fieldName,i-1);	
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,fieldName,i-1,fieldName,i-1);
	    }
	    else {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,fieldName,i-1);
	    }
	  }
	  else if(*psamplingLatticeIr==1) {
	    // cross-section
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li",fieldName,pXMLString->c_str(),fieldName,i-1);
	      if(!samplingSpace(i-1)) {
		fprintf(outfile,"-_i%li",i-1);
	      }
	      if(!(i==nDims)) {
		fprintf(outfile,"-1");
	      }
	      fprintf(outfile,")");
	      for(j=i;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	    }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(i-1)) {
	      // narrow k-space window
	      if(i==nDims) {
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"\n");
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li) {\n",
		      dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice%li-_%s_lattice%li)",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	}\n");
	    }
	    else {
	      // coarse x-space sampling
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li/_%s_lattice%li",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		if(!(i==nDims)) {
		  fprintf(outfile,"-1");
		}
		fprintf(outfile,")");
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,mgFieldName,main2PostDim(i-1));
	    }
	  }

	  if(!(*psamplingLatticeIr==1)) {
	    // if not cross-section
	    for(j=0;j<i;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"}\n");
	  }

	  psamplingLatticeIr++;
        }

        psamplingLatticeIr++;
        
        const char* dimName = simulation()->field()->geometry()->dimension(0)->name.c_str();

        fprintf(outfile,"\n");
        if(*psamplingLatticeIr==0) {
	  // integration
	  if(2==nDims) {
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<2;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
		      fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
	  }
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk0;\n",dimName,fieldName);
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",dimName,fieldName,fieldName);	
	  for(unsigned long j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice0*_%s_dk0;\n",dimName,fieldName,fieldName);
	}
        else if(*psamplingLatticeIr==1) {
	  // cross-section
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    for(j=0;j<2;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice0",fieldName,pXMLString->c_str(),fieldName);
	    if(!(2==nDims)) {
	      fprintf(outfile,"-1");
	    }
	    fprintf(outfile,")");
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	}
        else {
	  // normal sampling
	  // narrow k-space window
	  if(2==nDims) {
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<2;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
		      fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	    }
	    fprintf(outfile,"\n");
	  }
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk0;\n",dimName,fieldName);
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk0) {\n",
		  dimName,mgFieldName,main2PostDim(0),fieldName);
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    for(j=0;j<2;j++) {
		fprintf(outfile,"	");
	    }
	    fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice0-_%s_lattice%li)",
		    fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(0));
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	  for(j=0; j<2; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk0;\n",dimName,mgFieldName,main2PostDim(0),fieldName);
	  for(j=0;j<2;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	}\n");
	}

        if(!(*psamplingLatticeIr==1)) {
	  // if not cross-section
	  for(j=0;j<2;j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"}\n");
	}

        psamplingLatticeIr--;
                
        fprintf(outfile,"\n");
        if(*psamplingLatticeIr==0) {
	  // integration
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"k%s +=	_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"if(k%s>((_%s_lattice1-1)/2 + 0.1)*_%s_dk1)\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);	
	  for(unsigned long j=0; j<1; j++) {
	      fprintf(outfile,"	");
	  }
	  fprintf(outfile,"	k%s -= _%s_lattice1*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName,fieldName);
	}
        else if(*psamplingLatticeIr==1) {
	  // cross-section
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    fprintf(outfile,"	");
	    fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice1",fieldName,pXMLString->c_str(),fieldName);
	    fprintf(outfile,"-1)");
	    fprintf(outfile,"*_%s_lattice0",fieldName);
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	}
        else {
	  // normal sampling
	  // narrow k-space window
	  fprintf(outfile,"	");
	  fprintf(outfile,"k%s +=	_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),fieldName);
	  fprintf(outfile,"	");
	  fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk1) {\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),fieldName);
	  for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	    fprintf(outfile,"	");
	    fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice1-_%s_lattice%li)",
		    fieldName,pXMLString->c_str(),fieldName,mgFieldName,main2PostDim(1));
	    fprintf(outfile,"*_%s_lattice0",fieldName);
	    for(j=2;j<nDims;j++) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	    }
	    fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	  }
	  fprintf(outfile,"	");
	  fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk1;\n",simulation()->field()->geometry()->dimension(1)->name.c_str(),mgFieldName,main2PostDim(1),fieldName);
	  fprintf(outfile,"	");
	  fprintf(outfile,"	}\n");
	}

        // even if cross-section
        fprintf(outfile,"	");
        fprintf(outfile,"}\n");
    }
    else {
        for(i=nDims;i>0;i--) {
        
	  const char* dimName = simulation()->field()->geometry()->dimension(i-1)->name.c_str();
        
	  fprintf(outfile,"\n");
	  if(*psamplingLatticeIr==0) {
	    // integration
	    if(i==nDims) {
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	    }
	    if(samplingSpace(i-1)) {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,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",dimName,fieldName,i-1,fieldName,i-1);	
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,fieldName,i-1,fieldName,i-1);
	    }
	    else {
	      for(unsigned long j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,fieldName,i-1);
	    }
	  }
	  else if(*psamplingLatticeIr==1) {
	    // cross-section
	    for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li",fieldName,pXMLString->c_str(),fieldName,i-1);
	      if(!samplingSpace(i-1)) {
		fprintf(outfile,"-_i%li",i-1);
	      }
	      if(!(i==nDims)) {
		fprintf(outfile,"-1");
	      }
	      fprintf(outfile,")");
	      for(j=i;j<nDims;j++) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,j);
	      }
	      fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	    }
	  }
	  else {
	    // normal sampling
	    if(samplingSpace(i-1)) {
	      // narrow k-space window
	      if(i==nDims) {
		for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		  for(j=0;j<i;j++) {
		      fprintf(outfile,"	");
		  }
		  fprintf(outfile,"_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
			  fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
		}
		fprintf(outfile,"\n");
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"k%s +=	_%s_dk%li;\n",dimName,fieldName,i-1);
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"if(k%s>((_%s_lattice%li-1)/2 + 0.1)*_%s_dk%li) {\n",
		      dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"	_%s_%s_index_pointer += (_%s_lattice%li-_%s_lattice%li)",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      for(j=0; j<i; j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",dimName,mgFieldName,main2PostDim(i-1),fieldName,i-1);
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"	}\n");
	    }
	    else {
	      // coarse x-space sampling
	      for(list<XMLString>::const_iterator pXMLString = myVectorNamesList.begin(); pXMLString != myVectorNamesList.end(); pXMLString++) {
		for(j=0;j<i;j++) {
		    fprintf(outfile,"	");
		}
		fprintf(outfile,"_%s_%s_index_pointer += (_%s_lattice%li/_%s_lattice%li",
			fieldName,pXMLString->c_str(),fieldName,i-1,mgFieldName,main2PostDim(i-1));
		if(!(i==nDims)) {
		  fprintf(outfile,"-1");
		}
		fprintf(outfile,")");
		for(j=i;j<nDims;j++) {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,j);
		}
		fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
	      }
	      fprintf(outfile,"\n");
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"%s += _%s_dx%li;\n",dimName,mgFieldName,main2PostDim(i-1));
	    }
	  }
        
	  if(i==1&simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	    // even if cross-section
	    fprintf(outfile,"	");
	    fprintf(outfile,"}\n");
	  }
	  else
	    if(!(*psamplingLatticeIr==1)) {
	      // if not cross-section
	      for(j=0;j<i;j++) {
		  fprintf(outfile,"	");
	      }
	      fprintf(outfile,"}\n");
	    }
        
	  psamplingLatticeIr++;
	}
    }
  }

  fprintf(outfile,"\n");

  fprintf(outfile,"_%s_x0[_%s_sample_pointer]=%s;\n",mgFieldName,mgFieldName,simulation()->parameters()->propDimName.c_str());
  fprintf(outfile,"\n");

  if(simulation()->parameters()->nPaths==1) {
    if(simulation()->parameters()->usempi) {
      fprintf(outfile,"printf(\"Rank[%%i]  Sampled field (for moment group #%li) at %s	= %%e\\n\",rank,%s);\n",
	      myGroupNumber+1,simulation()->parameters()->propDimName.c_str(),simulation()->parameters()->propDimName.c_str());
    }
    else {
      fprintf(outfile,"printf(\"Sampled field (for moment group #%li) at %s	= %%e\\n\",%s);\n",
	      myGroupNumber+1,simulation()->parameters()->propDimName.c_str(),simulation()->parameters()->propDimName.c_str());
    }
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"_%s_sample_pointer++;\n",mgFieldName);
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");

  // ***************************************
  // ******** process routine  *************
  // ***************************************

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_process() {\n",mgFieldName);
  fprintf(outfile,"\n");
  fprintf(outfile,"double *_to_field;\n");
  if(simulation()->parameters()->nPaths>1) {
    fprintf(outfile,"double *_to_field_sd;\n");
  }

  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    if(complexRawVector) {
      fprintf(outfile,"complex %s;\n",pXMLString->c_str());
    }
    else {
      fprintf(outfile,"double %s;\n",pXMLString->c_str());
    }
  }

  if(myPostCode.length()>0) {
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
      fprintf(outfile,"double %s;\n",pXMLString->c_str());
    }
  }

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      if(complexRawVector) {
        fprintf(outfile,"complex* _temp_vector;\n");
      }
      else {
        fprintf(outfile,"double* _temp_vector;\n");
      }
      fprintf(outfile,"unsigned long _length;\n\n");

      fprintf(outfile,"_length = _%s_size*_%s_fullstep_ncomponents;\n",mgFieldName,mgFieldName);
      if(complexRawVector) {
        fprintf(outfile,"_temp_vector = new complex[_length];\n");
      }
      else {
        fprintf(outfile,"_temp_vector = new double[_length];\n");
      }
    
      if(complexRawVector) {
        fprintf(outfile,"MPI_Reduce(_%s_raw,_temp_vector,2*_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n\n",mgFieldName);
      }
      else {
        fprintf(outfile,"MPI_Reduce(_%s_raw,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n\n",mgFieldName);
      }
    
      fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
      fprintf(outfile,"        _%s_raw[_i0]=_temp_vector[_i0];\n",mgFieldName);
      fprintf(outfile,"\n");
      fprintf(outfile,"delete[] _temp_vector;\n\n");
    }

  fprintf(outfile,"\n");

  /* if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic)
     {
     fprintf(outfile,"if(rank==0) {\n");
     } */

  // assign pointers
  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"if(_half_step) {\n");
    fprintf(outfile,"	_to_field = _%s_halfstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"	_to_field_sd = _%s_halfstep_sd;\n",mgFieldName);
    }
    fprintf(outfile,"	}\n");
    fprintf(outfile,"else {\n");
    fprintf(outfile,"	_to_field = _%s_fullstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"	_to_field_sd = _%s_fullstep_sd;\n",mgFieldName);
    }
    fprintf(outfile,"	}\n");
  }
  else {
    fprintf(outfile,"_to_field = _%s_fullstep;\n",mgFieldName);
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"_to_field_sd = _%s_fullstep_sd;\n",mgFieldName);
    }
  }
  fprintf(outfile,"\n");

  vectors2space(outfile,myPostSpace,myPostProcessingVectorNamesList,"");

  fprintf(outfile,"unsigned long _%s_raw_index_pointer=0;\n",mgFieldName);
  fprintf(outfile,"unsigned long _%s_processed_index_pointer=0;\n",mgFieldName);
  fprintf(outfile,"\n");

  if(postSpace(0)) {
    fprintf(outfile,"double k%s = 0;\n",geometry()->dimension(0)->name.c_str());
  }

  fprintf(outfile,
	  "for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",mgFieldName);
  fprintf(outfile,"\n");

  if(!postSpace(0)) {
    fprintf(outfile,"%s = _%s_x0[_i0];\n",geometry()->dimension(0)->name.c_str(),mgFieldName);
  }
  fprintf(outfile,"\n");

  for(i=1; i<myNDims; i++) {
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	"); 
    }
    if(postSpace(i)) {
      fprintf(outfile,"double k%s = 0;\n",geometry()->dimension(i)->name.c_str());
    }
    else {
      fprintf(outfile,"double %s = _%s_xmin%li;\n",geometry()->dimension(i)->name.c_str(),mgFieldName,i);
    }
    fprintf(outfile,"\n");
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li; _i%li++) {\n",i,i,mgFieldName,i,i);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"\n");

  // kernal of loop struct
  i=0;
  for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin();pXMLString != mySamplingMomentsList.end(); pXMLString++) {
    for(j=0;j<myNDims;j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"%s=_%s_raw[_%s_raw_index_pointer+%li];\n",pXMLString->c_str(),mgFieldName,mgFieldName,i);
    i++;
  }
  fprintf(outfile,"\n");

  char indent[64];
  for(i=0;i<myNDims;i++) {
      indent[i]=0x09;
  }
  indent[myNDims]=0;

  const char* dotre="";
  if(complexRawVector) {
    dotre=".re";
  }

  i=0;
  if(myPostCode.length()==0) {
    for(list<XMLString>::const_iterator pXMLString = mySamplingMomentsList.begin(); pXMLString != mySamplingMomentsList.end(); pXMLString++) {
      if(simulation()->parameters()->nPaths>1) {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] += %s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre);
	fprintf(outfile,"%s_to_field_sd[_%s_processed_index_pointer + %li] += %s%s*%s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre,pXMLString->c_str(),dotre);
      }
      else {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] = %s%s;\n",
		indent,mgFieldName,i,pXMLString->c_str(),dotre);
      }
      i++;
    }
  }
  else {
    fprintf(outfile,"// ********** Post propagation code *************\n");
    fprintf(outfile,"%s\n",myPostCode.c_str());
    fprintf(outfile,"// **********************************************\n");
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end(); pXMLString++) {
      if(simulation()->parameters()->nPaths>1) {
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] += %s;\n",indent,mgFieldName,i,pXMLString->c_str());
	fprintf(outfile,"%s_to_field_sd[_%s_processed_index_pointer + %li] += %s*%s;\n",indent,mgFieldName,i,pXMLString->c_str(),pXMLString->c_str());
      }
      else
	fprintf(outfile,"%s_to_field[_%s_processed_index_pointer + %li] = %s;\n",indent,mgFieldName,i,pXMLString->c_str());
      i++;
    }
  }
  fprintf(outfile,"\n");

  fprintf(outfile,"%s_%s_raw_index_pointer += _%s_raw_ncomponents;\n",indent,mgFieldName,mgFieldName);
  fprintf(outfile,"%s_%s_processed_index_pointer += _%s_fullstep_ncomponents;\n",indent,mgFieldName,mgFieldName);

  for(i=myNDims; i>1; i--) {
    fprintf(outfile,"\n");
    if(postSpace(i-1)) {
      for(unsigned long j=0; j<i; j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"k%s += _%s_dk%li;\n",geometry()->dimension(i-1)->name.c_str(),mgFieldName,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",
	      geometry()->dimension(i-1)->name.c_str(),mgFieldName,i-1,mgFieldName,i-1);	
      for(unsigned long j=0; j<i; j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"	k%s -= _%s_lattice%li*_%s_dk%li;\n",
	      geometry()->dimension(i-1)->name.c_str(),mgFieldName,i-1,mgFieldName,i-1);
    }
    else {
      for(unsigned long j=0; j<i; j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"%s += _%s_dx%li;\n",geometry()->dimension(i-1)->name.c_str(),mgFieldName,i-1);
    }
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"}\n");
  }

  if(postSpace(0)) {
    fprintf(outfile,"	k%s += _%s_dk0;\n",geometry()->dimension(0)->name.c_str(),mgFieldName);
    fprintf(outfile,
	    "	if(k%s>((_%s_lattice0-1)/2 + 0.1)*_%s_dk0)\n",geometry()->dimension(0)->name.c_str(),
	    mgFieldName,mgFieldName);	
    fprintf(outfile,"		k%s -= _%s_lattice0*_%s_dk0;\n",geometry()->dimension(0)->name.c_str(),mgFieldName,mgFieldName);
  }

  fprintf(outfile,"	}\n");

  /*  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic)
      fprintf(outfile,"  }\n");  */

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

  // ***************************************
  // ********  write_out routine  **********
  // ***************************************

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _%s_write_out(\n",mgFieldName);
  fprintf(outfile,"	FILE *_outfile) {\n");
  fprintf(outfile,"\n");
  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"double _max_step_error=0;\n");
  }
  fprintf(outfile,"\n");

  fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_size*_%s_fullstep_ncomponents; _i0++) {\n",mgFieldName,mgFieldName);

  if(simulation()->parameters()->nPaths>1) {
    fprintf(outfile,"	_%s_fullstep[_i0] *= 1/(double)_n_paths;\n",mgFieldName);
    fprintf(outfile,"	_%s_fullstep_sd[_i0] *= 1/(double)_n_paths;\n",mgFieldName);
    fprintf(outfile,"	_%s_fullstep_sd[_i0] -= _%s_fullstep[_i0]*_%s_fullstep[_i0];\n",mgFieldName,mgFieldName,mgFieldName);
    fprintf(outfile,"	if(_%s_fullstep_sd[_i0]>0)\n",mgFieldName);
    fprintf(outfile,"		_%s_fullstep_sd[_i0]=sqrt(_%s_fullstep_sd[_i0]/_n_paths);\n",mgFieldName,mgFieldName);
    fprintf(outfile,"	else\n");
    fprintf(outfile,"		_%s_fullstep_sd[_i0]=0;\n",mgFieldName);
    fprintf(outfile,"\n");
    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"	_%s_halfstep[_i0] *= 1/(double)_n_paths;\n",mgFieldName);
      fprintf(outfile,"	_%s_halfstep_sd[_i0] *= 1/(double)_n_paths;\n",mgFieldName);
      fprintf(outfile,"	_%s_halfstep_sd[_i0] -= _%s_halfstep[_i0]*_%s_halfstep[_i0];\n",mgFieldName,mgFieldName,mgFieldName);
      fprintf(outfile,"	if(_%s_halfstep_sd[_i0]>0)\n",mgFieldName);
      fprintf(outfile,"		_%s_halfstep_sd[_i0]=sqrt(_%s_halfstep_sd[_i0]/_n_paths);\n",mgFieldName,mgFieldName);
      fprintf(outfile,"	else\n");
      fprintf(outfile,"		_%s_halfstep_sd[_i0]=0;\n",mgFieldName);
      fprintf(outfile,"\n");
    }
  }

  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"	double _step_error=_%s_fullstep[_i0] - _%s_halfstep[_i0];\n",mgFieldName,mgFieldName);
    fprintf(outfile,"	if(fabs(_step_error)>_max_step_error)\n");
    fprintf(outfile,"		_max_step_error=fabs(_step_error);\n");
  }

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


  if(simulation()->parameters()->errorCheck) {
    if(simulation()->parameters()->nPaths>1) {
      fprintf(outfile,"printf(\"maximum step error in moment group %li means was %%e\\n\",_max_step_error);\n",myGroupNumber+1);
    }
    else {
      fprintf(outfile,"printf(\"maximum step error in moment group %li was %%e\\n\",_max_step_error);\n",myGroupNumber+1);
    }
    fprintf(outfile,"\n");
  }

  // this is where the moment groups are written
  fprintf(outfile,"fprintf(_outfile,\"\\n\");\n");
  fprintf(outfile,"fprintf(_outfile,\"<XSIL Name=\\\"moment_group_%li\\\">\\n\");\n",myGroupNumber+1);

  unsigned long nVariables = myNDims-1;

  if(nSamples>1) {
    nVariables++;
  }

  fprintf(outfile,"fprintf(_outfile,\"	<Param Name=\\\"n_independent\\\">%li</Param>\\n\");\n",nVariables);
  fprintf(outfile,"fprintf(_outfile,\"	<Array Name=\\\"variables\\\" Type=\\\"Text\\\">\\n\");\n");

  nVariables += myPostMomentsList.size();
  if(simulation()->parameters()->errorCheck) {
    nVariables += myPostMomentsList.size();
  }
  if(simulation()->parameters()->nPaths>1) {
    nVariables += myPostMomentsList.size();
  }

  fprintf(outfile,"fprintf(_outfile,\"		<Dim>%li</Dim>\\n\");\n",nVariables);
  fprintf(outfile,"fprintf(_outfile,\"		<Stream><Metalink Format=\\\"Text\\\" Delimiter=\\\" \\\\n\\\"/>\\n\");\n");
  fprintf(outfile,"fprintf(_outfile,\"");
  if(nSamples>1) {
    if(postSpace(0)) {
	fprintf(outfile,"k");
    }
    fprintf(outfile,"%s ",simulation()->parameters()->propDimName.c_str());
  }
  for(i=1;i<myNDims;i++) {
    if(postSpace(i)) {
	fprintf(outfile,"k");
    }
    fprintf(outfile,"%s ",geometry()->dimension(i)->name.c_str());
  }
  if(simulation()->parameters()->nPaths>1) {
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end();pXMLString++) {
      fprintf(outfile,"mean_%s ",pXMLString->c_str());
    }
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end();pXMLString++) {
      fprintf(outfile,"sd_%s ",pXMLString->c_str());
    }
  }
  else {
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end();pXMLString++) {
      fprintf(outfile,"%s ",pXMLString->c_str());
    }
  }

  if(simulation()->parameters()->errorCheck) {
    for(list<XMLString>::const_iterator pXMLString = myPostMomentsList.begin(); pXMLString != myPostMomentsList.end();pXMLString++) {
      fprintf(outfile,"error_%s ",pXMLString->c_str());
    }
  }
  fprintf(outfile,"\\n\");\n");
  fprintf(outfile,"fprintf(_outfile,\"			</Stream>\\n\");\n");
  fprintf(outfile,"fprintf(_outfile,\"		</Array>\\n\");\n");
  fprintf(outfile,"fprintf(_outfile,\"	<Array Name=\\\"data\\\" Type=\\\"double\\\">\\n\");\n");
  if(nSamples>1) {
    fprintf(outfile,"fprintf(_outfile,\"		<Dim>%li</Dim>\\n\");\n",nSamples);
  }
  for(i=1;i<myNDims;i++) {
    fprintf(outfile,"fprintf(_outfile,\"		<Dim>%li</Dim>\\n\");\n",geometry()->dimension(i)->lattice);
  }
  fprintf(outfile,"fprintf(_outfile,\"		<Dim>%li</Dim>\\n\");\n",nVariables);
  // Should the output to the .xsil file be binary or ascii?
  if (simulation()->parameters()->binaryOutput) {
    // Are we using single or double precision output?
    string precisionString, precision;  // the "precision" variable is for the xsil output, and for matlab
    string castString;
    if (simulation()->parameters()->useDouble) {
      precisionString = "double";
      precision = "double";
      castString = "";
    }
    else {
      precisionString = "float";
      precision = "single";
      castString = "(float)";
    }

    fprintf(outfile,
	    "\nstring encodingStr;\n"
	    "if (CPU_IS_BIG_ENDIAN) {\n"
	    "    // big endian\n"
	    "    encodingStr = \"BigEndian\";\n"
	    "}\n"
	    "else if (CPU_IS_LITTLE_ENDIAN) {\n"
	    "    // little endian\n"
	    "    encodingStr = \"LittleEndian\";\n"
	    "}\n"
	    "else {\n"
	    "    // dunno what the byte order is\n"
	    "    cout << \"I don't know what kind of byte ordering you're using\\n\";\n"
	    "    cout << \"Using \\\"Binary\\\" as encoding\\n\";\n"
	    "    encodingStr = \"Binary\";\n"
	    "}\n"
	    "\nstring ulongTypeStr;\n"
	    "if (SIZEOF_UNSIGNED_LONG==4) {\n"
	    "    // 32 bit\n"
	    "    ulongTypeStr = \"uint32\";\n"
	    "}\n"
	    "else if (SIZEOF_UNSIGNED_LONG==8) {\n"
	    "    // 64 bit\n"
	    "    ulongTypeStr = \"uint64\";\n"
	    "}\n"
	    "else {\n"
	    "    // dunno what the default size of the unsigned long is\n"
	    "    ulongTypeStr = \"ulong\";\n"
	    "}\n"
	    "fprintf(_outfile,\"	<Stream><Metalink Format=\\\"Binary\\\" UnsignedLong=\\\"%%s\\\" \",ulongTypeStr.c_str());\n"
	    "fprintf(_outfile,\"precision=\\\"%s\\\" Type=\\\"Remote\\\" Encoding=\\\"%%s\\\"/>\",encodingStr.c_str());\n",precision.c_str());

    // get the xsil filename, rip off the extension if it equals '.xsil'
    // otherwise, leave it alone
    XMLString xsilFilename = simulation()->output()->getOutputFileName();
    XMLString xsilExtension;
    xsilExtension += ".xsil";
    XMLString testExtension;
    xsilFilename.subString(testExtension,xsilFilename.length()-5,xsilFilename.length());
    if (testExtension == xsilExtension) {
      xsilFilename.deleteData(xsilFilename.length()-5,5);
    }

    fprintf(outfile, "fprintf(_outfile,\"\\n%s_mg%li%s\\n\");\n",xsilFilename.c_str(),myGroupNumber+1,".dat");
    fprintf(outfile, "FILE *fpBinary;\n");
    fprintf(outfile, "if ((fpBinary = fopen(\"%s_mg%li%s\",\"wb\")) == NULL) {\n",xsilFilename.c_str(),myGroupNumber+1,".dat");
    fprintf(outfile, "	printf(\"Unable to open output file %s_mg%li%s\\n\");\n",xsilFilename.c_str(),myGroupNumber+1,".dat");
    fprintf(outfile, "	printf(\"Chucking a spack....\\n\");\n");
    fprintf(outfile, "  exit(255);}\n");
    fprintf(outfile, "vector<%s> ",precisionString.c_str());
    for (i=0;i<nVariables-1;i++) {
      fprintf(outfile, "_outVar%li, ",i);
    }
    fprintf(outfile, "_outVar%li;\n",nVariables-1);
    if (nSamples>1) {
      for (i=0;i<myNDims;i++) {
	fprintf(outfile, "_outVar%li.resize(_%s_lattice%li,0);\n",i,mgFieldName,i);
      }
    }
    else {
      for (i=0;i<myNDims-1;i++) {
	fprintf(outfile, "_outVar%li.resize(_%s_lattice%li,0);\n",i,mgFieldName,i+1);
      }
    }

    if(postSpace(0)) {
      fprintf(outfile,"for(long _j0=-(_%s_lattice0)/2;_j0<(_%s_lattice0 + 1)/2;_j0++) {\n",mgFieldName,mgFieldName);
      fprintf(outfile,"\n");
      fprintf(outfile,"	long _i0 = _j0;\n");
      fprintf(outfile,"	if(_i0<0)\n");
      fprintf(outfile,"		_i0 += _%s_lattice0;\n",mgFieldName);
      fprintf(outfile,"	double _k0 = _j0*_%s_dk0;\n",mgFieldName);
      fprintf(outfile,"\n");
    }
    else {
      fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",mgFieldName);
      fprintf(outfile,"\n");
    }

    for(i=1;i<myNDims;i++) {
      if(postSpace(i)) {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"long signed int _loopIter%li = -1;\n",i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"for(long _j%li=-(_%s_lattice%li)/2;_j%li<(_%s_lattice%li + 1)/2;_j%li++) {\n",i,mgFieldName,i,i,mgFieldName,i,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"       _loopIter%li++;\n",i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	long _i%li = _j%li;\n",i,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	if(_i%li<0)\n",i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"		_i%li += _%s_lattice%li;\n",i,mgFieldName,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	double _k%li = _j%li*_%s_dk%li;\n",i,i,mgFieldName,i);
	fprintf(outfile,"\n");
      }
      else {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"long signed int _loopIter%li = -1;\n",i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"double _x%li=_%s_xmin%li;\n",i,mgFieldName,i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,i,i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"       _loopIter%li++;\n",i);
	fprintf(outfile,"\n");
      }
    }

    for(i=0;i<myNDims;i++) {
	fprintf(outfile,"	");  // code indentation
    }
    fprintf(outfile,"unsigned long _%s_fullstep_index_pointer=0;\n",mgFieldName);
    for(i=0;i<myNDims;i++) {
      for(j=0;j<myNDims;j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_fullstep_index_pointer += _i%li",mgFieldName,i);
      for(j=i+1;j<myNDims;j++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,j);
      }
      fprintf(outfile,"*_%s_fullstep_ncomponents;\n",mgFieldName);
    }
    fprintf(outfile,"\n");

    unsigned long int dimsCount=0;
    if(nSamples>1) {
      fprintf(outfile,"_outVar%li[_i0] = ",dimsCount);
      dimsCount++;
      if(postSpace(0)) {
	fprintf(outfile,"_k0");
      }
      else {
	fprintf(outfile,"%s_%s_x0[_i0]",castString.c_str(),mgFieldName);
      }
      fprintf(outfile, ";\n");
    }

    for(i=1;i<myNDims;i++){
      // fprintf(outfile, "_outVar%li[_i%li] = ",dimsCount,i);
      fprintf(outfile, "_outVar%li[_loopIter%li] = ",dimsCount,i);
      dimsCount++;
      if(postSpace(i)){
	fprintf(outfile,"%s_k%li",castString.c_str(),i);
      }
      else{
	fprintf(outfile,"%s_x%li",castString.c_str(),i);}
      fprintf(outfile, ";\n");
    }

    if(simulation()->parameters()->errorCheck) {
      for(i=0;i<myPostMomentsList.size();i++){
	fprintf(outfile, "_outVar%li.push_back(",dimsCount);
	dimsCount++;
	fprintf(outfile,"%s_%s_halfstep[_%s_fullstep_index_pointer + %li]",castString.c_str(),mgFieldName,mgFieldName,i);
	fprintf(outfile,");\n");
      }
      if(simulation()->parameters()->nPaths>1) {
	for(i=0;i<myPostMomentsList.size();i++){
	  fprintf(outfile, "_outVar%li.push_back(",dimsCount);
	  dimsCount++;
	  fprintf(outfile,"%s_%s_halfstep_sd[_%s_fullstep_index_pointer + %li]",castString.c_str(),mgFieldName,mgFieldName,i);
	  fprintf(outfile,");\n");
	}	
      }
      for(i=0;i<myPostMomentsList.size();i++){
	fprintf(outfile,"double _temp%li = ",i);
	fprintf(outfile,"_%s_fullstep[_%s_fullstep_index_pointer + %li]-_%s_halfstep[_%s_fullstep_index_pointer + %li];\n",
		mgFieldName,mgFieldName,i,mgFieldName,mgFieldName,i);
	fprintf(outfile,"_outVar%li.push_back(_temp%li",dimsCount,i);
	dimsCount++;
	fprintf(outfile,");\n");
      }
    }
    else {
      for(i=0;i<myPostMomentsList.size();i++) {
	fprintf(outfile, "_outVar%li.push_back(",dimsCount);
	dimsCount++;
	fprintf(outfile,"%s_%s_fullstep[_%s_fullstep_index_pointer + %li]",castString.c_str(),mgFieldName,mgFieldName,i);
	fprintf(outfile,");\n");
      }
      if(simulation()->parameters()->nPaths>1) {
	for(i=0;i<myPostMomentsList.size();i++){
	  fprintf(outfile,"_outVar%li.push_back(",dimsCount);
	  dimsCount++;
	  fprintf(outfile,"%s_%s_fullstep_sd[_%s_fullstep_index_pointer + %li]",castString.c_str(),mgFieldName,mgFieldName,i);
	  fprintf(outfile,");\n");
	}
      }
    }

    for(i=myNDims;i>1;i--) {
      fprintf(outfile,"\n");
      if(!postSpace(i-1)) {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"_x%li += _%s_dx%li*_%s_lattice%li/_%s_lattice%li;\n",
		i-1,simulation()->field()->name()->c_str(),post2MainDim(i-1),
		simulation()->field()->name()->c_str(),post2MainDim(i-1),mgFieldName,i-1);
      }
      for(j=0;j<i;j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"}\n");
    }

    fprintf(outfile,"	}\n");

    fprintf(outfile, "unsigned long int ii, ");
    for (i=0;i<nVariables-1;i++) {
      fprintf(outfile, "_outVarLen%li, ",i);
    }
    fprintf(outfile, "_outVarLen%li;\n",nVariables-1);
    for (i=0; i<nVariables; i++) {
      fprintf(outfile, "_outVarLen%li = _outVar%li.size();\n",i,i);
    }
    for (i=0; i<nVariables; i++) {
      fprintf(outfile, "fwrite(&_outVarLen%li,sizeof(unsigned long int),1,fpBinary);\n",i);
      fprintf(outfile, "for (ii=0; ii<_outVarLen%li; ii++) {\n",i);
      fprintf(outfile, "fwrite(&_outVar%li[ii],sizeof(%s),1,fpBinary);\n",i,precisionString.c_str());
      fprintf(outfile, "}\n\n");
    }

    // close the binary file
    fprintf(outfile, "fclose(fpBinary);\n");
    fprintf(outfile, "fprintf(_outfile,\"                </Stream>\\n\");\n");
  }
  else {  // use this bit if not using binary output
    fprintf(outfile,"fprintf(_outfile,\"		<Stream><Metalink Format=\\\"Text\\\" Delimiter=\\\" \\\\n\\\"/>\\n\");\n");
    fprintf(outfile,"\n");

    if(postSpace(0)) {
      fprintf(outfile,"for(long _j0=-(_%s_lattice0)/2;_j0<(_%s_lattice0 + 1)/2;_j0++) {\n",mgFieldName,mgFieldName);
      fprintf(outfile,"\n");
      fprintf(outfile,"	long _i0 = _j0;\n");
      fprintf(outfile,"	if(_i0<0)\n");
      fprintf(outfile,"		_i0 += _%s_lattice0;\n",mgFieldName);
      fprintf(outfile,"	double _k0 = _j0*_%s_dk0;\n",mgFieldName);
      fprintf(outfile,"\n");
    }
    else {
      fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_lattice0; _i0++) {\n",mgFieldName);
      fprintf(outfile,"\n");
    }

    for(i=1;i<myNDims;i++) {
      if(postSpace(i)) {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	"); 
	}
	fprintf(outfile,"for(long _j%li=-(_%s_lattice%li)/2;_j%li<(_%s_lattice%li + 1)/2;_j%li++) {\n",i,mgFieldName,i,i,mgFieldName,i,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	long _i%li = _j%li;\n",i,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	if(_i%li<0)\n",i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"		_i%li += _%s_lattice%li;\n",i,mgFieldName,i);
	fprintf(outfile,"\n");
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"	double _k%li = _j%li*_%s_dk%li;\n",i,i,mgFieldName,i);
	fprintf(outfile,"\n");
      }
      else {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"double _x%li=_%s_xmin%li;\n",i,mgFieldName,i);
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"for(unsigned long _i%li=0;_i%li<_%s_lattice%li;_i%li++) {\n",i,i,mgFieldName,i,i);
	fprintf(outfile,"\n");
      }
    }

    for(i=0;i<myNDims;i++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"unsigned long _%s_fullstep_index_pointer=0;\n",mgFieldName);
    for(i=0;i<myNDims;i++) {
      for(j=0;j<myNDims;j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"_%s_fullstep_index_pointer += _i%li",mgFieldName,i);
      for(j=i+1;j<myNDims;j++) {
	fprintf(outfile,"*_%s_lattice%li",mgFieldName,j);
      }
      fprintf(outfile,"*_%s_fullstep_ncomponents;\n",mgFieldName);
    }
    fprintf(outfile,"\n");

    for(i=0;i<myNDims;i++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"fprintf(_outfile,\"%%e");
    for(i=0;i<nVariables-1;i++) {
      fprintf(outfile," %%e");
    }
    fprintf(outfile,"\\n\"");

    if(nSamples>1) {
      if(postSpace(0)) {
	fprintf(outfile,",_k0");
      }
      else {
	fprintf(outfile,",_%s_x0[_i0]",mgFieldName);
      }
    }

    for(i=1;i<myNDims;i++) {
      if(postSpace(i)) {
	fprintf(outfile,",_k%li",i);
      }
      else {
	fprintf(outfile,",_x%li",i);
      }
    }

    if(simulation()->parameters()->errorCheck) {
      for(i=0;i<myPostMomentsList.size();i++) {
	fprintf(outfile,",_%s_halfstep[_%s_fullstep_index_pointer + %li]",mgFieldName,mgFieldName,i);
      }
      if(simulation()->parameters()->nPaths>1) {
	for(i=0;i<myPostMomentsList.size();i++) {
	  fprintf(outfile,",_%s_halfstep_sd[_%s_fullstep_index_pointer + %li]",mgFieldName,mgFieldName,i);
	}
      }
      for(i=0;i<myPostMomentsList.size();i++) {
	fprintf(outfile,",_%s_fullstep[_%s_fullstep_index_pointer + %li]-_%s_halfstep[_%s_fullstep_index_pointer + %li]",
		mgFieldName,mgFieldName,i,mgFieldName,mgFieldName,i);
      }
    }
    else {
      for(i=0;i<myPostMomentsList.size();i++) {
	fprintf(outfile,",_%s_fullstep[_%s_fullstep_index_pointer + %li]",mgFieldName,mgFieldName,i);
      }
      if(simulation()->parameters()->nPaths>1) {
	for(i=0;i<myPostMomentsList.size();i++) {
	  fprintf(outfile,",_%s_fullstep_sd[_%s_fullstep_index_pointer + %li]",mgFieldName,mgFieldName,i);
	}
      }
    }

    fprintf(outfile,");\n");
    for(i=myNDims;i>1;i--) {
      fprintf(outfile,"\n");
      if(!postSpace(i-1)) {
	for(j=0;j<i;j++) {
	    fprintf(outfile,"	");
	}
	fprintf(outfile,"_x%li += _%s_dx%li*_%s_lattice%li/_%s_lattice%li;\n",
		i-1,simulation()->field()->name()->c_str(),post2MainDim(i-1),
		simulation()->field()->name()->c_str(),post2MainDim(i-1),mgFieldName,i-1);
      }
      for(j=0;j<i;j++) {
	  fprintf(outfile,"	");
      }
      fprintf(outfile,"}\n");
    }

    fprintf(outfile,"	}\n");

    fprintf(outfile,"fprintf(_outfile,\"			</Stream>\\n\");\n");
  }  // end of stuff distinguishing between binary and ascii output
  fprintf(outfile,
	  "fprintf(_outfile,\"		</Array>\\n\");\n"
	  "fprintf(_outfile,\"	</XSIL>\\n\");\n"
	  "}\n"
	  "\n");
};

// ******************************************************************************
bool xmdsMomentGroup::samplingSpace(
				    const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::samplingSpace\n");
  }

  if(index>=simulation()->field()->geometry()->nDims()) {
    throw xmdsException("Internal range error in xmdsMomentGroup::samplingSpace()");
  }

  return (mySamplingSpace >> index)&1;
};

// ******************************************************************************
bool xmdsMomentGroup::postSpace(
				const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::postSpace\n");
  }

  if(index>=geometry()->nDims()) {
    throw xmdsException("Internal range error in xmdsMomentGroup::postSpace()");
  }

  return (myPostSpace >> index)&1;
};

// ******************************************************************************
unsigned long xmdsMomentGroup::main2PostDim(
					    const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::main2PostDim\n");
  }

  if(index>=myMain2PostDimList.size()) {
    throw xmdsException("Internal range error in xmdsMomentGroup::main2PostDim()");
  }

  list<unsigned long>::const_iterator pULong = myMain2PostDimList.begin();
  for(unsigned long i=0; i<index; i++) {
    pULong++;
  }
  return *pULong;
};

// ******************************************************************************
unsigned long xmdsMomentGroup::post2MainDim(
					    const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsMomentGroup::post2MainDim\n");
  }

  if(index>=myPost2MainDimList.size()) {
    throw xmdsException("Internal range error in xmdsMomentGroup::post2MainDim()");
  }

  list<unsigned long>::const_iterator pULong = myPost2MainDimList.begin();
  for(unsigned long i=0; i<index; i++) {
    pULong++;
  }
  return *pULong;
};