#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "mucal.h"
#include "panel.h"

float CromerWaber(int ZZ, float SL)
{
/*
#UD  Parametrization of f0 (the non-dispersive part of the atomic 
#UD  scattering factor) vs sin(theta/lambda).
#UD  f0, the non-dispersive part of the atomic scattering factor, is a
#UD  function of the selected element and of sin(theta)/lambda, where
#UD  lambda is the photon wavelengh and theta is incident angle.
#UD  This function can be approximated by a function:
#UD
#UD    f0[SL] = c + [SUM a_i*EXP(-b_i*(SL^2)) ]
#UD                i=1,4
#UD
#UD  where SL = sin(theta) / lambda and c, a_i and b_i
#UD  are the so called Cromer-Mann coefficients (see International
#UD  Tables vol. 4 or vol C; in vol. C refer to pg 500-502) and
#UD  tabulated in this file (in columns) :  a1 a2 a3 a4 c b1 b2 b3 b4
*/
   FILE *in;
   static int first=1, Znew[100];
   static float Cromer[100][9];
   char iline[200], *lineptr, lastchar;
   int i, len, Z;
   float c, a1, a2, a3, a4, b1, b2, b3, b4, SL2, form;

   if (first) {
      in = fopen("f0_CromerMann.dat", "rt");
      if (in == NULL) {
         printf("Could not open 'f0_CromerMann.dat'\n");
         exit(0);
      }

next:
      fgets(iline, 200, in);
      if (feof(in)) goto got_it;
        
      lineptr = strstr(iline, "#S");
      if (lineptr == NULL) goto next;
      len = strlen(iline);
      lastchar = iline[len-2];
      lineptr++;
      lineptr++;
      sscanf(lineptr, "%d", &Z);
      if ((lastchar == ' ') && (Znew[Z] == 0)) { 
         iline[len-1] = 0;
/*         printf("Z: %2d, lastchar %1c,  %s  ", Z, lastchar, iline);
*/
         Znew[Z] = Z;
         fgets(iline, 200, in);
         if (feof(in)) goto got_it;
         fgets(iline, 200, in);
         if (feof(in)) goto got_it;
         fgets(iline, 200, in);
         if (feof(in)) goto got_it;
/*         printf("%s", iline);
*/
         sscanf(iline, "%f %f %f %f %f %f %f %f %f",
            &Cromer[Z][0], &Cromer[Z][1], &Cromer[Z][2], &Cromer[Z][3], &Cromer[Z][4], 
            &Cromer[Z][5], &Cromer[Z][6], &Cromer[Z][7], &Cromer[Z][8]);
/*         printf("                                ");   
         printf("%8.5f   %8.5f   %8.5f   %8.5f   %8.5f   %8.5f   %8.5f   %8.5f   %8.5f\n",  
            Cromer[Z][0], Cromer[Z][1], Cromer[Z][2], Cromer[Z][3], Cromer[Z][4], 
            Cromer[Z][5], Cromer[Z][6], Cromer[Z][7], Cromer[Z][8]);
*/
      }
      goto next;
   	
got_it:
      fclose(in);
   }
   
   a1 = Cromer[ZZ][0];
   a2 = Cromer[ZZ][1];
   a3 = Cromer[ZZ][2];
   a4 = Cromer[ZZ][3];
   c = Cromer[ZZ][4];
   b1 = Cromer[ZZ][5];
   b2 = Cromer[ZZ][6];
   b3 = Cromer[ZZ][7];
   b4 = Cromer[ZZ][8];
   
   SL2 = SL * SL;
   form = c + a1 * exp(-b1*SL2) + a2 * exp(-b2*SL2) 
            + a3 * exp(-b3*SL2) + a4 * exp(-b4*SL2);

   return(form);
}