#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "math.h"

typedef struct{
  int m;  /* 行 */
  int n;  /* 列 */
  double *x;
  double **y;
}DATA;

void Init(FILE *cn, FILE *coeff, DATA *rdata, DATA *qdata,
	  DATA *cdata){
  char buf[1024];
  int i;
  int rsize, qsize, pairsize;

  while (fgets(buf, sizeof(buf), cn)){
    if (buf[0] != '#') break;
  }
  rsize = 1;
  while (fgets(buf, sizeof(buf), cn)){
    rsize++;
  }
  fgets(buf, sizeof(buf), coeff);
  sscanf(buf, "%d", &qsize);
  fgets(buf, sizeof(buf), coeff);
  sscanf(buf, "%d", &pairsize);

  rdata->m = rsize;
  rdata->n = pairsize;
  rdata->x = malloc(rsize * sizeof(double));
  rdata->y = malloc(rsize * sizeof(double*));
  for (i=0; i<rsize; i++) {
    rdata->y[i] = malloc(pairsize * sizeof(double));
  }
  qdata->m = qsize;
  qdata->n = pairsize;
  qdata->x = malloc(qsize * sizeof(double));
  qdata->y = malloc(qsize * sizeof(double*));
  cdata->m = qsize;
  cdata->n = pairsize;
  cdata->x = malloc(qsize * sizeof(double));
  cdata->y = malloc(qsize * sizeof(double*));
  for (i=0; i<qsize; i++) {
    qdata->y[i] = malloc(pairsize * sizeof(double));
    cdata->y[i] = malloc(pairsize * sizeof(double));
  }
  rewind(cn);
  rewind(coeff);
}

void SetData(FILE *cn, FILE *coeff, DATA *rdata, DATA *cdata){
  int i, j;
  char buf[1024];
  char *tok;
  
  while (fgets(buf, sizeof(buf), cn)){
    if (buf[0] != '#') break;
  }
  for (i=0; i<rdata->m; i++){
    tok = strtok(buf, " ");  // 半角スペース区切りでパース
    rdata->x[i] = atof(tok);
    for (j=0; j<rdata->n; j++) {
      tok = strtok(NULL, " ");
      rdata->y[i][j] = atof(tok);
    }
    fgets(buf, sizeof(buf), cn);
  }
  /* 積算配位数からC(r)に変換 */
  for (i=rdata->m-1; i>0; i--){
    for (j=0; j<rdata->n; j++){
      rdata->y[i][j] = rdata->y[i][j] - rdata->y[i-1][j];
    }
  }
  /* 原子散乱因子 */
  for (i=0; i<3; i++){
    fgets(buf, sizeof(buf), coeff);
  }
  for (i=0; i<cdata->m; i++){
    tok = strtok(buf, " ");  // 半角スペース区切りでパース
    cdata->x[i] = atof(tok);
    for (j=0; j<cdata->n; j++) {
      tok = strtok(NULL, " ");
      cdata->y[i][j] = atof(tok);
    }
    fgets(buf, sizeof(buf), coeff);
  }
}


int main(int argn, char **argv){
  FILE *cn, *coeff;
  DATA rdata, qdata, cdata;
  int i, j, k;
  double q, r;
  double Rc;

  cn = fopen(argv[1], "r");
  coeff = fopen(argv[2], "r");
  Init(cn, coeff, &rdata, &cdata, &qdata);
  SetData(cn, coeff, &rdata, &cdata);
  for (i=0; i<cdata.m; i++) {
    qdata.x[i] = cdata.x[i];
    for (j=0; j<cdata.n; j++) {
      qdata.y[i][j] = 0.0;
    }
  }
  Rc = rdata.x[rdata.m-1];
  
  #pragma omp parallel for private(j, k, q, r)
  for (i=0; i<cdata.m; i++){
    for (j=0; j<cdata.n; j++){
      for (k=0; k<rdata.m; k++){
	q = cdata.x[i];
	r = rdata.x[k];
	qdata.y[i][j] = qdata.y[i][j] + rdata.y[k][j] * cdata.y[i][j] * sin(q*r)/q/r
	                * sin(M_PI*r/Rc)/(M_PI*r/Rc);
      }
    }
  }
  for (i=0; i<qdata.m; i++){
    for (j=0; j<qdata.n; j++){
      printf("%24.16e", qdata.y[i][j]);
    }
    printf("\n");
  }
}