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

#define MAXRADIUS   25. // massimo raggio dell'ammasso ammesso
#define STARMIN   1000
#define STARMAX 100000
#define SIZE       100
//#define _VARIANTE_1
//#define _VARIANTE_2
#define _VARIANTE_3

float getNumber(char *msg, float min, float max);
float z();  
float theta();
void fillDistribution(int xdistrib[], int n, float x);
void fillRadiusDistribution(int xdistrib[], int n, float r, float x);

int main() {
  int i;
  int R[SIZE] = {0};
  float radius = getNumber("Inserisci il raggio medio dell'ammasso: ", 0., MAXRADIUS);
  int N = 0;
  float n = 0.;
  do {
    n = getNumber("# Inserisci il numero di stelle: ", STARMIN, STARMAX);
    if (n != floor(n)) {
      printf("#    Errore: il numero dev'essere intero.\n");
      n = 0.;
    }
  } while (n <= 0.);
  N = n;
  char fname[255];
  printf("# Inserisci il nome del file di dati: ");
  scanf("%s", fname);
  FILE *fs = fopen(fname, "w");
  for (i = 0; i < N; i++) {
    float r = radius*z();
    float t = theta();
    float x = r*cos(t);
    float y = r*sin(t);
#ifdef _VARIANTE_3
    x = radius*z();
    y = radius*z();
    r = sqrt(x*x + y*y);
#endif
    fprintf(fs, "%f %f\n", x, y);
#ifdef _VARIANTE_1
    fillDistribution(R, SIZE, x);
#endif
#ifdef _VARIANTE_2
    fillRadiusDistribution(R, SIZE, r, x);
#endif
#ifdef _VARIANTE_3
    fillRadiusDistribution(R, SIZE, r, x);
#endif
  }
  fclose(fs);
  char dfname[255] = "";
  sprintf(dfname, "distrib%s", fname);
  fs = fopen(dfname, "w");
  for (int i = 0; i < SIZE; i++) {
    fprintf(fs, "%d %d\n", i-SIZE/2, R[i]);
  }
  fclose(fs);
}

void fillRadiusDistribution(int xdistrib[], int n, float r, float x) {
  int index = r;
  if (x < 0.) {
    index = n/2 - r;
  } else {
    index += n/2;
  }
  if ((index > 0) && (index < n)) {
    xdistrib[index]++;
  }
}

void fillDistribution(int xdistrib[], int n, float x) {
  int index = x;
  if (x < 0.) {
    index = n/2 + x;
  } else {
    index += n/2;
  }
  if ((index > 0) && (index < n)) {
    xdistrib[index]++;
  }
}

float z() {
  float u1 = (float)lrand48()/RAND_MAX;
  float u2 = (float)lrand48()/RAND_MAX;
  return sqrt(-2.*log(u1))*cos(2.*M_PI*u2);
}

float theta() {
  return 2.*M_PI*lrand48()/RAND_MAX;
}

float getNumber(char *msg, float min, float max) {
  float ret;
  int n = -1;
  do {
    printf("# %s", msg);
    n = scanf("%f", &ret);
    if ((ret < min) || (max > max)) {
      printf("#    Errore: inserisci %f < x < %f\n", min, max);
      n = -1;
    }
    if (n != 1) {
      char buffer[255];
      scanf("%s", buffer);
    }
  } while (n != 1);
  return ret;
}