/*
 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: xmdsintegrateark45ip.cc,v 1.3 2005/05/02 02:35:58 sebwuester Exp $
*/

/*! @file xmdsintegrateark45ip.cc
  @brief Integrate element parsing classes and methods; fourth order Runge-Kutta in the interaction picture

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45IP public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateARK45IPs=0; //!< The number of xmds integrate ARK45IP objects

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

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

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


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


  xmdsIntegrate::processElement(yourElement);
  xmdsIntegrateARK45::processElement(yourElement);
  xmdsIntegrateIP::processElement(yourElement);
  
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45IP protected
// ******************************************************************************



// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateARK45IP private
// ******************************************************************************


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

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARKIP::writeResetRoutine: 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,"// segment %li (ARK45IP) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");

  xmdsIntegrate::writePrototypes(outfile);
  xmdsIntegrateIP::writePrototypes(outfile); 
  
  fprintf(outfile,"// integrate (ARK45IP) prototypes\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"void _segment%li_reset(%s* _reset_to,double _step);\n",segmentNumber,typeName);
  
  xmdsIntegrateARK45::writePrototypes(outfile);
};


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

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

  xmdsIntegrate::writeRoutines(outfile);
  xmdsIntegrateIP::writeRoutines(outfile);
  writeResetRoutine(outfile);   
  xmdsIntegrateARK45::writeRoutines(outfile);
};



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

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

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrateARKIP::writeResetRoutine: 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_reset(%s* _reset_to,double _step){\n",segmentNumber,typeName);
  fprintf(outfile,"\n");
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    fprintf(outfile,"for(long _i1=0; _i1<total_local_size*_%s_main_ncomponents; _i1++)\n",fieldName);
  }
  else {
    fprintf(outfile,"for(unsigned long _i1=0; _i1<_%s_size*_%s_main_ncomponents; _i1++)\n",fieldName,fieldName);
  }
  fprintf(outfile,"    _%s_main[_i1] = _reset_to[_i1];\n",fieldName);
  fprintf(outfile,"\n");	
  
  if(usesKOperators()){    
     fprintf(outfile,"_active_%s_main=_%s_main;\n",fieldName,fieldName);	 
     if(!Smallmemory())
       fprintf(outfile,"_segment%li_k_propagate(-1);\n",segmentNumber);
      else
       fprintf(outfile,"_segment%li_k_propagate(-1.0/5*_step);\n",segmentNumber);       
     }
  fprintf(outfile,"}\n");	
  fprintf(outfile,"\n");				     
					     
}

// ******************************************************************************
void xmdsIntegrateARK45IP::writeSingleStepCode(
					     FILE *const outfile,
					     const stepCaseEnum& stepCase) const {
  if(debugFlag) {
    printf("xmdsIntegrateARK45IP::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 unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();

  const char* indent = "	";

  const char* noiseVector="";
  if((simulation()->parameters()->stochastic)&&(!noNoises())) {
    noiseVector=",_noise_vector";
  }
  if(usesKOperators()){  
    if(!Smallmemory())   
         fprintf(outfile,"%s_segment%li_calculate_k_operator_field(_step);\n",indent,segmentNumber);
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent);
     }
    else
     simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);    

  fprintf(outfile,"%s// a_k=y1\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	akfield_%s_main[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
       fprintf(outfile,"%s// a_i=D(a_2*dt)[y1]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(1);\n",indent,segmentNumber);
       fprintf(outfile,"\n");     
     }else{
       fprintf(outfile,"%s// a_i=D(a_2*dt)[y1]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(1.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
     }   

  fprintf(outfile,"%s// a_i=y2=y1\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_check[_i1] = _%s_main[_i1];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s_active_%s_main = akfield_%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a_k=G[a_k,t]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent);     
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
       fprintf(outfile,"%s// a_i=D(a_2*dt)[a_k]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(1);\n",indent,segmentNumber);
       fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_i=D(a_2*dt)[y1]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(1.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
     }  

  fprintf(outfile,"%s// y1=y1+c_1*a_k\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	_%s_main[_i1] += c[1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_1*a_k\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	_%s_check[_i1] += cs[1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_k = a_i + b_21*a_k\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	akfield_%s_main[_i1] = aifield_%s_main[_i1] + b[2][1]*akfield_%s_main[_i1];\n",indent,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[2]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s// a_k=G[a_k,t+aa_2*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent); 
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s// c_2==cs_2==0\n",indent);
  
    fprintf(outfile,"%s// a_j = d_1*a_i + d_2*y1 + d_3*a_k \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	ajfield_%s_main[_i1] = d[1]*aifield_%s_main[_i1] + d[2]*_%s_main[_i1]+d[3]*akfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s%s += a[3]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s_active_%s_main = ajfield_%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_j=D((a_3-a_2)*dt)[a_j]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(2);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_j=D((a_3-a_2)*dt)[a_j]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(1.0/10*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    }  

  fprintf(outfile,"%s// a_j=G[a_j,t+aa_3*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent); 
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_j=D(-(a_3-a_2)*dt)[a_j]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(-2);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_j=D(-(a_3-a_2)*dt)[a_j]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(-1.0/10*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    }  

  fprintf(outfile,"%s// a_l = e_1*a_i + e_2*y1 + e_3*a_k + e_4*a_j\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	alfield_%s_main[_i1] = e[1]*aifield_%s_main[_i1] + e[2]*_%s_main[_i1] + e[3]*akfield_%s_main[_i1] + e[4]*ajfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// y1=y1+c_3*a_j\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	_%s_main[_i1] += c[3]*ajfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_3*a_j\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	_%s_check[_i1] += cs[3]*ajfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s%s += a[4]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  fprintf(outfile,"%s_active_%s_main = alfield_%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n"); 
 
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D((a_4-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(3);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D((a_4-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(2.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    }  

  fprintf(outfile,"%s// a_l=G[a_l,t+aa_4*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent); 
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D(-(a_4-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(-3);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D(-(a_4-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(-2.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    }  

  fprintf(outfile,"%s// y1=y1+c_4*a_j\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	_%s_main[_i1] += c[4]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_4*a_j\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	_%s_check[_i1] += cs[4]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_l = f_1*a_i + f_2*y1 + f_3*a_k + f_4*a_j +f_5*a_l;\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	alfield_%s_main[_i1] = f[1]*aifield_%s_main[_i1] + f[2]*_%s_main[_i1] + f[3]*akfield_%s_main[_i1] + f[4]*ajfield_%s_main[_i1] + f[5]*alfield_%s_main[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[5]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");

  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D((a_5-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(4);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D((a_5-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(4.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    }  

  fprintf(outfile,"%s// a_l=G[a_l,t+aa_5*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent); 
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D(-(a_5-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(-4);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D(-(a_5-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(-4.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    } 
  
  fprintf(outfile,"%s// c_5==0\n",indent);
  
  fprintf(outfile,"%s// y2=y2+cs_5*a_l\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	_%s_check[_i1] += cs[5]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);

  fprintf(outfile,"%s// a_l= g_1*a_i + g_2*a_k + g_3*a_j + g_4*y_1 +g_5*a_l +g_6*y_2;\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	alfield_%s_main[_i1] = g[1]*aifield_%s_main[_i1] + g[2]*akfield_%s_main[_i1] + g[3]*ajfield_%s_main[_i1] + g[4]*_%s_main[_i1] + g[5]*alfield_%s_main[_i1] + g[6]*_%s_check[_i1] ;\n",indent,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName,fieldName);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s%s += a[6]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D((a_6-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(5);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D((a_6-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(27.0/40*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    } 

  fprintf(outfile,"%s// a_l=G[a_l,t+aa_6*dt]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  if(usesKOperators())     
     simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,indent); 
  fprintf(outfile,"\n");
  
  if(usesKOperators()) 
    if(!Smallmemory()){
    fprintf(outfile,"%s// a_l=D(-(a_6-a_2)*dt)[a_l]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(-5);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// a_l=D(-(a_6-a_2)*dt)[a_l]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(-27.0/40*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    } 
  
  fprintf(outfile,"%s// y1=y1+c_6*a_l\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	_%s_main[_i1] += c[6]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// y2=y2+cs_6*a_l\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	_%s_check[_i1] += cs[6]*alfield_%s_main[_i1];\n",indent,fieldName,fieldName);
  
  fprintf(outfile,"%s// t->t+dt\n",indent); 
  fprintf(outfile,"%s%s -= a[6]*_step;\n",indent,propDim);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s_active_%s_main = _%s_check;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");
  
// these exponents are "accidentally" the same as for "4"
if(usesKOperators()) 
 if(!Smallmemory()){  
    fprintf(outfile,"%s// y2=D((1-a_2)*dt)[y2]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(4);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s//  y2=D((1-a_2)*dt)[y2]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(4.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    } 
  
  fprintf(outfile,"%s_active_%s_main = _%s_main;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");  
if(usesKOperators()) 
 if(!Smallmemory()){  
    fprintf(outfile,"%s// y1=D((1-a_2)*dt)[y1]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(4);\n",indent,segmentNumber);
    fprintf(outfile,"\n");
    }else{
       fprintf(outfile,"%s// y1=D((1-a_2)*dt)[y1]\n",indent);
       fprintf(outfile,"%s_segment%li_k_propagate(4.0/5*_step);\n",indent,segmentNumber);
       fprintf(outfile,"\n"); 
    } 
  
};