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

/*

      program estat

      compute electrostatic potential energy for N atomic nuclei in 3D space


*/

#define NBATOMS   1000
#define NTH       4


int     bounds[NTH+1];
double  pos[NBATOMS][3];


void    make_bounds (int [], int, int);
void    make_pos    (double [][3], int);
double  esfunc      (double [][3], int, int, int []);


int main (int argc, char* argv[])
{
    double  energy;
    int     i;
    int nbcomp;

    if (argc != 2) {
        fprintf(stderr, "Syntax: %s <nbcomp>\n", argv[0]);
        exit(1);
    }
    nbcomp = atoi(argv[1]);

    /* build bound array to evenly slice the job: */
    make_bounds (bounds, NBATOMS, NTH);

    /* build position array: */
    make_pos (pos, NBATOMS);

    /* do it several time for ease of time measurement: */
    for (i=1; i<= nbcomp; i++) {
        energy = esfunc (pos, NBATOMS, NTH, bounds);
    }

    printf (" energy = %f\n", energy);

    return 0;
}


void make_bounds (int bounds[], int nbatoms, int nth)
{
    double  x, rn;
    int     i;

    bounds[0]   = 0;
    bounds[nth] = nbatoms-1;

    x  = 0.0;
    rn = 1.0 / (double) nth;

    for (i=1; i < nth; i++) {
        x = sqrt(rn + x*x);
        bounds[i] = (int) (x * nbatoms);
    }

    for (i=0; i <= nth; i++) {
        printf (" bounds[%d] = %d\n", i, bounds[i]);
    }
}


void make_pos (double pos[][3], int nbatoms)
{
    double  rr3, t;
    int     i;

    rr3 = 1.0 / sqrt(3.0);

    for (i=0; i < nbatoms; i++) {
        t = rr3 * (double) i;

        pos[i][0] = t;
        pos[i][1] = t;
        pos[i][2] = t;
    }
}


double esfunc (double pos[][3], int nbatoms, int nth, int bounds[])
{
    double  sum, sum_arr[64];
    double  partsum, e, e1, e2;
    int     l1, l2, is, i, j;

    for (is=0; is < nth; is++) {
        l1 = bounds[is]+1;
        l2 = bounds[is+1];

        partsum = 0.0;

        for (i=l1; i<=l2; i++) {
            for (j=0; j < i; j++) {
                double  dx, dy, dz;

                dx = pos[i][0] - pos[j][0];
                dy = pos[i][1] - pos[j][1];
                dz = pos[i][2] - pos[j][2];

                e1 = dx*dx + dy*dy + dz*dz;
                e1 = 1.0 / e1;
                partsum += sqrt(e1);

            }
        }

        /* printf ("es: is=%d  l1=%d  l2=%d  partsum=%f\n", is, l1, l2, partsum); */
        sum_arr[is] = partsum;
    }

    sum = 0.0;
    for (is=0; is < nth; is++) {
        sum += sum_arr[is];
    }

    return sum;
}