#include <stdlib.h>
#include <stdio.h>
#include <math.h>


void spheremaker(float r, float* c, int dim, float* mat);
void centerofmass(int dim, float* mat, float* cm);
void xyzprinter(float* matrix, int dim, FILE* file);


int main(){

int i;
FILE* fp;
float* matrix;
float* com;

int np=10000;
float c[]={0.0,0.0,0.0};
float rc=10.0;

/*INITIALIZING*/
fp=fopen("data.dat","w+");
matrix=(float*)malloc(np*3*sizeof(float));
com=(float*)malloc(3*sizeof(float));

/*CORE*/
    spheremaker(rc,c,np,matrix);
    centerofmass(np,matrix,com);
    xyzprinter(matrix,np,fp);
    xyzprinter(com,1,fp);

/*OUTITIALIZING*/
fclose(fp);
free((void*)matrix);
free((void*)com);
exit(0);
}




void spheremaker(float r, float* c, int dim, float* matrix){
float x,y,z,norm;
int i=0;

while(i<dim)
    {

    x=2.0*rand()/RAND_MAX-1;
    y=2.0*rand()/RAND_MAX-1;
    z=2.0*rand()/RAND_MAX-1;

    norm=pow(x,2)+pow (y,2)+pow(z,2);
    norm=sqrt(norm);
    if(norm<1){
        matrix[0+3*i]=r*x+c[0];
        matrix[1+3*i]=r*y+c[1];
        matrix[2+3*i]=r*z+c[2];
        i=i+1;
        }
    }
}




void centerofmass(int dim, float* matrix, float* com){
int i,j;

for(j=0;j<3;j=j+1){
    com[j]=0;
    for(i=0;i<dim;i=i+1){
        com[j]=com[j]+matrix[j+3*i];
        }
    com[j]=com[j]/dim;
    }
}




void xyzprinter(float* matrix, int dim, FILE* file){
int i;

for(i=0;i<dim;i=i+1){
    fprintf(file, "%f %f %f \n", matrix[0+3*i],matrix[1+3*i],matrix[2+3*i]);
    }

fprintf(file,"\n");
}