#include "dump2sq.h"
/* ####################################################################### */
/* ヘッダーの出力 */
/* ####################################################################### */
void OutputHead(COMMAND *com, HEAD *head, TYPE *type, FILE *f, char *x_string){
  int i;

  /* ヘッダ */
  fprintf(f, "#");
  for (i=0; i<80; i++) fprintf(f, "-");
  fprintf(f, "\n");
  fprintf(f, "# Input file: %s\n", com->infile);
  fprintf(f, "# R min: %g\n", com->GRmin);
  fprintf(f, "# R max: %g\n", com->GRmax);
  fprintf(f, "# R delta: %g\n", com->GRdel);
  fprintf(f, "# Q min: %g\n", com->SQmin);
  fprintf(f, "# Q max: %g\n", com->SQmax);
  fprintf(f, "# Q delta: %g\n", com->SQdel);
  fprintf(f, "# Total steps: %d\n", com->total_step);
  
  fprintf(f, "#");
  for (i=0; i<80; i++) fprintf(f, "-");
  fprintf(f, "\n");			 
  fprintf(f, "# Box length: (%g %g %g)\n", head->lx, head->ly, head->lz);
  fprintf(f, "# Box angle: (%g %g %g)\n", head->alpha, head->beta, head->gamma);
  fprintf(f, "# Atomic number: %d\n", head->atoms);
  fprintf(f, "# Number of density: %g\n", head->rho);
  fprintf(f, "# ALL ATOMIC TYPE\n");
  for (i=0; i<type->types; i++){
    fprintf(f, "# %4s: %d\n", type->elem[i], type->ids[i]);
  }
  fprintf(f, "#");			 
  for (i=0; i<80; i++) fprintf(f, "-");
  fprintf(f, "\n");
  fprintf(f, "# ALL ATOMIC PAIR\n");
  for (i=0; i<type->pairs; i++){
    fprintf(f, "# %2d: %-15s\n", i+1, type->elemij[i]);
  }
  /* 1行 */
  fprintf(f, "#");
  for (i=0; i<80; i++) fprintf(f, "-");
  fprintf(f, "\n");			 
  fprintf(f, "# %12s", x_string);
  for (i=0; i<type->pairs; i++){
    fprintf(f, "%15s", type->elemij[i]);
  }
}


/* ####################################################################### */
/* ファイルに出力 */
/* ####################################################################### */
void OutputCN(COMMAND *com, HEAD *head, TYPE *type, DATASET *cn){
  int i, j;
  FILE *f;
  char outfile[1030];
  double **cn_;

  cn_ = malloc(sizeof(double*) * cn->n);
  for (i=0; i<cn->n; i++) cn_[i] = malloc(sizeof(double) * cn->m);
 
  for (i=1; i<cn->n; i++){
    for (j=0; j<cn->m; j++){
      cn_[i][j] = cn_[i-1][j] + cn->y[i][j];
    }
  }
  for (j=0; j<cn->m; j++) cn_[0][j] = cn->y[0][j];

  sprintf(outfile, "%s.cn", com->outbase);
  f = fopen(outfile, "w");
  fprintf(f, "# Numebr of Coordination\n");
  /* ヘッダ */
  OutputHead(com, head, type, f, "r");
  fprintf(f, "%15s\n", "Total");  
  
  for (i=0; i<cn->n; i++){
    fprintf(f, "%24.16e", cn->x[i]);
    for (j=0; j<cn->m; j++) {
      fprintf(f, " %24.16e", cn_[i][j]);
    }
    fprintf(f, "\n");
  }
  printf("%s was generated.\n", outfile);
  fclose(f);
}

void OutputGr(COMMAND *com, HEAD *head, TYPE *type, DATASET *gr){
  int i, j;
  FILE *f;
  char outfile[1030];

  sprintf(outfile, "%s.gr", com->outbase);
  f = fopen(outfile, "w");
  fprintf(f, "# Radial distribution function\n");
  /* ヘッダ */
  OutputHead(com, head, type, f, "r");
  fprintf(f, "%15s\n", "Total");
  
  for (i=0; i<gr->n; i++){
    fprintf(f, "%24.16e", gr->x[i]);
    for (j=0; j<gr->m; j++) {
      fprintf(f, " %24.16e", gr->y[i][j]);
    }
    fprintf(f, "\n");
  }
  printf("%s was generated.\n", outfile);
  fclose(f);
}

/* ####################################################################### */
/* ファイルに出力 */
/* ####################################################################### */
void OutputSQ(COMMAND *com, HEAD *head, TYPE *type, DATASET *sq){
  int i, j;
  char outfile[1030];
  FILE *f;

  sprintf(outfile, "%s.sq", com->outbase);
  f = fopen(outfile, "w");
  /* S(q)データ */
  fprintf(f, "# XRD and ND S(Q)\n");
  OutputHead(com, head, type, f, "Q");
  fprintf(f, "%15s", "Sq-XRD");
  fprintf(f, "%15s", "Sq-ND");
  fprintf(f, "\n");
  /* データ */
  for (i=0; i<sq->n; i++){
    fprintf(f, "%24.16e", sq->x[i]);
    for (j=0; j<sq->m; j++){
	fprintf(f, " %24.16e", sq->y[i][j]);
    }
    fprintf(f, "\n");
  }
  printf("%s was generated.\n", outfile);
  fclose(f);
}

void OutputCoeff(COMMAND *com, HEAD *head, TYPE *type, DATASET *sq){
  int i, j;
  FILE *f;
  char outfile[1030];

  sprintf(outfile, "%s.coeff", com->outbase);
  f = fopen(outfile, "w");
  fprintf(f, "# Neutron scattering length\n");
  fprintf(f, "# see: https://www.ncnr.nist.gov/resources/n-lengths/\n");
  for (i=0; i<type->types; i++){
    fprintf(f, "%s %lf\n", type->elem[i], type->NDcoeff[i]);
  }
  printf("\n");
  fprintf(f, "# X-ray scattering factor\n");
  fprintf(f, "# see: xcoeff on https://www.szfki.hu/~nphys/rmc++/downloads.html\n");
  fprintf(f, " %12s", "Q");
  for (i=0; i<type->types; i++){
    fprintf(f, " %12s", type->elem[i]);
  }
  fprintf(f, "\n");
  for (i=0; i<sq->n; i++){
    fprintf(f, " %24.16e", sq->x[i]);
    for (j=0; j<type->types; j++){
      fprintf(f, " %24.16e", type->XRDcoeff[i][j]);
    }
    fprintf(f, "\n");
  }
  fclose(f);
  printf("%s was generated.\n", outfile);
}