/*************************************************************************
 * Programa para calcular las energias liberadas en reacciones nucleares *
 * segun la ecuacion de Weizsaecker.                                     *
 *                                                                       *
 * Compilar con: gcc -o weizsaecker weizsaecker.c -lm -lreadline         *
 *                                                                       *
 * Ecuaciones y valores tomadas de:                                      *
 *  "Nuclei & Particles", Emilio Segre, ISBN 84-291-4170-7               *
 *  "Kerphysik", Theo Mayer-Kuckuk, ISBN 3-519-13223-0                   *
 *                                                                       *
 * (C) 2005, Jaime Anguiano Olarra <jaime@gnome.org>                     *
 *************************************************************************/
#include <stdio.h>
#include <math.h>
#include <readline/readline.h>

#define elem_name(x) x->name
#define elem_energy(x) x->energie
#define elem_z(x) *(x->val)
#define elem_a(x) *((x->val) + 1)
#define me_energy 0.511

int DEBUG = 1;
int EXACT_ZA = 0;

const double u2mev = 931.494;
const double mn = 939.565;
const double mp = 938.272;
double me = me_energy;

const double a_V = 15.85;     /* Kondensationsenergie */
const double a_S = 18.34;     /* Oberflaechenenergie */
const double a_C = 0.71;      /* Coulomb-Energie */
const double a_A = 92.86;     /* Asymmetrie-Energie */
const double a_p_abs = 11.46; /* Paarungsenergie */

typedef struct element { int* val; const char* name; double energie; } elem;
elem* new_elem (int z, int a);
const char* elem_get_name (int z, int a);
void print_elem (elem* e);
double weizsaecker (int z, int a);
int run_interactive (void);
double exact_fission (int z, int a);

int
main (int argc, char** argv)
{
  int input[2];
  double w;
  double q;
  elem* e;
  
  input[0] = 0;
  input[1] = 0;
  q = 0;
  
  if (argc < 2) {
  __help:
    printf ("weizsaecker <Z> <A>\n"
	    "weizsaecker cmd\t\tModo interactivo\n"
	    "Para modo depuracion: entrar en modo interactivo, luego, entrar DEBUG\n"
	    "o !DEBUG segun se quiera activar o desactivar.\n\n"
	    "(C) Jaime Anguiano Olarra <jaime@gnome.org>\n\n");
    return -1;
  } else if (argc == 2){
    if (!(strcmp ("cmd", argv[1])))
      run_interactive ();
    else goto __help;
  } else if (argc == 3) {
    sscanf (argv[1], "%d", &input[0]);
    sscanf (argv[2], "%d", &input[1]);
    e = new_elem (input[0], input[1]);
    printf ("%s%d,%d: %e [~ %d MeV]\n", elem_name(e), elem_z(e),elem_a(e), elem_energy(e), (int) elem_energy(e));
    free (e);
  }
  return 0;
}

double 
weizsaecker (int z, int a)
{
  double m;
  double delta;

  /* Look for neutrons */
  if ((z == 0) && (a == 1)) {
    if (DEBUG)
      printf ("Neutron Energy: %e (MeV) [~ %d]\n", mn, (int) mn);
    return mn;
  }
  /* Look for protons */
  else if ((z == 1) && (a == 1)) {
    if (DEBUG)
      printf ("Proton Energy: %e (MeV) [~ %d]\n", mp, (int) mp);
    return mp;
  }
  else if ((z%2 == 0) && (a%2 == 0)) 
    delta = (-1) * a_p_abs * (1/sqrt(a));
  else if ((z%2 != 0) && (a%2 != 0)) 
    delta = a_p_abs * (1/sqrt(a));
  else
    delta = 0;
  if (DEBUG) {
    printf ("For element %d, %d:\n", z, a);
    if (EXACT_ZA)
      printf ("(a - z)*mn + z*(mp + me): %d MeV\n", 0);
    else
      printf ("(a - z)*mn + z*(mp + me): %e MeV\n", (a - z)*mn + z*(mp + me));
    printf ("- a_V * a: %e MeV\n", (-1) * a_V * a);
    printf ("+ a_S * pow(a, 2/3): %e (MeV)\n", a_S * pow(a, 2/3));
    printf ("a_A * pow(a/2 - z, 2)/a %e (MeV)\n", a_A * pow(a/2 - z, 2)/a);
    printf ("a_C * z*z/cbrt(a): %e MeV\n", a_C * z*z/cbrt(a));
    printf ("delta: %e MeV\n", delta);
  }
  if (EXACT_ZA)
    m = (double) 0;
  else
    m = (a - z)*mn + z*(mp + me);
  m -= a_V * a /* Kondensationsenergie */
    + a_S * pow(a, 2/3) /* Oberflaechenenergie */
    + a_A * pow(a/2 - z, 2)/a /* Asymmetrie-Energie */
    + a_C * z*z/cbrt(a) /* Coulomb-Energie */
    + delta; /* Paarungsenergie */
  if (DEBUG)
    printf ("Energy: %e [~ %d]\n", m, (int) m);
  return m;
}

