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

/*! @file xmdsintegrateip.cc
  @brief Integrate element parsing classes and methods; interaction picture

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateIP public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateIPs=0;  //!< Number of integrate IP objects

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

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

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

  // ************************************
  // parse code for operators

  const xmdsVector* mainVector;

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

  XMLString* theCode = propagationCode();

  XMLString nextOperatorName;
  XMLString nextComponentName;
  unsigned long start=0;
  unsigned long end=0;

  while(end<theCode->length()) {
    if(findNextcoPair(nextOperatorName,nextComponentName,start,end)) {

      unsigned long nextComponentNumber;

      if(!mainVector->getComponent(nextComponentName,nextComponentNumber)) {
	sprintf(errorMessage(),"[%s] is not a component of the main vector",nextComponentName.c_str());
	throw xmdsException(yourElement,errorMessage());
      }

      unsigned long nextOperatorNumber;

      if(!getKOperator(nextOperatorName,nextOperatorNumber)) {
	sprintf(errorMessage(),"'%s' was not defined in <k_operators>",nextOperatorName.c_str());
	throw xmdsException(yourElement,errorMessage());
      }


      unsigned long dummycoKey;

      if(getcoKey(nextComponentNumber,nextOperatorNumber,dummycoKey)) {
	sprintf(errorMessage(),"'%s[%s]' used multiply",nextOperatorName.c_str(),nextComponentName.c_str());
	throw xmdsException(yourElement,errorMessage());
      }

      if(verbose()) {
	printf("adding operator-component pair: %s[%s]\n",nextOperatorName.c_str(),nextComponentName.c_str());
      }

      addcoPair(nextComponentNumber,nextOperatorNumber);

      theCode->replaceData(start,end-start,"0");

      start = start+1;
    }
    else {
      start = end;
    }
  }
       
};

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

  if(usesKOperators()) {
   if(AdaptiveIP()){
   
      if(Smallmemory()){
      fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
      fprintf(outfile,"	const double& _step);\n");
      fprintf(outfile,"\n");
      }
      else{
      fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step);\n",segmentNumber);
      fprintf(outfile,"\n");
      fprintf(outfile,"void _segment%li_k_propagate(int exponent);\n",segmentNumber);
      fprintf(outfile,"\n");             
         } 

   }else{
    if(constantK()) {
      fprintf(outfile,"void _segment%li_calculate_k_operator_field();\n",segmentNumber);
      fprintf(outfile,"\n");
      fprintf(outfile,"void _segment%li_k_propagate();\n",segmentNumber);
      fprintf(outfile,"\n");
    }
    else {
      fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
      fprintf(outfile,"	const double& _step);\n");
      fprintf(outfile,"\n");
    }
   }
  }
};

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

  if(usesKOperators()) {
   if(AdaptiveIP()){
      if(Smallmemory())
        writeTimeDepkPropagateRoutine(outfile);  // use conventional routine
      else{
         if(constantK())                     // optimized for speed, but needs lots of memory
             writeARK45CalculatekOperatorFieldRoutine(outfile);     
            else                         // optimized for speed, but needs even more memory
             writeARK45TimeDepkCalculatekOperatorFieldRoutine(outfile);
         writeARK45kPropagateRoutine(outfile);              
         }   
   }else{  
    if(constantK()) {
      writeCalculatekOperatorFieldRoutine(outfile);
      writeConstkPropagateRoutine(outfile);
    }
    else
      writeTimeDepkPropagateRoutine(outfile);
    }
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateIP private
// ******************************************************************************
// ******************************************************************************

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

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

  unsigned long i;

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_k_operator_field() {\n",segmentNumber);
  fprintf(outfile,"\n");
  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");

  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
      fprintf(outfile,"_segment%li_k_operator_field = new complex[total_local_size*_segment%li_nkoperators];\n",
	      segmentNumber,segmentNumber);
      fprintf(outfile,"\n");
    }

  fprintf(outfile,"double _step=(%s/(double)%li)/2;\n",interval()->c_str(),lattice());
  fprintf(outfile,"\n");

  if(simulation()->parameters()->errorCheck) {
    fprintf(outfile,"if(_half_step)\n");
    fprintf(outfile,"  _step=_step/2;\n");  
    fprintf(outfile,"\n");
  }

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");

  fprintf(outfile,"unsigned long _k_operator_index_pointer=0;\n");

  simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());

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

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"%s_segment%li_k_operator_field[_k_operator_index_pointer + %li].re = exp(%s.re*_step)*cos(%s.im*_step);\n",
	    indent,segmentNumber,i,KOperator(i)->c_str(),KOperator(i)->c_str());
    fprintf(outfile,"%s_segment%li_k_operator_field[_k_operator_index_pointer + %li].im = exp(%s.re*_step)*sin(%s.im*_step);\n",
	    indent,segmentNumber,i,KOperator(i)->c_str(),KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
  simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());

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


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

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

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

  const xmdsVector* mainVector;

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_k_propagate() {\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"double _temp;\n");
  fprintf(outfile,"unsigned long _segment%li_kop_index_pointer=0;\n",segmentNumber);
  fprintf(outfile,"unsigned long _active_%s_index_pointer=0;\n",fieldName);
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,"");

  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,"\n");

  for(unsigned long i=0; i<mainVector->nComponents(); i++) {

    const coStruct* thecoStruct;

    if(getcoStruct(i,thecoStruct)) {

      for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {

	fprintf(outfile,"	_temp = _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
		segmentNumber,segmentNumber,*pULong);
	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].re\n",fieldName,fieldName,i);
	fprintf(outfile,"	      - _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,segmentNumber,*pULong);
	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].im;\n",fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"	_active_%s_main[_active_%s_index_pointer + %li].im =\n",fieldName,fieldName,i);
	fprintf(outfile,"	     _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].re\n",
		segmentNumber,segmentNumber,*pULong);
	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].im\n",fieldName,fieldName,i);
	fprintf(outfile,"	   + _segment%li_k_operator_field[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,segmentNumber,*pULong);
	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].re;\n",fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"	_active_%s_main[_active_%s_index_pointer + %li].re=_temp;\n",fieldName,fieldName,i);
	fprintf(outfile,"\n");
      }
    }
  }
  fprintf(outfile,"  _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
  fprintf(outfile,"  _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
  fprintf(outfile,"	};\n");
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

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

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

  const xmdsVector* mainVector;

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_k_propagate(\n",segmentNumber);
  fprintf(outfile,"	const double& _step) {\n");
  fprintf(outfile,"\n");
  for(unsigned long i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"double _temp;\n");
  fprintf(outfile,"double _temp2;\n");
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");

  simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());

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

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");

  for(unsigned long i=0;i<nKOperators();i++) {
    fprintf(outfile,"%s_temp2 = exp(%s.re*_step);\n",indent,KOperator(i)->c_str());
    fprintf(outfile,"%s_temp = _temp2*cos(%s.im*_step);\n",indent,KOperator(i)->c_str());
    fprintf(outfile,"%s%s.im = _temp2*sin(%s.im*_step);\n",indent,KOperator(i)->c_str(),KOperator(i)->c_str());
    fprintf(outfile,"%s%s.re = _temp;\n",indent,KOperator(i)->c_str());
    fprintf(outfile,"\n");
  }
  for(unsigned long i=0; i<mainVector->nComponents(); i++) {

    const coStruct* thecoStruct;

    if(getcoStruct(i,thecoStruct)) {

      for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {
	
	fprintf(outfile,"%s_temp =  %s.re*_active_%s_main[_%s_main_index_pointer + %li].re\n",indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
	fprintf(outfile,"%s       - %s.im*_active_%s_main[_%s_main_index_pointer + %li].im;\n",indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li].im =  %s.re*_active_%s_main[_%s_main_index_pointer + %li].im\n",
		indent,fieldName,fieldName,i,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
	fprintf(outfile,"%s	                                            + %s.im*_active_%s_main[_%s_main_index_pointer + %li].re;\n",
		indent,KOperator(*pULong)->c_str(),fieldName,fieldName,i);
	fprintf(outfile,"\n");
	fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li].re=_temp;\n",indent,fieldName,fieldName,i);
	fprintf(outfile,"\n");
      }
    }
  }

  simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());

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

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

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

  unsigned long i;

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step) {\n",segmentNumber);
  fprintf(outfile,"\n");
  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"double factor1=1.0/5;   // a_raw[2]\n");
  fprintf(outfile,"double factor2=1.0/10;  // a_raw[3]-a_raw[2]\n");
  fprintf(outfile,"double factor3=2.0/5;   // a_raw[4]-a_raw[2]\n");
  fprintf(outfile,"double factor4=4.0/5;   // a_raw[5]-a_raw[2]\n");
  fprintf(outfile,"double factor5=27.0/40; // a_raw[6]-a_raw[2]\n");  
  

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");

  fprintf(outfile,"unsigned long _k_operator_index_pointer=0;\n");

  simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());

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

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");
  
  for(int n=1;n<6;n++){
  for(i=0;i<nKOperators();i++) {
fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].re = exp(%s.re*factor%i*_step)*cos(%s.im*factor%i*_step);\n",
	    indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
    fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].im = exp(%s.re*factor%i*_step)*sin(%s.im*factor%i*_step);\n",
	    indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
  }
  fprintf(outfile,"\n");  
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
  simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());

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


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

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

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

  const xmdsVector* mainVector;

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_k_propagate(int exponent) {\n",segmentNumber);
  fprintf(outfile,"\n");
  fprintf(outfile,"double _temp;\n");
  fprintf(outfile,"unsigned long _segment%li_kop_index_pointer=0;\n",segmentNumber);
  fprintf(outfile,"unsigned long _active_%s_index_pointer=0;\n",fieldName);
  fprintf(outfile,"\n");

  simulation()->field()->vectors2space(outfile,fullSpace,tempVectorNamesList,"");

  fprintf(outfile,"if(exponent>0){\n");
  fprintf(outfile,"switch(exponent){\n");
  
  for(int n=1;n<6;n++){
      fprintf(outfile,"case %i:\n",n);  

      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,"\n");

      for(unsigned long i=0; i<mainVector->nComponents(); i++) {

        const coStruct* thecoStruct;

        if(getcoStruct(i,thecoStruct)) {

          for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {

    	fprintf(outfile,"	_temp = _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].re\n",
    		segmentNumber,n,segmentNumber,*pULong);
    	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].re\n",fieldName,fieldName,i);
    	fprintf(outfile,"	      - _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,n,segmentNumber,*pULong);
    	fprintf(outfile,"				*_active_%s_main[_active_%s_index_pointer + %li].im;\n",fieldName,fieldName,i);
    	fprintf(outfile,"\n");
    	fprintf(outfile,"	_active_%s_main[_active_%s_index_pointer + %li].im =\n",fieldName,fieldName,i);
    	fprintf(outfile,"	     _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].re\n",
		segmentNumber,n,segmentNumber,*pULong);
    	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].im\n",fieldName,fieldName,i);
    	fprintf(outfile,"	   + _segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li].im\n",
		segmentNumber,n,segmentNumber,*pULong);
    	fprintf(outfile,"			*_active_%s_main[_active_%s_index_pointer + %li].re;\n",fieldName,fieldName,i);
    	fprintf(outfile,"\n");
    	fprintf(outfile,"	_active_%s_main[_active_%s_index_pointer + %li].re=_temp;\n",fieldName,fieldName,i);
    	fprintf(outfile,"\n");
          }
        }
      }
      fprintf(outfile,"  _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
      fprintf(outfile,"  _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
      fprintf(outfile,"	};\n");
      fprintf(outfile,"	break;\n");
  }
  fprintf(outfile,"}\n");
  fprintf(outfile,"} else {\n");
  fprintf(outfile,"exponent=abs(exponent);\n");
  fprintf(outfile,"switch(exponent){\n");
  
  for(int n=1;n<6;n++){
      fprintf(outfile,"case %i:\n",n);  

      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,"\n");

      for(unsigned long i=0; i<mainVector->nComponents(); i++) {

        const coStruct* thecoStruct;

        if(getcoStruct(i,thecoStruct)) {

          for(list<unsigned long>::const_iterator pULong = thecoStruct->operatorNumbersList.begin(); pULong != thecoStruct->operatorNumbersList.end(); pULong++) {

    	fprintf(outfile,"	_active_%s_main[_active_%s_index_pointer + %li]/=_segment%li_k_operator_field_%i[_segment%li_kop_index_pointer + %li];\n",fieldName,fieldName,i,segmentNumber,n,segmentNumber,*pULong);
    	fprintf(outfile,"\n");
          }
        }
      }
      fprintf(outfile,"  _segment%li_kop_index_pointer+=_segment%li_nkoperators;\n",segmentNumber,segmentNumber);
      fprintf(outfile,"  _active_%s_index_pointer+=_%s_main_ncomponents;\n",fieldName,fieldName);
      fprintf(outfile,"	};\n");
      fprintf(outfile,"	break;\n");
  }
  fprintf(outfile,"}\n");
  fprintf(outfile,"}\n");  
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

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

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const char *const propDim = simulation()->parameters()->propDimName.c_str();  
  const unsigned long fullSpace = simulation()->field()->geometry()->fullSpace();

  unsigned long i;

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_k_operator_field(double _step) {\n",segmentNumber);
  fprintf(outfile,"\n");
  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"complex %s;\n",KOperator(i)->c_str());
  }
  fprintf(outfile,"\n");
  fprintf(outfile,"double factor1=1.0/5;   // a_raw[2]\n");
  fprintf(outfile,"double factor2=1.0/10;  // a_raw[3]-a_raw[2]\n");
  fprintf(outfile,"double factor3=2.0/5;   // a_raw[4]-a_raw[2]\n");
  fprintf(outfile,"double factor4=4.0/5;   // a_raw[5]-a_raw[2]\n");
  fprintf(outfile,"double factor5=27.0/40; // a_raw[6]-a_raw[2]\n");  
  fprintf(outfile,"double _store_time=t;\n");  

  fprintf(outfile,"\n");
  fprintf(outfile,"double _store_%s=%s; \n",propDim,propDim);  
  fprintf(outfile,"double _d%s1=0.0;   \n",propDim); 
  fprintf(outfile,"double _d%s2=3.0/10*_step;   \n",propDim);  
  fprintf(outfile,"double _d%s3=3.0/10*_step;  \n",propDim);  
  fprintf(outfile,"double _d%s4=2.0/5*_step;   \n",propDim);  
  fprintf(outfile,"double _d%s5=-1.0/8*_step;   \n",propDim); 
  fprintf(outfile,"\n"); 

  simulation()->field()->vectors2space(outfile,fullSpace,*KVectorNamesList(),"");
  
  fprintf(outfile,"unsigned long _k_operator_index_pointer;\n");
  
for(int n=1;n<6;n++){
  fprintf(outfile,"_k_operator_index_pointer=0;\n");
  fprintf(outfile,"%s+=_d%s%i;\n",propDim,propDim,n);
  fprintf(outfile,"\n");
  fprintf(outfile,"{\n"); //make counters in loops local

  simulation()->field()->openLoops(outfile,fullSpace,*KVectorNamesList());

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

  fprintf(outfile,"\n");
  fprintf(outfile,"// ********** Code from k_operators *************\n");
  fprintf(outfile,"%s\n",KOperatorsCode()->c_str());
  fprintf(outfile,"// **********************************************\n");
  fprintf(outfile,"\n");
  

  for(i=0;i<nKOperators();i++) {
    fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].re = exp(%s.re*factor%i*_step)*cos(%s.im*factor%i*_step);\n",
	    indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
    fprintf(outfile,"%s_segment%li_k_operator_field_%i[_k_operator_index_pointer + %li].im = exp(%s.re*factor%i*_step)*sin(%s.im*factor%i*_step);\n",
	    indent,segmentNumber,n,i,KOperator(i)->c_str(),n,KOperator(i)->c_str(),n);
  }

  fprintf(outfile,"\n");
  fprintf(outfile,"%s_k_operator_index_pointer += _segment%li_nkoperators;\n",indent,segmentNumber);
  simulation()->field()->closeLoops(outfile,fullSpace,*KVectorNamesList());
  fprintf(outfile,"}\n");// end-make counters in loops local

  fprintf(outfile,"\n");  
  }

  fprintf(outfile,"%s=_store_%s;\n",propDim,propDim);  
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};