/*
 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: xmdsintegraterk4ip.cc,v 1.15 2005/05/17 06:01:25 sebwuester Exp $
*/

/*! @file xmdsintegraterk4ip.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>

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateRK4IP public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrateRK4IPs=0; //!< The number of xmds integrate RK4IP objects

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

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

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

  if(verbose()) {
    printf("Processing integrate RK4IP element ...\n");
  }
  
  xmdsIntegrate::processElement(yourElement);
  xmdsIntegrateRK4::processElement(yourElement);
  xmdsIntegrateIP::processElement(yourElement);
  
    if(!constantK()) 
    {
    printf("\n");
    printf("WARNING: RK4IP does not work as 4th order algorithm with time dependent k-operators.\n");
    printf("\n");
    }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrateRK4IP private
// ******************************************************************************
// ******************************************************************************

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

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

  xmdsIntegrate::writePrototypes(outfile);
  xmdsIntegrateIP::writePrototypes(outfile);
  xmdsIntegrateRK4::writePrototypes(outfile);
};

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

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

  xmdsIntegrate::writeRoutines(outfile);
  xmdsIntegrateIP::writeRoutines(outfile);
  xmdsIntegrateRK4::writeRoutines(outfile);
};

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

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

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

  const char* kStep="";
  if(!constantK()) {
    kStep="_step/2";
  }

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

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

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

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

  if(usesKOperators()) {
    fprintf(outfile,"%s// a=D[a]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
    fprintf(outfile,"\n");
    simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);
  }


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

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

  if(usesKOperators()) {
    fprintf(outfile,"%s// a_k=D[a_k]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
    fprintf(outfile,"\n");
    simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);
  }

  fprintf(outfile,"%s// a=a+a_k/6\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] += akfield_%s_main[_i1]/6;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

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

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

  fprintf(outfile,"%s// a_k=G[a_k,t+h/2]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/3\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] += akfield_%s_main[_i1]/3;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

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

  fprintf(outfile,"%s// a_k=G[a_k,t+h/2]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

  fprintf(outfile,"%s// a=a+a_k/3\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] += akfield_%s_main[_i1]/3;\n",indent,fieldName,fieldName);
  fprintf(outfile,"\n");

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

  if(usesKOperators()) {
    fprintf(outfile,"%s// a_k=D[a_k]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
    fprintf(outfile,"\n");
    simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);
  }

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

  fprintf(outfile,"%s// a_k=G[a_k,t+h]\n",indent);
  fprintf(outfile,"%s_segment%li_calculate_delta_a(_step,cycle%s);\n",indent,segmentNumber,noiseVector);
  fprintf(outfile,"\n");

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

  if(usesKOperators()) {
    fprintf(outfile,"%s// a=D[a]\n",indent);
    fprintf(outfile,"%s_segment%li_k_propagate(%s);\n",indent,segmentNumber,kStep);
    fprintf(outfile,"\n");
    simulation()->field()->vectors2space(outfile,0,tempVectorNamesList,indent);
  }

  fprintf(outfile,"%s// a=a+a_k/6\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] += akfield_%s_main[_i1]/6;\n",indent,fieldName,fieldName);
};