130 lines
3.0 KiB
C
130 lines
3.0 KiB
C
#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");
|
|
}
|
|
}
|