#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpc.h>
#define DELTA (0.5)

typedef struct {
    double len; 
    double wid; 
    double hei; 
    double mass;
} rail;

nettype HeteroNet(int n, double v[n]) {
  coord I=n;
  node {I>=0: v[I];};
  parent [0];
};

double Density(double x, double y, double z) {
  return 6.0*sqrt( exp( sin( sqrt( x*y*z ) ) ) );
}

double RailMass(double len, double wid, double hei, double delta) {
  double mass, x, y, z;
  for(mass=0., x=0.; x<len; x+=delta)
    for(y=0.; y<wid; y+=delta)
      for(z=0.; z<hei; z+=delta)
        mass += Density(x,y,z);
  return mass*delta*delta*delta;
}

int [*]main(int [host]argc, char **[host]argv) {
  repl N=3;

  if(argc>1)
    N = [host]atoi(argv[1]);

  if(N>0) {
    static rail [host]steel_hedgehog[[host]N];
    repl double volumes[N], [host]start;
    int [host]i;
    repl j;

    for(i=0; i<[host]N; i++) {
      steel_hedgehog[i].len = 200.*(i+1);
      steel_hedgehog[i].wid = 5.*(i+1);
      steel_hedgehog[i].hei = 10.*(i+1);
    }

    start = [host]MPC_Wtime();

    for(j=0; j<N; j++)
      volumes[j] =
          steel_hedgehog[j].len * steel_hedgehog[j].wid * steel_hedgehog[j].hei;

    recon RailMass(20., 4., 5., 0.5);
    {
      net HeteroNet(N, volumes) mynet;
      [mynet]: {
        rail myrail;

        myrail = steel_hedgehog[];
        myrail.mass =
                   RailMass(myrail.len, myrail.wid, myrail.hei, DELTA);
        MPC_Printf("Rail #%d is %gcm x %gcm x%gcm and weights %g kg\n",
                   I coordof mynet, myrail.len, myrail.wid,
                   myrail.hei, myrail.mass/1000.0);
        [host]printf("The steel hedgehog weights %g kg\n",
               [host]((myrail.mass)[+])/1000.0);
      }
    }
    [host]printf("\nIt took %.1f seconds to run the program.\n",
                 [host]MPC_Wtime()-start);
  }
  else
    [host]printf("Wrong input (N=%d)\n", [host]N);
}