/*
 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: xmdsintegrateark45.cc,v 1.6 2005/08/12 05:46:58 sebwuester Exp $
*/

/*! @file xmdsintegrateark45.cc
  @brief Integrate element parsing classes and methods; fourth order Runge-Kutta

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45 public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateARK45s=0;   //!< Number of xmds integrate ARK45 objects

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

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

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

  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    printf("\n");
    printf("WARNING: ARK45 methods may not always yield correct stochastic convergence.\n");
    printf("\n");
  }
  
  list<XMLString> myXMLStringList;
  list<unsigned long> myULongList;
  
  // ************************************
  // find 'tolerance'

  getAssignmentStrings(yourElement,"tolerance",1,1,myXMLStringList);

  myTolerance=*myXMLStringList.begin();

  if(verbose()) {
    printf("integration tolerance = %s\n",myTolerance.c_str());
  }
  
  // ************************************
  // find 'cutoff'

  getAssignmentStrings(yourElement,"cutoff",0,1,myXMLStringList);


  if(myXMLStringList.size()==1) {
     myCutoff=*myXMLStringList.begin();

    if(verbose()) {
    printf("cutoff  = %s\n",myCutoff.c_str());
      
    }
  }
  else {
      printf("cutoff defaulting to 1e-3 \n");
      myCutoff="1e-3";
       }
       
  // ************************************
  // find 'maximum iterations'

  getAssignmentULongs(yourElement,"max_iterations",0,1,myULongList);

  if(myULongList.size()==1) {
     myMaxIterations=*myULongList.begin();
     
    if(myMaxIterations==0) 
      throw xmdsException(yourElement,"Maximum Iterations must be >= 1 !");

    if(verbose()) {
    printf("Maximum iterations  = %li\n",myMaxIterations);
      
    }
  }
  else {
      if(verbose()) 
          printf("Maximum iterations = infinity \n");
      myMaxIterations=0;
      // this means the feature is disabled 
       }
  
};



// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45 protected
// ******************************************************************************
// ******************************************************************************
const XMLString* xmdsIntegrateARK45::tolerance() const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45::tolerance\n");
  }

  return &myTolerance;
};

// ******************************************************************************
const XMLString* xmdsIntegrateARK45::cutoff() const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45::cutoff\n");
  }

  return &myCutoff;
};

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

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARK45::writePrototypes: cannot find 'main' vector");
  }
  
  const char* typeName="";
  if(mainVector->vectorType()==COMPLEX) {
    typeName="complex";
  }
  else if(mainVector->vectorType()==DOUBLE) {
    typeName="double";
  }
  
  fprintf(outfile,"// integrate (ARK45) prototypes\n");
  fprintf(outfile,"\n");
   
  fprintf(outfile,"double _segment%li_timestep_error(%s* _checkfield);\n",segmentNumber,typeName);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"double _segment%li_setup_sampling(bool* _next_sample_flag,unsigned long* _next_sample_counter);\n",segmentNumber);
  fprintf(outfile,"\n");

  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 xmdsIntegrateARK45::writeRoutines(
				     FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45::writeRoutines\n");
  }

  writeTimestepErrorRoutine(outfile);
  writeSetupSamplingRoutine(outfile);
  writeMainIntegrateRoutine(outfile);

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45 private
// ******************************************************************************
// ******************************************************************************
void xmdsIntegrateARK45::writeTimestepErrorRoutine(
						 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45::writeTimestepErrorRoutine\n");
  }

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

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARK45::writeTimestepErrorRoutine: 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,"double _segment%li_timestep_error(%s* _checkfield) {\n",segmentNumber,typeName);
  fprintf(outfile,"\n");
  fprintf(outfile,"double _error=1e-24;\n");
  fprintf(outfile,"double _result[_%s_main_ncomponents];\n",fieldName);
  
  if(simulation()->field()->geometry()->nDims()>0){
     // finds one peak value for each component of field
     fprintf(outfile,"double _peak[_%s_main_ncomponents];\n",fieldName);
     fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_main_ncomponents;_i0++){\n",fieldName);
     fprintf(outfile,"     _peak[_i0]=0.0;\n");
     fprintf(outfile,"     _result[_i0]=0.0;\n");
     fprintf(outfile,"     }\n");
     fprintf(outfile,"\n");
  }
  fprintf(outfile,"double _temp_error=0.0;\n");
  fprintf(outfile,"double _temp_mod=0.0;\n");
  fprintf(outfile,"\n");
  fprintf(outfile,"unsigned long _%s_main_index_pointer=0;\n",fieldName);
  fprintf(outfile,"\n");

  if(simulation()->field()->geometry()->nDims()>0){
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
         fprintf(outfile,"for(unsigned long _i0=0; _i0<total_local_size; _i0++){\n");
      }
      else {
         fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_size; _i0++){\n",fieldName);
      }
      fprintf(outfile,"   for(unsigned long _i1=0;_i1<_%s_main_ncomponents;_i1++){\n",fieldName);
       if(typeName=="complex")  
          fprintf(outfile,"    _temp_mod=mod2(_%s_main[_%s_main_index_pointer + _i1]);\n",fieldName,fieldName);   
      else 
          fprintf(outfile,"    _temp_mod=fabs(_%s_main[_%s_main_index_pointer + _i1]);\n",fieldName,fieldName);    
    // this construction will add zero if field(x)<peak and correct peak to its new value if field(x)>peak

      fprintf(outfile,"       _peak[_i1]+= 0.5*(_temp_mod-_peak[_i1] + fabs(_temp_mod-_peak[_i1])  );\n");
      fprintf(outfile,"    }\n");  
  
      fprintf(outfile,"       _%s_main_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
      fprintf(outfile,"    }\n");  

      fprintf(outfile,"\n");

     //now the peak value is multiplied with the cutoff so that it is now in fact the amplitude-threshold for error determination
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    
         fprintf(outfile,"MPI_Allreduce(&_peak,&_result,_%s_main_ncomponents,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);\n",fieldName);
         fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_main_ncomponents;_i0++){\n",fieldName);
         fprintf(outfile,"    _peak[_i0]=_result[_i0]*(%s);\n",cutoff()->c_str());
         fprintf(outfile,"    _result[_i0]=0;\n");
         fprintf(outfile,"}\n");
      }
      else {
         fprintf(outfile,"for(unsigned long _i0=0;_i0<_%s_main_ncomponents;_i0++)\n",fieldName); 
         fprintf(outfile,"    _peak[_i0]*=%s;\n",cutoff()->c_str());
      }
  
     fprintf(outfile,"\n");
     fprintf(outfile,"_%s_main_index_pointer=0;\n",fieldName);
  }// end if ndims>0
  fprintf(outfile,"\n");

if(simulation()->field()->geometry()->nDims()>0){ 

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
     fprintf(outfile,"for(long _i0=0; _i0<total_local_size; _i0++){\n");
  }
  else {
     fprintf(outfile,"for(unsigned long _i0=0; _i0<_%s_size; _i0++){\n",fieldName);
  }
  fprintf(outfile,"   for(unsigned long _i1=0;_i1<_%s_main_ncomponents;_i1++)\n",fieldName);

  if(typeName=="complex"){
      fprintf(outfile,"       if(mod2(_%s_main[_%s_main_index_pointer + _i1])>_peak[_i1]){\n",fieldName,fieldName);    
      fprintf(outfile,"            _temp_error=mod(_%s_main[_%s_main_index_pointer + _i1]-_checkfield[_%s_main_index_pointer + _i1])/mod(_%s_main[_%s_main_index_pointer + _i1]);\n",fieldName,fieldName,fieldName,fieldName,fieldName);
      }
      else{
      fprintf(outfile,"       if(fabs(_%s_main[_%s_main_index_pointer + _i1])>_peak[_i1]){\n",fieldName,fieldName);
      fprintf(outfile,"            _temp_error=fabs(_%s_main[_%s_main_index_pointer + _i1]-_checkfield[_%s_main_index_pointer + _i1])/fabs(_%s_main[_%s_main_index_pointer + _i1]);\n",fieldName,fieldName,fieldName,fieldName,fieldName);  
      }  
  fprintf(outfile,"            _error+= 0.5*(_temp_error - _error +fabs(_temp_error - _error) );\n");
  fprintf(outfile,"       }\n");  
  fprintf(outfile,"       _%s_main_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
  fprintf(outfile,"    }\n");
  
 }else{// if ndims==0  
 
  fprintf(outfile,"   for(unsigned long _i1=0;_i1<_%s_main_ncomponents;_i1++){\n",fieldName);
 
  if(typeName=="complex"){   
      fprintf(outfile,"       _temp_error=mod(_%s_main[_i1]-_checkfield[_i1])/mod(_%s_main[_i1]);\n",fieldName,fieldName,fieldName,fieldName,fieldName);
      }
      else{
      fprintf(outfile,"       _temp_error=fabs(_%s_main[_i1]-_checkfield[_i1])/fabs(_%s_main[_i1]);\n",fieldName,fieldName,fieldName,fieldName,fieldName);  
      }  
  fprintf(outfile,"       _error+= 0.5*(_temp_error - _error +fabs(_temp_error - _error) );\n");
  fprintf(outfile,"    }\n");
 }
    
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      fprintf(outfile,"MPI_Allreduce(&_error,&_result,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);\n");
  }
  else {
  fprintf(outfile,"_result[0]=_error;\n");
  }
  fprintf(outfile,"return(_result[0]);\n");        
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n"); 
   
  }

// ******************************************************************************
  
void xmdsIntegrateARK45::writeSetupSamplingRoutine(
						 FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45::writeSetupSamplingRoutine\n");
  }
  
  const char *const propDim = simulation()->parameters()->propDimName.c_str();

  fprintf(outfile,"/* **************************************************/\n"); 
  fprintf(outfile,"double _segment%li_setup_sampling(bool* _next_sample_flag,unsigned long* _next_sample_counter) {\n",segmentNumber);
  fprintf(outfile,"\n");
  // will contain the numbers of all the moment groups that need to be sampled at the next sampling point. The N+1'th entry means "reached end of integration interval". 
  fprintf(outfile,"unsigned long _number_next_mg[%li];\n",simulation()->output()->nMomentGroups()+1);
  // number of mg's that are sampled at the next sampling point
  fprintf(outfile,"unsigned long _number_minima=1;\n");
  
  // The check if n (T_tot/N_samp) = (or<) m (T_tot/M_samp) for two momentgroups, where n,m are _next_sample_counters
  // T_tot is the integration interval and N_samp (M_samp) are the numbers of sampling points, will be replaced
  // by n M_samp = (or<) m N_samp to avoid floating point precision problems.
  
  fprintf(outfile,"unsigned long _previous_m=1;\n");
  fprintf(outfile,"unsigned long _previous_M=1;\n");
  fprintf(outfile,"\n");
  fprintf(outfile,"double _%s_break_next=(double)%s;\n",propDim,interval()->c_str(),lattice());
  fprintf(outfile,"_number_next_mg[0]=%li;\n",simulation()->output()->nMomentGroups());
  // initialize all flags to false
  fprintf(outfile,"for(unsigned long _i0=0; _i0<%li; _i0++)\n",simulation()->output()->nMomentGroups()+1);
  fprintf(outfile,"    _next_sample_flag[_i0]=false;\n");
  fprintf(outfile,"\n");

  // check if moment group needs sampling at the same time as another, already discovered sample (or the final time). If so, add this moment group to the to-be-sampled-list. If moment group demands sampling earlier than all previously noted mg's erase all previous ones from list and set the brekpoint-time to this earlier one.
   for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) 
    if(samples(i)!=0){
      fprintf(outfile,"if(_next_sample_counter[%i]*_previous_M==_previous_m*%li){\n",i,samples(i));
      fprintf(outfile,"     _number_next_mg[_number_minima]=%i;\n",i);
      fprintf(outfile,"     _number_minima++;\n");
      fprintf(outfile,"     }\n");  
      fprintf(outfile,"else if(_next_sample_counter[%i]*_previous_M<_previous_m*%li){\n",i,samples(i));
      fprintf(outfile,"    _%s_break_next=_next_sample_counter[%i]*%.23e;\n",propDim,i,atof(interval()->c_str())/samples(i));
      fprintf(outfile,"    _number_minima=1;\n");
      fprintf(outfile,"    _number_next_mg[0]=%i;\n",i);
      fprintf(outfile,"    _previous_M=%li;\n",samples(i));
      fprintf(outfile,"    _previous_m=_next_sample_counter[%i];\n",i);
      fprintf(outfile,"     }\n");  
      }
  fprintf(outfile,"\n"); 
  //Values of _number_next_mg until _number_minima contain now the complete list of mg's that need to be sampled at the next breakpoint. Set their flags to true.
  fprintf(outfile,"for(unsigned long _i0=0;_i0<_number_minima;_i0++)\n");
  fprintf(outfile,"   _next_sample_flag[_number_next_mg[_i0]]=true;\n");
  fprintf(outfile,"\n"); 
  fprintf(outfile,"return(_%s_break_next);\n",propDim); 
  fprintf(outfile,"}\n"); 
  fprintf(outfile,"\n"); 
    
  }
  

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

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

  const xmdsVector* mainVector;

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

  const char* typeName="";
  if(mainVector->vectorType()==COMPLEX) {
    typeName="complex";
  }
  else if(mainVector->vectorType()==DOUBLE) {
    typeName="double";
  }
  
  bool max_iter=true;
   
  if(myMaxIterations==0 ){ 
      max_iter=false;
      }

  fprintf(outfile,"/* ******************************************** */\n");
  fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
  fprintf(outfile,"\n");
  if((simulation()->parameters()->usempi)&!(simulation()->parameters()->stochastic)){
    fprintf(outfile,"%s *akfield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
    fprintf(outfile,"%s *aifield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
    fprintf(outfile,"%s *ajfield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
    fprintf(outfile,"%s *alfield_%s_main = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
    fprintf(outfile,"\n");
    fprintf(outfile,"%s *_%s_check = new %s[total_local_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName);
  }
  else {
    fprintf(outfile,"%s *akfield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
    fprintf(outfile,"%s *aifield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
    fprintf(outfile,"%s *ajfield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
    fprintf(outfile,"%s *alfield_%s_main = new %s[_%s_size*_%s_main_ncomponents];\n",typeName,fieldName,typeName,fieldName,fieldName);
    fprintf(outfile,"\n");
    
    //_main array will contain 5th order solution and _check array 4yh order
    fprintf(outfile,"%s *_%s_check = 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,"// Cash-Karp coefficients\n");
  fprintf(outfile,"double a_raw[7];\n");
  fprintf(outfile,"double a[7];\n");
  fprintf(outfile,"double b[7][7];\n");
  fprintf(outfile,"double c[7];\n");
  fprintf(outfile,"double cs[7];\n");
  fprintf(outfile,"// linear combinations for the (k_i)s\n");
  fprintf(outfile,"double d[4];\n");
  fprintf(outfile,"double e[5];\n");
  fprintf(outfile,"double f[6];\n");
  fprintf(outfile,"double g[7];\n");
  fprintf(outfile,"\n");
  
  // fractions of _step where the intermediate points are taken
  fprintf(outfile,"a_raw[0]=0.0;\n");
  fprintf(outfile,"a_raw[1]=0.0;\n");
  fprintf(outfile,"a_raw[2]=1.0/5;\n");
  fprintf(outfile,"a_raw[3]=3.0/10;\n");
  fprintf(outfile,"a_raw[4]=3.0/5;\n");
  fprintf(outfile,"a_raw[5]=1.0;\n");
  fprintf(outfile,"a_raw[6]=7.0/8.0;\n");
  fprintf(outfile,"\n");

  // timestep increments calculated from the above
  fprintf(outfile,"a[0]=0.0;\n");
  fprintf(outfile,"a[1]=0.0;\n");
  fprintf(outfile,"for(unsigned long _i0=2;_i0<7;_i0++)\n");
  fprintf(outfile,"  a[_i0]=a_raw[_i0]-a_raw[_i0-1];\n");
  fprintf(outfile,"\n");

  // Cash-Karp coefficient, see numerical recipes
  fprintf(outfile,"b[2][1]=1.0/5;\n");
  fprintf(outfile,"b[3][1]=3.0/40;\n");
  fprintf(outfile,"b[3][2]=9.0/40;\n");
  fprintf(outfile,"b[4][1]=3.0/10;\n");
  fprintf(outfile,"b[4][2]=-9.0/10;\n");
  fprintf(outfile,"b[4][3]=6.0/5;\n");
  fprintf(outfile,"b[5][1]=-11.0/54;\n");
  fprintf(outfile,"b[5][2]=5.0/2;\n");
  fprintf(outfile,"b[5][3]=-70.0/27;\n");
  fprintf(outfile,"b[5][4]=35.0/27;\n");
  fprintf(outfile,"b[6][1]=1631.0/55296;\n");
  fprintf(outfile,"b[6][2]=175.0/512;\n");
  fprintf(outfile,"b[6][3]=575.0/13824;\n");
  fprintf(outfile,"b[6][4]=44275.0/110592;\n");
  fprintf(outfile,"b[6][5]=253.0/4096;\n");
  fprintf(outfile,"\n");

  // for the 5th order solution
  fprintf(outfile,"c[0]=0.0;\n");
  fprintf(outfile,"c[1]=37.0/378;\n");
  fprintf(outfile,"c[2]=0.0;\n");
  fprintf(outfile,"c[3]=250.0/621;\n");
  fprintf(outfile,"c[4]=125.0/594;\n");
  fprintf(outfile,"c[5]=0.0;\n");
  fprintf(outfile,"c[6]=512.0/1771;\n");
  fprintf(outfile,"\n");

  // for the 4th order solution
  fprintf(outfile,"cs[0]=0.0;\n");
  fprintf(outfile,"cs[1]=2825.0/27648;\n");
  fprintf(outfile,"cs[2]=0.0;\n");
  fprintf(outfile,"cs[3]=18575.0/48384;\n");
  fprintf(outfile,"cs[4]=13525.0/55296;\n");
  fprintf(outfile,"cs[5]=277.0/14336;\n");
  fprintf(outfile,"cs[6]=1.0/4;\n");
  fprintf(outfile,"\n");

  // In order to save memory the intermediate results are not 
  // using the b[i][j] and kjs directly. These coefficients are those
  // of the linearcombinations of aifield,ajfield,...,main,check ...
  
  // at the 3rd intermediate step
  fprintf(outfile,"d[0]=0.0;\n");
  fprintf(outfile,"d[1]=1.0-b[3][1]/c[1];\n");
  fprintf(outfile,"d[2]=b[3][1]/c[1];\n");
  fprintf(outfile,"d[3]=b[3][2];\n");
  fprintf(outfile,"\n");

  // at the 4th intermediate step
  fprintf(outfile,"e[0]=0.0;\n");
  fprintf(outfile,"e[1]=1.0-b[4][1]/c[1];\n");
  fprintf(outfile,"e[2]=b[4][1]/c[1];\n");
  fprintf(outfile,"e[3]=b[4][2];\n");
  fprintf(outfile,"e[4]=b[4][3];\n");
  fprintf(outfile,"\n");

  // at the 5th intermediate step
  fprintf(outfile,"f[0]=0.0;\n");
  fprintf(outfile,"f[1]=1.0-b[5][1]/c[1];\n");
  fprintf(outfile,"f[2]=b[5][1]/c[1];\n");
  fprintf(outfile,"f[3]=b[5][2];\n");
  fprintf(outfile,"f[4]=b[5][3]-b[5][1]/c[1]*c[3];\n");
  fprintf(outfile,"f[5]=b[5][4]-b[5][1]/c[1]*c[4];\n");
  fprintf(outfile,"\n");

  // at the 6th intermediate step
  fprintf(outfile,"double _den=c[1]*cs[4]-cs[1]*c[4];\n");
  fprintf(outfile,"g[0]=0.0;\n");
  fprintf(outfile,"g[1]=( b[6][4]*(cs[1]-c[1]) + b[6][1]*(c[4]-cs[4]) )/_den + 1.0;\n");
  fprintf(outfile,"g[2]=  b[6][2];\n");
  fprintf(outfile,"g[3]=( b[6][4]*(cs[1]*c[3] - c[1]*cs[3]) + b[6][1]*(cs[3]*c[4] - c[3]*cs[4]) )/_den + b[6][3];\n");
  fprintf(outfile,"g[4]=( b[6][1]*cs[4]-b[6][4]*cs[1] )/_den;\n");
  fprintf(outfile,"g[5]=  b[6][5] + cs[5]*( b[6][1]*c[4]-b[6][4]*c[1] )/_den;\n");
  fprintf(outfile,"g[6]=( -b[6][1]*c[4]+b[6][4]*c[1] )/_den;\n");
  fprintf(outfile,"\n");

  
  fprintf(outfile,"double _step = %s/(double)%li;\n",interval()->c_str(),lattice());
  fprintf(outfile,"double _min_step=%s;\n",interval()->c_str());
  fprintf(outfile,"double _max_step=0.0;\n");
  fprintf(outfile,"double _tolerance=%s;\n",tolerance()->c_str());
  
  if(simulation()->parameters()->errorCheck) {
  fprintf(outfile,"if(_half_step)\n");
  fprintf(outfile,"    _tolerance=_tolerance/16.0;\n");
  fprintf(outfile,"\n");
  }
  fprintf(outfile,"double _error;\n");
  fprintf(outfile,"bool _discard=false;\n"); 
  fprintf(outfile,"bool _break_next=false;\n"); 
  fprintf(outfile,"bool _next_sample_flag[%li];\n",simulation()->output()->nMomentGroups()+2);
  fprintf(outfile,"for(unsigned long _i0=0;_i0<%li;_i0++)\n",simulation()->output()->nMomentGroups()+2);
  fprintf(outfile,"    _next_sample_flag[_i0]=false;\n");
  fprintf(outfile,"unsigned long _next_sample_counter[%li];\n",simulation()->output()->nMomentGroups());
  fprintf(outfile,"for(unsigned long _i0=0;_i0<%li;_i0++)\n",simulation()->output()->nMomentGroups());
  fprintf(outfile,"    _next_sample_counter[_i0]=1;\n");
  fprintf(outfile,"\n");
  fprintf(outfile,"const double %s_ini=%s;\n",propDim,propDim); 
  fprintf(outfile,"\n");
  fprintf(outfile,"double _%s_break_next=_segment%li_setup_sampling(_next_sample_flag,_next_sample_counter);\n",propDim,segmentNumber);
  fprintf(outfile,"if((%s-%s_ini+_step)>=_%s_break_next){\n",propDim,propDim,propDim);    
  fprintf(outfile,"    _break_next=true;\n");
  fprintf(outfile,"    _step=_%s_break_next-%s+%s_ini;\n",propDim,propDim,propDim);
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");

  if(max_iter){
    fprintf(outfile,"unsigned long _step_counter=0;\n");
    fprintf(outfile,"unsigned long _max_steps;\n");
    if(simulation()->parameters()->errorCheck) {
       fprintf(outfile,"if(_half_step)\n");
       fprintf(outfile,"    _max_steps=%li;\n",2*myMaxIterations);
       fprintf(outfile,"  else\n");
       fprintf(outfile,"    _max_steps=%li;\n",myMaxIterations);
    }else{
       fprintf(outfile,"_max_steps=%li;\n",myMaxIterations);
    }
    fprintf(outfile,"\n");
    }
  
  fprintf(outfile,"do{\n    do{\n");
  fprintf(outfile,"\n");
  
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
  
  if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"    _make_noises(_gen1,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"\n");
    }else{
      fprintf(outfile,"    _make_noises(_gen,_var/_step,_noise_vector,_%s_size*_n_noises);\n",fieldName);
      fprintf(outfile,"\n");        
    }
   }		
   writeSingleStepCode(outfile,FULLSTEP);
   
  fprintf(outfile,"\n");
   
   fprintf(outfile,"    _error=_segment%li_timestep_error(_%s_check);\n",segmentNumber,fieldName);
   fprintf(outfile,"\n");
   fprintf(outfile,"    if(_error<_tolerance){;\n");
   fprintf(outfile,"        if(_step>_max_step)\n");
   fprintf(outfile,"             _max_step=_step;\n");
   fprintf(outfile,"        if(!_break_next)\n");
   fprintf(outfile,"             if(_step<_min_step)\n");
   fprintf(outfile,"                 _min_step=_step;\n");
   fprintf(outfile,"        _step*=0.92*pow(fabs(_tolerance/_error),0.2);\n");
   fprintf(outfile,"        _discard=false;\n");
   fprintf(outfile,"        }\n");
   fprintf(outfile,"    else{\n");
   fprintf(outfile,"        %s-=_step;\n",propDim);
   fprintf(outfile,"        _step*=0.92*pow(fabs(_tolerance/_error),0.25);\n");
   fprintf(outfile,"        _discard=true;\n");
   fprintf(outfile,"        _break_next=false;\n");  
   fprintf(outfile,"        _segment%li_reset(aifield_%s_main,_step);\n",segmentNumber,fieldName);
   fprintf(outfile,"        }\n");

   if(max_iter){
    fprintf(outfile,"\n");
    fprintf(outfile,"    _step_counter++;\n");
    fprintf(outfile,"    if(_step_counter>= _max_steps){\n");
    fprintf(outfile,"        _discard=false;\n");
    fprintf(outfile,"        _break_next=true;\n");
    fprintf(outfile,"        _next_sample_flag[%li]=true;\n",simulation()->output()->nMomentGroups()+1);
    fprintf(outfile,"    }\n");
    }
   
   fprintf(outfile,"    }while(_discard);\n");


   fprintf(outfile,"\n");
   fprintf(outfile,"    if(_break_next){\n");
   for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) 
    if(samples(i)!=0){    
      fprintf(outfile,"        if(_next_sample_flag[%li]){\n",i);       
      fprintf(outfile,"    		   _mg%li_sample();\n",i);                   
      fprintf(outfile,"    		   _next_sample_counter[%li]++;\n",i);
      fprintf(outfile,"        }\n");
      }
   fprintf(outfile,"        if(_next_sample_flag[%li])\n",simulation()->output()->nMomentGroups());       
   fprintf(outfile,"    		   _next_sample_flag[%li]=true;\n",simulation()->output()->nMomentGroups()+1);  
   fprintf(outfile,"             else{\n");     
   fprintf(outfile,"               _break_next=false;\n");
   fprintf(outfile,"               _%s_break_next=_segment%li_setup_sampling(_next_sample_flag,_next_sample_counter);\n",propDim,segmentNumber);
   fprintf(outfile,"                 }\n");    
   fprintf(outfile,"    }\n");
   fprintf(outfile,"    if((%s-%s_ini+_step)>_%s_break_next){\n",propDim,propDim,propDim);    
   fprintf(outfile,"        _break_next=true;\n");
if(simulation()->parameters()->usempi) 
   fprintf(outfile,"    	if(rank==0)\n");
   
   // calculate the previous timestep which was actually accepted while the current
   // value of _step is only the guess for the next iteration.
   fprintf(outfile,"    	    printf(\"Current timestep: %%e\\n\",_step/(0.92*pow(fabs(_tolerance/_error),0.2)));\n");  
   fprintf(outfile,"        _step=_%s_break_next-%s+%s_ini;\n",propDim,propDim,propDim);
   fprintf(outfile,"    }\n");
   fprintf(outfile,"}while(!_next_sample_flag[%li]);\n",simulation()->output()->nMomentGroups()+1);
   fprintf(outfile,"\n"); 
   fprintf(outfile,"printf(\"Segment %li: minimum timestep: %%e maximum timestep: %%e \\n\",_min_step,_max_step);\n",segmentNumber);      
    fprintf(outfile,"\n");
    
if(max_iter){
    fprintf(outfile,"if(_step_counter>= _max_steps){\n");
    fprintf(outfile,"    printf(\" \\n \");\n");
    if(simulation()->parameters()->errorCheck) {
        fprintf(outfile,"    if(_half_step)\n");
        fprintf(outfile,"        printf(\"Reached %li iterations, exiting at %s = %%e \\n \",%s);\n",2*myMaxIterations,propDim,propDim);
        fprintf(outfile,"    else\n");
        fprintf(outfile,"        printf(\"Reached %li iterations, exiting at %s = %%e \\n \",%s);\n",myMaxIterations,propDim,propDim);
        }else{
        fprintf(outfile,"    printf(\"Reached %li iterations, exiting at %s = %%e \\n \",%s);\n",myMaxIterations,propDim,propDim);
        }
    fprintf(outfile,"    printf(\"Last error: %%e \\n \",_error);\n");
    fprintf(outfile,"    printf(\"Last planned timestep: %%e \\n \",_step);\n");
    fprintf(outfile,"    printf(\" \\n \");\n");
    fprintf(outfile,"}\n"); 
}//end max iter   

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

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

// ******************************************************************************
void xmdsIntegrateARK45::writeCalculateCrossFieldRoutine(
						       FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45IP::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() {\n",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 xmdsIntegrateARK45::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 xmdsIntegrateARK45::writeCalculateCrossFieldRoutine: cannot find cross vector");
    }

    crossVectorList.push_back(crossVector);

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

    fprintf(outfile,"%s *_%s_%s_K = 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_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 K and 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_K[_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,"%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,"%s	}\n",indent);
    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// ********** step 1 ***************\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s	{\n",indent);
  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 k1\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	}\n",indent);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// add to K 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_K[_j] += _%s_%s_d[_j]/6;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// create next cross 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[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j]/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// create midpoint main vectors for steps 2 and 3\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 for steps 2 and 3\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,"%s// ********** step 2 ***************\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s	{\n",indent);
  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 k1\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	}\n",indent);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// add to K 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_K[_j] += _%s_%s_d[_j]/3;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// create next cross 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[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j]/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// ********** step 3 ***************\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s	{\n",indent);
  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 k1\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	}\n",indent);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// add to K 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_K[_j] += _%s_%s_d[_j]/3;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// create next cross 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[_%s_%s_index_pointer_begin + _j] = _%s_%s_I[_j] + _%s_%s_d[_j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }

  fprintf(outfile,"%s// move main vectors to next lattice point for step 4\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[_%s_cross_size*_%s_%s_ncomponents + _j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),
	    fieldName,fieldName,pXMLString->c_str());
    fprintf(outfile,"\n");
  }
  fprintf(outfile,"%s// move cross dim to next lattice point for step 4\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,"%s// ********** step 4 ***************\n",indent);
  fprintf(outfile,"\n");
  fprintf(outfile,"%s	{\n",indent);
  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 k1\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	}\n",indent);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// add to K 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_K[_j] += _%s_%s_d[_j]/6;\n",indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str());
    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 K cross vector into next 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 + _%s_cross_size*_%s_%s_ncomponents + _j] = _%s_%s_K[_j];\n",
	    indent,fieldName,pXMLString->c_str(),fieldName,pXMLString->c_str(),fieldName,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);

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

    fprintf(outfile,"\n");
    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");
  }


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

    fprintf(outfile,"    delete[] _%s_%s_old;\n", fieldName, pXMLString->c_str());
    fprintf(outfile,"\n");
  }
	
  for (list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {

    fprintf(outfile, "delete[] _%s_%s_K;\n", fieldName, pXMLString->c_str());
    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");
  fprintf(outfile,"\n");
};