Integración de Monte Carlo en MPI

Sea 4 dimensiones: x, y, z, t. Hay un toro de tamaño variable en el tiempo:

(R(t) - sqrt(x²+y²) )² +z² <= r(t)²

 
R(t): distancia del centro del tubo del toro al centro del toro (0,0,0)
r(t): radio del tubo

Dentro del toro se genera energia, con una potencia de P(x,y,z,t) Watts (=Joules/s)
Diga cuanta energia se genera desde el tiempo 0 al tiempo 10 para los siguientes R,r,P:

R(t)=100-t; r(t)=(t+1)²/2; P(x,y,z,t)=abs(z)  

#include <stdio.h >
#include <stdlib.h >
#include <time.h >
#include <math.h >
#include "mpi.h"

double R(double t){
  return 100.0 - 1.0*t;
}

double r(double t){
  return ((t+1)*(t+1))/2.0;
}

double P(double x, double y, double z, double t){
  return (double)abs(z);
}

int main(int argc, char *argv[])
{
  int nprocs, myproc;
  int i, n;
  double acum, sum_acum;
  double E, w, x, y, z, t, xx, yy, R_prima, r_prima, W;
  double t_i, t_f;
  double t0, t1;
  double x_max, x_min;
  double y_max, y_min;
  double z_max, z_min;
  double avg;

  /* Init MPI */
  MPI_Init( &argc, &argv );
  MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
  MPI_Comm_rank( MPI_COMM_WORLD, &myproc );

  if( argc == 2 ) {      
    sscanf( argv[1], "%d", &n );
  }
  else {
    if (myproc == 0)
    {
      printf("Input the number of random points: (0 quit) ");fflush(stdout);
      scanf("%d",&n);     
    }
 }

  if( myproc == 0 ) printf("Calculating E using %d points\n", n);

  MPI_Bcast( &n, 2, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Barrier( MPI_COMM_WORLD );
  t0 = MPI_Wtime();

  E = 0.0;
  sum_acum = 0;
  acum = 0;
  t_i = 0.0;
  t_f = 10.0;
  x_min = y_min = -150.5; //-(R(t) + r(t))
  x_max = y_max = 150.5;  //R(t) + r(t)
  z_min = -60.5;  // -r(t)
  z_max = 60.5;   // r(t)
  srand( (unsigned)time( NULL )+myproc );
  
  for( i=myproc; i<n; i+=nprocs ) {
    t = ((double) rand()/RAND_MAX)*(t_f - t_i);    
    x = (((double) rand()/RAND_MAX)* (x_max - x_min)) - x_max;
    y = (((double) rand()/RAND_MAX)* (y_max - y_min)) - y_max;
    z = (((double) rand()/RAND_MAX)* (z_max - z_min)) - z_max;
    r_prima = r(t);
    R_prima = R(t);
    xx = ((double) R_prima - (double)sqrt(x*x + y*y));
    if(xx*xx + z*z <= r_prima*r_prima)
      acum += P(x, y, z, t); 
  }
  MPI_Allreduce( &acum, &sum_acum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
  avg = (double)sum_acum / (double)n; 
  
  E = (t_f - t_i) * (x_max - x_min) * (y_max - y_min) * (z_max - z_min) * avg;

  t1 = MPI_Wtime();
  if( myproc == 0 ) {
    printf("E =  %f Joules\n", E);
    printf("Time = %f s.\n", t1-t0);
  }
  MPI_Finalize();
  return 0;
}