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

/* Program created on 8 Nov 2007
   Modified and in working form on 13 Nov 2007

   This program will read in a text file from the Audi-Wapstra database
   (modified for the author's purposes) and offers the option to calculate
   Q-values for beta minus, beta plus, electron capture, and alpha decays. 
   To some extent it can calculate the supposed log ft value for the 
   transition but it should be cautioned that since they are the semi-
   empirical ones from Moszkowski et al, 1955, they should be checked
   against experimentally determined values in the primary literature
   or from the NNDC database in the ENSDF or XUNDL as desired.
   
   The calculations are done from mass defects, as outlined in Ehmann & Vance,
   or Krane. Q = mass defect(parent) - mass defect(daughter(s)) + other corrections
   as needed (e.g. for beta plus we need to add 1.022 MeV).
   
   The kinetic correction for detected alpha particles is also done.
   
   TODO: Add ability to "Fall through" to calculations using the SEMBE
   if one or both of the proposed isotopes do not have a mass defect available
   from Audi-Wapstra. */

// Time to declare function prototypes

int alphadecay(char* parent, char* daughter);
int betadecay(char* parent, char* daughter, int mode);
int discard_input(void);
int checknum(char *variable);
float find_element(char* element);

main()
{
      float def_parent = 0.0;     // Mass defect of the parent
      float def_daughter = 0.0;   // Mass defect of the daughter
      char parent[5];             // Parent isotope, in the form A(name)
      char daughter[5];           // Ditto the daughter, e.g. 99Tc
      char decay[2];              // Mode of decay
      char x[40];                 // Temporary line input variable
      char yesno;                 // Temp variable
      int status = 1;             // Whether to continue or not

      while (status == 1)
      {
       printf("Decay Q-value Calculator\n\n");
       printf("Enter the parent: ");
       scanf("%s", &parent);
       printf("Enter the daughter: ");
       scanf("%s", &daughter);

       printf("\nDesired mode of decay? (Enter a for alpha, b- for beta minus,\nb+ for beta plus, or ec for electron capture)\n-> ");
       scanf("%s", &decay);

       if (!strcmp(decay, "a")) alphadecay(parent, daughter);
       if (!strcmp(decay, "b-")) betadecay(parent, daughter, 1);
       if (!strcmp(decay, "b+")) betadecay(parent, daughter, 2);
       if (!strcmp(decay, "ec")) betadecay(parent, daughter, 3);
      
       discard_input();

       printf("\n\nWould you like another calculation (y/n) ? ");
       scanf("%c", &yesno);
       if (yesno == 'n') status = 0;
      }
}

/* Functions begin here */

int alphadecay(char* parent, char* daughter)
{
      const float def_alpha = 2424.911; // Mass defect of the alpha particle
      float qv = 0.0;             // Q-value result
      float ek = 0.0;             // Kinetic correction for the alpha energy
      float def_parent;           // Mass defect of the parent
      float def_daughter;         // Mass defect of the daughter
      int a_p = 0;                // A-value of the parent
      int a_d = 0;                // A-value of the daughter

      def_parent = find_element(parent);
      def_daughter = find_element(daughter);
      
      qv = (def_parent - (def_daughter + def_alpha)) / 1000.0;
      
      printf("\nAlpha decay Q-value: %.3f MeV", qv);

      a_p = checknum(parent);
      a_d = checknum(daughter);
      
      ek = ((float)a_d/(float)a_p) * qv;
      printf("\nOutgoing alpha particle kinetic energy: %.3f MeV", ek);
      
      return 0;
}

int betadecay(char* parent, char* daughter, int mode)
{
      float qv = 0.0;             // Q-value result
      float def_parent;           // Mass defect of the parent
      float def_daughter;         // Mass defect of the daughter
      int a_p = 0;                // A-value of the parent
      int a_d = 0;                // A-value of the daughter

      def_parent = find_element(parent);
      def_daughter = find_element(daughter);
      
      qv = (def_parent - def_daughter) / 1000.0;

      switch(mode) {
                  case 1:
                  printf("\nBeta minus decay");
                  break;
                  
                  case 2:
                  printf("\nBeta plus decay");
                  qv = qv - 1.022;
                  break;
                  
                  case 3:
                  printf("\nElectron capture");
                  break;
                        
                  default:
                  printf("Unspecified decay mode.\n");
                  abort();
                  }

      printf(" Q-value: %.3f MeV", qv);

      return 0;
}

int discard_input(void)
{
      int c;
      while ((c=getchar())!=EOF && c!='\n') ;
      return c==EOF;
} 

int checknum(char* variable)
{
      int i, x;
      char z[5];
      
      for (i = 0; i <= strlen(variable); i++) {
          x = isdigit(variable[i]);
          if (x == 0) { break; }
          }
          
      strncpy(z, variable, i);
      x = atoi(z);
      return x;
}

float find_element(char* element)
{
      FILE *fp;                   // Declares the file pointer
      float massdefect = 0.0;     // Mass defect of the element under consideration
      char x[40];                 // Comparison variable for the isotope to find
      int z;                      // Atomic number, currently discarded
      int status = 0;             // Allows us to return element not found

      fp = fopen("awm95s.dat", "r");
      if (fp == NULL) { 
             printf("Error opening file!");
             abort(); }

      while (!feof(fp))      
      {
            fscanf(fp, "%s %i %f", &x, &z, &massdefect);
            
            if (!strcmp(x, element)) { 
                            status = 1;
                            break; }
      }
      fclose(fp);
      
      /* TODO - fix this so if an element isn't found it can calculate the mass
      defect with a semi-empirical formula */
      
      if (status == 0) { printf("ERROR: Element not found!");
                      abort();
                      }

      return massdefect;
}