source: proiecte/NBody/PlummerModel/ModelParallel.c @ 164

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>

#include "Model.h"

body* bodies; //vector that contains all the bodies in the system

void setBody(body* b, double mass, double pos[3], double vel[3]){
                
        b->mass = mass;
        memcpy(b->pos, pos, 3 * sizeof(double));
        memcpy(b->vel, vel, 3 * sizeof(double));
}


void initializeNBody(int n){
        int i;
        double mass = 0;
        double pos[3] = {0, 0, 0};
        double vel[3] = {0, 0, 0};
        
        #pragma omp parallel for default(none) private(i, mass, pos, vel) shared(n, bodies)
        for (i = 0; i < n; i++){
                setBody(&bodies[i], mass, pos, vel);
                //printf("%f %f\n", bodies[i].mass, bodies[i].pos[2]);
        }
        #pragma omp barrier
}

//generate random float number between [low, high)
double frand(double low, double high)
{
        double temp;

        /* swap low & high around if the user makes no sense */
        if (low > high){
                temp = low;
                low = high;
                high = temp;
        }

        /* calculate the random number & return it */
        temp = (rand() / ( RAND_MAX + 1.0))
                * (high - low) + low;
        return temp;
}



void mkplummer(int n, unsigned int seed){
        int i;
        double mass;
        double pos[3];
        double vel[3];
        double radius, theta, phi;
        double x, y;
        double velocity;

        if(seed == 0)
                srand(time(NULL)); //different succession of results in the subsequent calls to rand
        else srand(seed);   //used for debugging

        //initializeNBody(n);
        
        #pragma omp parallel for default(none)\
                private(i, mass,radius,theta, phi, pos, vel, x, y, velocity) \
                shared(n, bodies)
        for(i = 0;i < n; i++){
                mass = 0;
                pos[0] = 0; pos[1] = 0; pos[2] = 0;
                vel[0] = 0; vel[1] = 0; vel[2] = 0;
                setBody(&bodies[i], mass, pos, vel);
                mass = 1.0 / n;
                radius = 1.0 / sqrt(pow(frand(0, 1),-(2.0 / 3.0)) - 1.0);
                theta = acos(frand(-1, 1));
                phi = frand(0, 2 * PI);
                pos[0] = radius * sin(theta) * cos(phi);
                pos[1] = radius * sin(theta) * sin(phi);
                pos[2] = radius * cos(theta);
                x = 0.0;
                y = 0.1;
                while( y > x * x * pow (1.0 - x*x, 3.5)){
                        x = frand(0, 1);
                        y = frand(0, 0.1);
                }
                velocity = x * sqrt(2.0) * pow(1.0 + radius * radius, -0.25);
                theta = acos(frand(-1, 1));
                phi = frand(0, 2 * PI);
                vel[0] = velocity * sin(theta) * cos(phi);
                vel[1] = velocity * sin(theta) * sin(phi);
                vel[2] = velocity * cos(theta);
                setBody(&bodies[i], mass, pos, vel);
        }

        #pragma omp barrier

}

void writeData(int n){
        char line[256];
        int i;  
        FILE *f = fopen(FILE_NAME, "w+");
        if(f == NULL){
                printf("Error opening file\n");
                return;
        }
        sprintf(line, "%d\n", n);
        fputs(line, f);
        sprintf(line, "%d\n", 0);
        fputs(line, f);
        for(i = 0; i < n; i++){
                //printf("here\n");
                sprintf(line, "%f %f %f %f %f %f %f\n", 
                        bodies[i].mass,
                        bodies[i].pos[0], bodies[i].pos[1], bodies[i].pos[2],
                        bodies[i].vel[0], bodies[i].vel[1], bodies[i].vel[2]);
                fputs(line, f);
        }
        fclose(f);

}

int main(int argv, char** argc){
        int n;
        unsigned int seed;
        if(argv != 3){
                printf("usage ./Model nr_bodies seed\n");
                return 0;
        }

        n = atoi(argc[1]);
        seed = atoi(argc[2]);

        bodies = (body*) malloc(n * sizeof(body));
        mkplummer(n, seed);
        /*int id;
        #pragma omp master{
                id = omp_get_thread_num();
                printf("thread %d\n", id);
                writeData(n);
        }*/
        //writeData(n);
        return 1;
}