/*
 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: xmdsintegratesiex.cc,v 1.19 2004/10/06 07:52:38 joehope Exp $
*/

/*! @file xmdsintegratesiex.cc
  @brief Integrate element parsing classes and methods; semi-implicit method in the explicit picture

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateSIEX public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateSIEXs=0;  //!< The number of xmds integrate SIEX objects

// ******************************************************************************
xmdsIntegrateSIEX::xmdsIntegrateSIEX(
				     const xmdsSimulation *const yourSimulation,
				     const bool& yourVerboseMode) :
  xmdsIntegrate(yourSimulation,yourVerboseMode),
  xmdsIntegrateEX(yourSimulation,yourVerboseMode),
  xmdsIntegrateSI(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsIntegrateSIEXs++;
    printf("xmdsIntegrateSIEX::xmdsIntegrateSIEX\n");
    printf("nxmdsIntegrateSIEXs=%li\n",nxmdsIntegrateSIEXs);
  }
};

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

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

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

  xmdsIntegrate::processElement(yourElement);
  xmdsIntegrateSI::processElement(yourElement);
  xmdsIntegrateEX::processElement(yourElement);
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateSIEX private
// ******************************************************************************
// ******************************************************************************

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

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// segment %li (SIEX) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");

  xmdsIntegrate::writePrototypes(outfile);
  xmdsIntegrateEX::writePrototypes(outfile);

  fprintf(outfile,"void _segment%li(unsigned long cycle);\n",segmentNumber);
  fprintf(outfile,"\n");

  if(crossVectorNamesList()->size() > 0) {
    fprintf(outfile,"void _segment%li_calculate_cross_field();\n",segmentNumber);
    fprintf(outfile,"\n");
  }
};

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

  fprintf(outfile,"// ********************************************************\n");
  fprintf(outfile,"// segment %li (SIEX) routines\n",segmentNumber);
  fprintf(outfile,"\n");

  xmdsIntegrate::writeRoutines(outfile);
  xmdsIntegrateEX::writeRoutines(outfile);

  writeMainIntegrateRoutine(outfile);

  if(crossVectorNamesList()->size() > 0) {
    writeCalculateCrossFieldRoutine(outfile);
  }
};

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

  const char *const fieldName = simulation()->field()->name()->c_str();

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateSIEX::writeMainIntegrateRoutine: cannot find 'main' vector");
  }

  const char* typeName="";
  if(mainVector->vectorType()==COMPLEX) {
    typeName="complex";
  }
  else if(mainVector->vectorType()==DOUBLE) {
    typeName="double";
  }

  fprintf(outfile,"/* ******************************************** */\n");
  fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s *aifield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
  fprintf(outfile,"\n");
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    fprintf(outfile,"const double _var = 1");
    for(unsigned long i=0;i<simulation()->field()->geometry()->nDims();i++) {
      fprintf(outfile,"/_%s_dx%li",fieldName,i);
    }
    fprintf(outfile,";\n");
    fprintf(outfile,"double *_noise_vector = new double[_%s_size*_n_noises];\n",fieldName);
    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"double *_noise_vector2 = new double[_%s_size*_n_noises];\n",fieldName);
    }
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"for(unsigned long _i0=0;_i0<%li;_i0++) {\n",lattice());
  fprintf(outfile,"\n");

  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"	if(_half_step) {\n");
    fprintf(outfile,"\n");
    fprintf(outfile,"		const double _step = 0.5*%s/(double)%li;\n",interval()->c_str(),lattice());
    fprintf(outfile,"\n");
    if((simulation()->parameters()->stochastic)&&(!noNoises())) {
      fprintf(outfile,"		_make_noises(_gen1,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"\n");
    }		
    writeSingleStepCode(outfile,FIRST_HALFSTEP);
    fprintf(outfile,"\n");
    if((simulation()->parameters()->stochastic)&&(!noNoises())) {
      fprintf(outfile,"		_make_noises(_gen2,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"\n");
    }		
    writeSingleStepCode(outfile,SECOND_HALFSTEP);
    fprintf(outfile,"		}\n");
    fprintf(outfile,"	else {\n");
    fprintf(outfile,"\n");
    fprintf(outfile,"		const double _step = %s/(double)%li;\n",interval()->c_str(),lattice());
    fprintf(outfile,"\n");
    if((simulation()->parameters()->stochastic)&&(!noNoises())) {
      fprintf(outfile,"		_make_noises(_gen1,_var/_step/2,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"		_make_noises(_gen2,_var/_step/2,_noise_vector2,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"		for(unsigned long _i1=0;_i1<_%s_size*_n_noises;_i1++)\n",fieldName);
      fprintf(outfile,"			_noise_vector[_i1] += _noise_vector2[_i1];\n");
      fprintf(outfile,"\n");
    }		
    writeSingleStepCode(outfile,FULLSTEP);
    fprintf(outfile,"		}\n");
  }
  else {
    fprintf(outfile,"\n");
    fprintf(outfile,"	const double _step = %s/(double)%li;\n",interval()->c_str(),lattice());
    fprintf(outfile,"\n");
    if((simulation()->parameters()->stochastic)&&(!noNoises())) {
      fprintf(outfile,"	_make_noises(_gen,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"\n");
    }		
    writeSingleStepCode(outfile,FULLSTEP);
  }

  for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) {
    if(samples(i)!=0){
      fprintf(outfile,"\n");
      fprintf(outfile,"	if(%li*((_i0+1)/%li)==(_i0+1))\n",lattice()/samples(i),lattice()/samples(i));
      fprintf(outfile,"		_mg%li_sample();\n",i);
      }
  }
  fprintf(outfile,"	}\n");

  fprintf(outfile,"    delete[] aifield_%s_main;",fieldName);
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    fprintf(outfile,"    delete[] _noise_vector;\n");
    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"    delete[] _noise_vector2;\n");
    }
    fprintf(outfile,"\n");
  }

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

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

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_cross_field() {",segmentNumber);
  fprintf(outfile,"\n");
  if(crossDimNumber()+1<nDims) {
    fprintf(outfile,"const unsigned long _%s_cross_size = _%s_lattice%li",fieldName,fieldName,crossDimNumber()+1);
    for(unsigned long i=crossDimNumber()+2; i<nDims; i++) {
      fprintf(outfile,"*_%s_lattice%li",fieldName,i);
    }
    fprintf(outfile,";\n");
  }
  else
    fprintf(outfile,"const unsigned long _%s_cross_size = 1;\n",fieldName);
  fprintf(outfile,"\n");

  // need to create a vector list for ONLY main vectors

  list<XMLString> myMainVectorNamesList;

  for(list<XMLString>::const_iterator pXMLString = vectorNamesList()->begin(); pXMLString != vectorNamesList()->end(); pXMLString++) {
    list<XMLString>::const_iterator pXMLString2 = crossVectorNamesList()->begin();
    while((pXMLString2 != crossVectorNamesList()->end()) && (*pXMLString2 != *pXMLString)) {
      pXMLString2++;
    }
    if(*pXMLString2 != *pXMLString) {
      myMainVectorNamesList.push_back(*pXMLString);
    }
  }

  const char* typeName;
  list<const xmdsVector*> mainVectorList;

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

    const xmdsVector* mainVector;

    if(!simulation()->field()->getVector(*pXMLString,mainVector)) {
      throw xmdsException("Internal error in xmdsIntegrateSIEX::writeCalculateCrossFieldRoutine: cannot find main vector");
    }

    mainVectorList.push_back(mainVector);

    if(mainVector->vectorType()==DOUBLE) {
      typeName="double";
    }
    else {
      typeName="complex";
    }

    fprintf(outfile,"%s *_%s_%s_old = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
	    typeName,fieldName,pXMLString->c_str(),typeName,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  list<const xmdsVector*> crossVectorList;

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

    const xmdsVector* crossVector;

    if(!simulation()->field()->getVector(*pXMLString,crossVector)) {
      throw xmdsException("Internal error in xmdsIntegrateSIEX::writeCalculateCrossFieldRoutine: cannot find cross vector");
    }

    crossVectorList.push_back(crossVector);

    if(crossVector->vectorType()==DOUBLE) {
      typeName="double";
    }
    else {
      typeName="complex";
    }

    fprintf(outfile,"%s *_%s_%s_I = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
	    typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,fieldName,crossVector->name()->c_str());
    fprintf(outfile,"%s *_%s_%s_d = new %s[_%s_cross_size*_%s_%s_ncomponents];\n",
	    typeName,fieldName,crossVector->name()->c_str(),typeName,fieldName,fieldName,crossVector->name()->c_str());
    fprintf(outfile,"\n");

    for(unsigned long i=0;i<crossVector->nComponents();i++) {
      fprintf(outfile,"%s d%s_d%s;\n",
	      typeName,crossVector->componentName(i)->c_str(),
	      simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str());
    }
    fprintf(outfile,"\n");
  }

  // add cross vectors to total vectors to use

  list<XMLString> myTotalVectorsList = myMainVectorNamesList;

  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    myTotalVectorsList.push_back(*pXMLString);
  }


  simulation()->field()->vectors2space(outfile,0,myTotalVectorsList,"");

  // open outer loops

  for(unsigned long i=0; i<crossDimNumber(); i++) {

    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,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,fieldName,i,i);
    fprintf(outfile,"\n");
  }

  for(unsigned long j=0; j<crossDimNumber(); j++) {
      fprintf(outfile,"	");
  }
  fprintf(outfile,"double %s = _%s_xmin%li;\n",
	  simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
  fprintf(outfile,"\n");

  for(unsigned long j=0; j<crossDimNumber(); j++) {
      fprintf(outfile,"	");
  }
  fprintf(outfile,"for(unsigned long _i%li=0; _i%li<_%s_lattice%li-1; _i%li++) {\n",
	  crossDimNumber(),crossDimNumber(),fieldName,crossDimNumber(),crossDimNumber());
  fprintf(outfile,"\n");

  char indent[64];
  for(unsigned long i=0;i<crossDimNumber()+1;i++) {
    indent[i]=0x09;
  }
  indent[crossDimNumber()+1]=0;

  char indent2[64];
  for(unsigned long i=0;i<nDims;i++) {
    indent2[i]=0x09;
  }
  indent2[nDims]=0;

  for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
    fprintf(outfile,"%sunsigned long _%s_%s_index_pointer_begin=0;\n",indent,fieldName,pXMLString->c_str());
    for(unsigned long i=0; i<crossDimNumber()+1; i++) {
      fprintf(outfile,"%s_%s_%s_index_pointer_begin += _i%li",indent,fieldName,pXMLString->c_str(),i);
      for(unsigned long j=i+1;j<nDims;j++) {
	fprintf(outfile,"*_%s_lattice%li",fieldName,j);
      }
      fprintf(outfile,"*_%s_%s_ncomponents;\n",fieldName,pXMLString->c_str());
    }
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// copy cross vectors into I vectors\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s	_%s_%s_I[_j] = _%s_%s[_%s_%s_index_pointer_begin + _j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// store main vectors into old\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
    fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s	_%s_%s_old[_j] = _active_%s_%s[_%s_%s_index_pointer_begin + _j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// create midpoint main vectors\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
    fprintf(outfile,"%sfor(unsigned long _j=_%s_%s_index_pointer_begin; _j < _%s_%s_index_pointer_begin + _%s_cross_size*_%s_%s_ncomponents; _j++)\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s	_active_%s_%s[_j] = (_active_%s_%s[_j]+_active_%s_%s[_%s_cross_size*_%s_%s_ncomponents + _j])/2;\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),
	    fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"%s// move cross dim to lattice midpoint\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s%s += _%s_dx%li/2;\n",
	  indent,simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
  fprintf(outfile,"\n");

  fprintf(outfile,"%sfor(unsigned long _j=0; _j<%li; _j++) {\n",indent,nIterations());
  fprintf(outfile,"\n");

  for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
    fprintf(outfile,"%s	unsigned long _%s_%s_index_pointer=_%s_%s_index_pointer_begin;\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
  }
  fprintf(outfile,"\n");

  fprintf(outfile,"%s	// calculate delta a\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    fprintf(outfile,"%s	unsigned long _%s_%s_index_pointer_local=0;\n",indent,fieldName,pXMLString->c_str());
  }
  fprintf(outfile,"\n");

  // open inner loops
  for(unsigned long i=crossDimNumber()+1; i<nDims; i++) {

    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"	double %s = _%s_xmin%li;\n",simulation()->field()->geometry()->dimension(i)->name.c_str(),fieldName,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,fieldName,i,i);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"// *********** cross_propagation code ***********\n");
  fprintf(outfile,"%s\n",crossPropagationCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  for(list<const xmdsVector*>::const_iterator pxmdsVector = crossVectorList.begin(); pxmdsVector != crossVectorList.end(); pxmdsVector++) {
    for(unsigned long i=0;i<(*pxmdsVector)->nComponents();i++) {
      fprintf(outfile,"%s	_%s_%s_d[_%s_%s_index_pointer_local + %li] = d%s_d%s*_%s_dx%li;\n",
	      indent2,fieldName,(*pxmdsVector)->name()->c_str(),fieldName,(*pxmdsVector)->name()->c_str(),i,
	      (*pxmdsVector)->componentName(i)->c_str(),
	      simulation()->field()->geometry()->dimension(crossDimNumber())->name.c_str(),fieldName,crossDimNumber());
    }
    fprintf(outfile,"\n");
  }

  //close inner loops
  if(crossDimNumber()+1<nDims) {
    for(list<XMLString>::const_iterator pXMLString = myTotalVectorsList.begin(); pXMLString != myTotalVectorsList.end(); pXMLString++) {
      fprintf(outfile,"%s	_%s_%s_index_pointer += _%s_%s_ncomponents;\n",
	      indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    }
    fprintf(outfile,"\n");
    for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
      fprintf(outfile,"%s	_%s_%s_index_pointer_local += _%s_%s_ncomponents;\n",
	      indent2,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    }
    fprintf(outfile,"\n");
  }

  for(unsigned long i=nDims; i>crossDimNumber()+1; i--) {
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"	%s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
    fprintf(outfile,"\n");
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"	}\n");
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s	// create next cross vectors\n",indent);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s	if(_j<%li) {\n",indent,nIterations()-1);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    fprintf(outfile,"%s		for(unsigned long _k=0; _k<_%s_cross_size*_%s_%s_ncomponents; _k++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s			_%s_%s[_%s_%s_index_pointer_begin + _k] = _%s_%s_I[_k] + _%s_%s_d[_k]/2;\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"%s		}\n",indent);
  fprintf(outfile,"%s	else {\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    fprintf(outfile,"%s		for(unsigned long _k=0; _k<_%s_cross_size*_%s_%s_ncomponents; _k++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s			_%s_%s[_%s_%s_index_pointer_begin + _%s_cross_size*_%s_%s_ncomponents + _k] = _%s_%s_I[_k] + _%s_%s_d[_k];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,fieldName,pXMLString->c_str(),
	    fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"%s		}\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s	}\n",indent);
  fprintf(outfile,"\n");


  fprintf(outfile,"%s// copy I cross vector back into old main cross vector\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s	_%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// copy old main vectors back into last lattice point\n",indent);
  fprintf(outfile,"\n");
  for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {
    fprintf(outfile,"%sfor(unsigned long _j=0; _j<_%s_cross_size*_%s_%s_ncomponents; _j++)\n",
	    indent,fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"%s	_active_%s_%s[_%s_%s_index_pointer_begin + _j] = _%s_%s_old[_j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  // close outer loops

  fprintf(outfile,"%s}\n",indent);
  fprintf(outfile,"\n");

  for(unsigned long i=crossDimNumber(); i>0; i--) {

    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"%s += _%s_dx%li;\n",simulation()->field()->geometry()->dimension(i-1)->name.c_str(),fieldName,i-1);
    fprintf(outfile,"\n");
    for(unsigned long j=0; j<i; j++) {
	fprintf(outfile,"	");
    }
    fprintf(outfile,"}\n");
    fprintf(outfile,"\n");
  }

  // Free allocated memory
  for(list<XMLString>::const_iterator pXMLString = myMainVectorNamesList.begin(); pXMLString != myMainVectorNamesList.end(); pXMLString++) {

    fprintf(outfile,"    delete[] _%s_%s_old;\n",fieldName,pXMLString->c_str());
  }

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

    fprintf(outfile,"    delete[] _%s_%s_I;\n",fieldName,pXMLString->c_str());
    fprintf(outfile,"    delete[] _%s_%s_d;\n",fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");

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


// ******************************************************************************
void xmdsIntegrateSIEX::writeSingleStepCode(
					    FILE *const outfile,
					    const stepCaseEnum& stepCase) const {
  if(debugFlag) {
    printf("xmdsIntegrateSIEX::writeSingleStepCode\n");
  }

  const char *const fieldName = simulation()->field()->name()->c_str();
  const char *const propDim = simulation()->parameters()->propDimName.c_str();

  list<XMLString> tempVectorNamesList;
  tempVectorNamesList.push_back("main");

  const char* indent = "	";
  if(simulation()->parameters()->errorCheck) {
    indent = "		";
  }

  const char* noiseVector = "";
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    noiseVector = ",_noise_vector";
  }

  simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);

  fprintf(outfile,"%s// a_i=a\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%sfor(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s	aifield_%s_main[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
  fprintf(outfile,"\n");

  fprintf(outfile,"%sfor(unsigned long _i1=0; _i1<%li; _i1++) {\n",indent,nIterations());
  fprintf(outfile,"\n");

  if(usesKOperators()) {
    fprintf(outfile,"%s	// calculate the co derivative terms\n",indent);
    fprintf(outfile,"%s	_segment%li_calculate_co_terms();\n",indent,segmentNumber);
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s	// delta a = G[a,t]\n",indent);
  fprintf(outfile,"%s	_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s	if(_i1<%li)\n",indent,nIterations()-1);
  fprintf(outfile,"%s		// a = ai + delta a/2\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%s		for(long _i2=0; _i2<total_local_size*_%s_main_ncomponents; _i2++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%s		for(unsigned long _i2=0; _i2<_%s_size*_%s_main_ncomponents; _i2++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s			_%s_main[_i2] = aifield_%s_main[_i2] + _%s_main[_i2]/2;\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"%s	else\n",indent);
  fprintf(outfile,"%s		// a = ai + delta a\n",indent);
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"%s		for(long _i2=0; _i2<total_local_size*_%s_main_ncomponents; _i2++)\n",indent,fieldName);
  }
  else {
    fprintf(outfile,"%s		for(unsigned long _i2=0; _i2<_%s_size*_%s_main_ncomponents; _i2++)\n",indent,fieldName,fieldName);
  }
  fprintf(outfile,"%s			_%s_main[_i2] = aifield_%s_main[_i2] + _%s_main[_i2];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"%s	}\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s%s += _step/2;\n",indent,propDim);
  fprintf(outfile,"\n");

};