#include"nanocpp.h"
using namespace std;

int main(int argc, char* argv[]){

  ////////////////////// INIT THE CODE /////////////////
  nanocpp n(argc,argv);
  cout<<n<<endl; //to screen


  /////////////////// MY CODE ////////////////////
  double force=0.;

  //create a volume
  double vol_coords[2*FD_N]={2.1e-6,2.2e-6,1.95e-6,2.05e-6};
  int vol_faces[FD_N]={CE,CE};
  volume<FD_N> vol;
  vol.init(vol_coords,n.get_grid(),vol_faces);

  //create a box
  double mycoords[2*FD_N];
  int myfaces[FD_N*FD_N*FD_N];
  for(int i=0;i<FD_N;i++){
    mycoords[2*i]=vol.get_gsta_coord(i)*n.get_grid().get_dx(i);
    mycoords[2*i+1]=vol.get_gend_coord(i)*n.get_grid().get_dx(i);
      }
  for(int i=0;i<FD_N*FD_N*FD_N;i++)
    myfaces[i]=CE;
  box<FD_N> mybox;
  mybox.init(mycoords,n.get_grid(),myfaces);

  //// Force arrays
  tensor<1,double> fx(&n.total_steps,0.);
  tensor<1,double> gsum(fx);

  //index for force arrays
  int i=0;

  ///////////////////// CYCLE /////////////////////
  double sta=MPI_Wtime();

  for(n.tcurr;n.tcurr<=n.total_steps;n.tcurr++){

    n.time=n.tcurr*n.dt;     //update time;
    n.single_step();    //single step


    /////////////////// MY CODE ////////////////////
    fx[i++]=n.get_steady_optical_force(mybox,0);
    /////////////////////////////////////////////////////

  }


  /////////////////// MY CODE ////////////////////

  //Reduce
  MPI_Reduce(&fx[0],&gsum[0],n.total_steps,MPI_DOUBLE,MPI_SUM,0,n.get_chart().get_comm());
  
  //Average
  if(n.get_chart().get_me()==0){
      string fname=n.get_rootdir()+"fx.bin";
      gsum.save(fname);
      force=0.5*gsum[0];
      for(i=1;i<fx.get_length()-1;i++)
	force+=gsum[i];
      force+=0.5*gsum[n.total_steps-1];
      cout<<"Fx= "<<force*1e12/(1.*n.total_steps)<<"pN"<<endl;
  }

  //////////////////// SAVE OUTPUT AND LOGOUT /////////////

  n.post_process();    //process data
  n.goodbye(sta);   //get time and print goodbye

  ////////////////////////////////////////////////////////


}