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

/* Program created on 8 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(int mode);
int discard_input(void);
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("Desired 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(1);
       if (!strcmp(decay, "b+")) betadecay(2);
       if (!strcmp(decay, "ec")) betadecay(3);
      
       discard_input();

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

/* Functions begin here */

int alphadecay(char* parent, char* daughter)
{
      float qv = 0.0;             // Q-value result
      const def_alpha = 2424.911; // Mass defect of the alpha particle
      float def_parent;           // Mass defect of the parent
      float def_daughter;         // Mass defect of the daughter

      def_parent = find_element(parent);
      def_daughter = find_element(daughter);
      
      qv = def_parent - (def_daughter + def_alpha);
      
      printf("\nAlpha decay Q-value: %d MeV", qv);
      return 0;
}

int betadecay(int mode)
{
      float def_parent = 0.0;     // Mass defect of the parent
      float def_daughter = 0.0;   // Mass defect of the daughter
      printf("\nNothing to see here for now, move along.\n");
      switch(mode) {
                  case 1:
                  printf("Beta minus decay.\n");
                  break;
                  
                  case 2:
                  printf("Beta plus decay.\n");
                  break;
                  
                  case 3:
                  printf("Electron capture.\n");
                  break;
                        
                  default:
                  printf("Unspecified decay mode.\n");
                  break;
                  }
      return 0;
}

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

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);
      
      if (status == 1) {
                 printf("Element found!");
                 }
                 else {
                      printf("ERROR: Not found!");
                      abort();
                      }

      return massdefect;
}