elem*
new_elem (int z, int a)
{
  elem* e;
  int* vals;

  e = (elem*) malloc(sizeof(elem));
  vals = (int*) malloc (sizeof(int) * 2);
  vals[0] = z;
  vals[1] = a;
  e->val = vals;
  e->name = (const char*) elem_get_name(z, a);
  e->energie = weizsaecker (z, a);
  return e;
}

void
print_elem (elem* e)
{
  if (e->name != NULL)
    printf ("%s%d,%d: %e [~ %d (MeV)]", elem_name(e), elem_z(e), elem_a(e), elem_energy(e), (int) elem_energy(e));
  else printf ("Z: %d  A: %d %e [~ %d (MeV)]", elem_z(e), elem_a(e), elem_energy(e), (int) elem_energy(e));
}

const char*
elem_get_name (int z, int a)
{
  switch (z) {
  case 0:
    return "n";
  case 1:
    return "H";
  case 2:
    return "He";
  case 82:
    return "Pb";
  case 86:
    return "Rn";
  case 88:
    return "Ra";
  case 90:
    return "Th";
  case 92:
    return "U";
  case 94:
    return "Pu";
  default:
    return NULL;
  }
}

double
berech_energie (elem** elem_array, int length)
{
  double enrg;
  int ctr;
  elem* m;

  enrg = 0;
  for (ctr = 0; ctr < length; ctr++) {
    m = *(elem_array + ctr);
    enrg += fabs(elem_energy(m));
    printf ("Nuclid [%d] %e MeV\n", ctr, elem_energy(m));
  }
  return enrg;
}

int
run_interactive (void)
{
  char* line;
  int reacts_ctr;
  int prods_ctr;
  int z, a, rs, ps, i;
  double energie, energie_ps;
  
  elem** reactives;
  elem** products;

  printf ("<<-- Weizsaecker 1.0 -->\n\nMDW Rechenmaschine\n\n");
 __interact_init:
  reacts_ctr = 0;
  prods_ctr = 0;
  
  reactives = (elem**) malloc (sizeof(elem*) * 4);
  products = (elem**) malloc (sizeof(elem*) * 6);
  
  while (reacts_ctr < 4) {   
    line = readline ("Enter reactive <Z> <A> (MAX 4, Enter to end): ");
    if (!strcmp("DEBUG", line)) {
      DEBUG = 1;
      printf ("Debugging mode activated...\n");
      continue;
    } else if (!strcmp("!DEBUG", line)) {
      DEBUG = 0;
      continue;
    } else if ((!strcmp("EXACT_ZA", line)) || (!strcmp("e", line))) {
      EXACT_ZA = 1;
      printf ("Exact reaction expected...\n");
      continue;
    } else if ((!strcmp("!EXACT_ZA", line)) || (!strcmp("!e", line))) {
      EXACT_ZA = 0;
      continue;
    } else if (!strcmp("!electrons", line)) {
      me = 0;
      continue;
    } else if (!strcmp("electrons", line)) {
      me = me_energy;
      continue;
    } else if (!strcmp("fission", line)) {
      line = readline ("Enter fissionable nuclid <Z> <A>: ");
      if (sscanf(line, "%d %d", &z, &a))
	printf ("Energy for exact fission: %e (MeV) [~ %d]\n", exact_fission (z, a), (int) exact_fission(z,a));
      goto __interact_end;
    }
 
    if (sscanf (line, "%d %d", &z, &a) > 0)
      reactives[reacts_ctr++] = new_elem (z, a);
    else {
      rs = reacts_ctr;
      reacts_ctr = 4;
    }
  }
  while (prods_ctr < 6) {
    line = readline ("Enter product <Z> <A> (MAX 6, Enter to end): ");
    if (sscanf (line, "%d %d", &z, &a) > 0)
      products[prods_ctr++] = new_elem (z, a);
    else {
      ps = prods_ctr;
      prods_ctr = 6;
    }
  }
  printf ("Reactives energy (RE): %e (MeV) [~ %d]\n", energie = berech_energie (reactives, rs), (int) energie);
  printf ("Products energy (PE): %e (MeV/c^2) [~ %d]\n", energie_ps = berech_energie (products, ps), (int) energie_ps);
  printf ("RE - PE: %e (MeV/c^2) [~ %d]\n", energie - energie_ps, (int) (energie - energie_ps));
  for (i = 0; i < rs; i++) 
    free (reactives[i]);
  for (i = 0; i < ps; i++)
    free (products[i]);
  free (reactives);
  free (products);
 __interact_end:
  if (!(strcmp("y", line = readline("Another reaction? (y/N): ")))) goto __interact_init;
} 

double
exact_fission (int z, int a)
{
  return ((-5.12 * pow (a, 2/3) + 0.284 * z * z / cbrt(a)) * u2mev);
}