/* ============================================================================================================ */
/* remove_par_file.c: write remove particle information to file                                                 */
/* This program was made in reference to accept.c                                                               */
/*     remove_par_file(ECS, "file name")                                                                        */
/* made by Kenichi Kaneta (email : kaneta@cics.jp)                                                */
/* ============================================================================================================ */

#include <stdio.h>
#include <math.h>
#include "elem.h"
#define CHUNKCOUNT 65536

enum state {PARNEW, PARKIL};

struct buffparinfo{
  enum state state;
  int ID;
  double  Wr[3],  GBr[3],  G,  time; /* start  informations of particles */
  double kWr[3], kGBr[3], kG, ktime; /* remove informations of particles */
};

struct remove_par_file_info{
  struct buffparinfo *buff;
  unsigned int lastpar;
  gpttransform outtf;                /* Transformation to output particles */
  char *file;
};

static int  remove_par_file_out_usr(double t, double *dt, const double *x, void *vinfo );
static void remove_par_file_end_usr(double tstart, double tend, double dt, const double *xstart, const double *xend, void *vinfo, void *stepinfo);
static void remove_par_file_ter(void *vinfo );

/* ============================================================================================================ */
/*  initial routine                                                                                             */
/* ============================================================================================================ */
void remove_par_file_init(gptinit *init){
  struct remove_par_file_info *info;

  gptbuildECS( init );

  if( gptgetargnum(init)!=1 ){
    gpterror( "Syntax: %s(ECS)\n", gptgetname(init) );
  }

  /* Allocate and initialize info structure */
  info = (struct remove_par_file_info *)gptmalloc( sizeof(struct remove_par_file_info) );
  info->buff    = (struct buffparinfo *)gptmallocitemarray(sizeof(struct buffparinfo),CHUNKCOUNT,1);
  info->lastpar = 0;
  gptconcattransform( &info->outtf, &init->axis->a, &init->e ) ;

  info->file = (char*) gptgetargstring(init,1);

  /* Register ODE functions */
  odeaddoutfunction(ODEFNC_USR, remove_par_file_out_usr, info);  /* Output fuction. Derivative information is already present. */
  odeaddendfunction(ODEFNC_USR, remove_par_file_end_usr, info);  /* Initialization after a successful timestep. */
  gptaddmainfunction(GPTMAINFNC_TER, remove_par_file_ter,info); /* Termination routines, this fuction will be called */
}

/* ============================================================================================================ */
/*   Store new particle                                                                                         */
/* ============================================================================================================ */
static int remove_par_file_out_usr( double t, double *dt, const double *x, void *vinfo ){
  struct remove_par_file_info *info = (struct remove_par_file_info *)vinfo;
  int i;

  for(i=0; i<numpar; i++){
    if( (unsigned int)pars[i].ID > info->lastpar ){
      info->buff = (struct buffparinfo *)gptmallocitem(info->buff,sizeof(struct buffparinfo),CHUNKCOUNT,info->lastpar);

      gpttoUCS(&info->outtf,pars[i].Wr ,info->buff[pars[i].ID-1].Wr );
      gptdirectiontoUCS(&info->outtf,pars[i].GBr,info->buff[pars[i].ID-1].GBr);
      info->buff[pars[i].ID-1].G     = pars[i].G;
      info->buff[pars[i].ID-1].time  = t;
      info->buff[pars[i].ID-1].state = PARNEW;
      info->buff[pars[i].ID-1].ID    = pars[i].ID;
      info->lastpar = pars[i].ID;
    }    
  }

  return(0);
}


/* ============================================================================================================ */
/*  Mark deleted particles and store particle infomation                                                        */
/* ============================================================================================================ */
static void remove_par_file_end_usr( double tstart, double tend, double dt, const double *xstart, const double *xend, void *vinfo, void *stepinfo ){
  struct remove_par_file_info *info = (struct remove_par_file_info *)vinfo;
  int i;

  for(i=0; i<numpar; i++){
    if(pars[i].tokill){
      gpttoUCS(&info->outtf,pars[i].Wr ,info->buff[pars[i].ID-1].kWr );
      gptdirectiontoUCS(&info->outtf,pars[i].GBr,info->buff[pars[i].ID-1].kGBr);
      info->buff[pars[i].ID-1].kG    = pars[i].G;
      info->buff[pars[i].ID-1].ktime = tend;
      info->buff[pars[i].ID-1].state = PARKIL;
    }
  }
}

/* ============================================================================================================ */
/*  write remove particle infomation to file                                                                    */
/* ============================================================================================================ */
static void remove_par_file_ter( void *vinfo ){
  struct remove_par_file_info *info = (struct remove_par_file_info *)vinfo;
  unsigned int i;
  FILE *fp;
  errno_t error;

  /*  if((fp = fopen(info->file, "w")) == NULL){*/
  if((error = fopen_s(&fp, info->file, "w")) != 0){
    printf("file open error!!  at  remove_par_file\n");
    return;
  }

  fprintf(fp, "#");
  fprintf(fp, "ID\t");
  fprintf(fp, "x\t");
  fprintf(fp, "y\t");
  fprintf(fp, "z\t");
  fprintf(fp, "beta_x\t");
  fprintf(fp, "beta_y\t");
  fprintf(fp, "beta_z\t");
  fprintf(fp, "G\t");
  fprintf(fp, "time\n");
    
  for(i=0 ; i<info->lastpar ; i++){
    if(info->buff[i].state==PARKIL){
      fprintf(fp, "%d\t", info->buff[i].ID);
      fprintf(fp, "%e\t", info->buff[i].kWr[0]);
      fprintf(fp, "%e\t", info->buff[i].kWr[1]);
      fprintf(fp, "%e\t", info->buff[i].kWr[2]);
      fprintf(fp, "%e\t", info->buff[i].kGBr[0]/info->buff[i].kG);
      fprintf(fp, "%e\t", info->buff[i].kGBr[1]/info->buff[i].kG);
      fprintf(fp, "%e\t", info->buff[i].kGBr[2]/info->buff[i].kG);
      fprintf(fp, "%e\t", info->buff[i].kG);
      fprintf(fp, "%e\n", info->buff[i].ktime);
    }
  }

  fclose(fp);

  gptfreeitemarray(info->buff,sizeof(struct buffparinfo),CHUNKCOUNT,info->lastpar);
  gptfree(info);

}