//=======================================
// Binary output of scalar field data
// Latest update: 2014.1.10
//=======================================
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc,char* argv[])
{
    double** point_coord; // three-dimensional coordinate (x,y,z)
    double* field_data; // scalar function f(x,y,z)
    double Lx,Ly,Lz;
    double cx,cy,cz;
    double dx,dy,dz;
    int nx,ny,nz;
    double x,y,z,r;
    int dim; // the number of spatial dimension
    int size; // the number of data points
    int i,j,k,d,count;
    FILE *fp;

    // ---------- Initialization ----------
    dim = 3;

    Lx = 5.0;
    Ly = Lx;
    Lz = Lx;

    cx = Lx/2.0;
    cy = Ly/2.0;
    cz = Lz/2.0;
	
    nx = 51;
    ny = nx;
    nz = nx;
    
    dx = Lx/(nx-1);
    dy = Ly/(ny-1);
    dz = Lz/(nz-1);

    size = nx*ny*nz;

    point_coord = (double **)malloc(size*sizeof(double*));
    for(i=0;i<size;i++){
        point_coord[i] = (double *)malloc(dim*sizeof(double));
    }
    field_data = (double *)malloc(size*sizeof(double));
    // ------------------------------------
    
    // ----- Generation of grid points and field data -----
    count = 0;
    z = 0.0;
    for(k=0;k<nz;k++){
        y = 0.0;
        for(j=0;j<ny;j++){
            x = 0.0;
            for(i=0;i<nx;i++){
                point_coord[count][0] = x; // set x
                point_coord[count][1] = y; // set y
                point_coord[count][2] = z; // set z
                
                r = sqrt((x-cx)*(x-cx)+(y-cy)*(y-cy)+(z-cz)*(z-cz));
                field_data[count] = exp(-r*r); // Set scalar function value. The center point is (Lx/2, Ly/2, Lz/2).
                count++;
                x = x + dx;
            }
            y = y + dy;
        }
        z = z + dz;
    }
    // ----------------------------------------------------
	
	// ---------- Binary output ----------
	fp = fopen("scalar_field_3d.dat","wb");
	if(fp==NULL){
		printf("File cannot be opened.\n");
		return -1;
	}
	
	for(i=0;i<size;i++){
		for(d=0;d<dim;d++){
			fwrite(&point_coord[i][d],sizeof(double),1,fp); // output coordinates
		}
		fwrite(&field_data[i],sizeof(double),1,fp); // output function value
	}
	// -----------------------------------
	
	// ------- Memory free --------
	free(field_data);
	for(i=0;i<size;i++){
		free(point_coord[i]);
	}
	free(point_coord);
	// ----------------------------
	
	return 0;
}