#include "transmission.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>


void usage()
{
  printf("\
Neutron cross section calculation.\n\
Usage:\n\
======\n\n\
Input parameters:\n\
-----------------\n\
 -w, --wavelength  lambda           neutron wavelength in Angstrom\n\
                   begin:step:end   neutron wavelength range in Angstrom\n\
 -e, --energy_meV  value            neutron energy in meV\n\
                   begin:step:end   neutron energy range in meV\n\n\
Unit cell input parameters:\n\
---------------------------\n\
 -f, --parameter_file  filename     nxs parameter file with all unit cell\n\
                                    information needed\n\
\n\
 In case any of the following arguments is used the information in the nxs\n\
 parameter file will be omitted. The argument values will be used instead.\n\
\n\
 -s, --space_group      number      the space group number\n\
 -a, --lattice_a        angstrom    lattice constant a in Angstrom\n\
 -b, --lattice_b        angstrom    lattice constant b in Angstrom\n\
 -c, --lattice_c        angstrom    lattice constant c in Angstrom\n\
 -A, --lattice_alpha    degree      lattice angle alpha in degrees\n\
 -B, --lattice_beta     degree      lattice angle alpha in degrees\n\
 -s, --lattice_gamma    degree      lattice angle alpha in degrees\n\
 -D, --debye_temp       kelvin      Debye temperature in Kelvin\n\
\n\n\
Output parameters:\n\
------------------\n\
 -o, --output [ARG1, ARG2, ...]\n\
         xsect_cohel       outputs neutron coherent elastic scattering\n\
                           cross section\n\
         xsect_cohinel     outputs neutron coherent inelastic scattering\n\
                           cross section\n\
         xsect_inel        outputs neutron total inelastic scattering\n\
                           cross section using a combination on of single and\n\
                           multi phonon calculations\n\
         xsect_inelBinder  outputs neutron total inelastic scattering\n\
                           cross section after BINDER (1970)\n\
         xsect_sph         outputs the single phonon contribition\n\
         xsect_mph         outputs the multi phonon contribition as proposed\n\
                           by ADIB (2007)\n\
         xsect_mphCassels  outputs the multi phonon contribition as proposed\n\
                           by CASSELS (1950)\n\
         xsect_mphFreund   outputs the multi phonon contribition as proposed\n\
                           by FREUND (1983)\n\
         xsect_incel       outputs neutron incoherent elastic scattering\n\
                           cross section\n\
         xsect_incinel     outputs neutron incoherent inelastic\n\
                           contribution after BINDER (1970)\n\
         xsect_abs         outputs neutron absorption cross section\n\
         xsect_total       outputs neutron total cross section as xsect_cohel\n\
                           + xsect_abs + xsect_inel + xsect_incel\n\
         transmission      outputs transmisson with the help of --length_cm\n\
         attenuation       outputs attenuation coefficient\n\
\n\
 If output is not set the default output option is 'attenuation'.\n\
\n\n\
Special options:\n\
----------------\n\
 -t, --temperature       kelvin     sample temperature in Kelvin,\n\
                                    default is 293.0 K\n\
 -d, --density           g/cm^3     sample density in g/cm^3\n\
 -m, --max_hkl_index     integer    defines the maximum hkl index\n\
 -c2, --mph_c2           value      a fitting parameter for the multi phonon\n\
                                    calculation; by default mph_c2 is\n\
                                    calculated after FREUND (1983), but can be\n\
                                    set differently here\n\
 -l, --length_cm         cm         transmission length through the sample in\n\
                                    cm, default is 1.0 cm\n\
\n\
Examples:\n\
---------\n");
  printf("   nxsCalc -parameter_file Fe.nxs -wavelength 3.5\n\n");
  printf("   nxsCalc -parameter_file Fe.nxs -wavelength 1.0:0.01:3.5 \n\
      -lattice_a 2.865 -debye_temp 435.0\n\n");
}





int main( int argc, char* argv[] )
{
  int i=0;
  int j=0;

  double *wavelength = (double*)malloc( sizeof(double) );
  int numWavelengths=0;
  short unsigned int energy = 0;
  char *tokptr = NULL;

  NXS_UnitCell uc;
  NXS_AtomInfo *atomInfoList = NULL;
  int numAtomInfos = 0;
  int maxHKL_index = 8;
  short unsigned int verbose = 0;
  double unitcell_temperature = 293.0;
  double lattice_a=0.0, lattice_b=0.0, lattice_c=0.0;
  double lattice_alpha=0.0, lattice_beta=0.0, lattice_gamma=0.0;
  double length_cm = 1.0;

  enum NXSCalc
  {
    NXS_Abs, NXS_Cohel, NXS_Cohinel, NXS_Inel, NXS_InelBinder, NXS_SPh, NXS_MPh,
    NXS_MPhCassels, NXS_MPhFreund, NXS_Incel, NXS_Incinel, NXS_Total, NXS_Atten,
    NXS_Trans
  };
  int *output = (int*)malloc( sizeof(enum NXSCalc) );
  int numOutputs = 1;

  double unitcell_density = -1.0;
  double mph_c2 = -1.0;

  if( argc < 2 )
  {
    usage();
    return -1;
  }


  /* parsing the arguments */
  for( i=1; i<argc; ++i )
  {
    if( argv[i][0]!='-' )
      continue;

    if( strcmp(argv[i], "--wavelength")==0 || strcmp(argv[i],"-w")==0 )
    {
      double start = 0.0;
      double step = 0.0;
      double stop = 0.0;
      tokptr = NULL;
      tokptr = strtok(argv[++i], ":");
      if( tokptr!=NULL )
      {
        start = atof( tokptr );
        tokptr = strtok(NULL, ":");
        if( tokptr!=NULL )
        {
          step = atof( tokptr );
          tokptr = strtok(NULL, ":");
          if( tokptr!=NULL )
          {
            stop = atof( tokptr );
          }
        }

        if( start>1E-6 && step>1E-6 && stop>1E-6)
        {
          numWavelengths = (int)(fabs(stop-start)/step)+1;
          wavelength = (double*)realloc( wavelength, numWavelengths*sizeof(double) );
          if( start>stop )
            step = -step;
          for( j=0; j<numWavelengths; j++ )
            wavelength[j] = step*j+start;
        }
        else if ( start>1E-6 )
        {
          wavelength = (double*)realloc( wavelength, sizeof(double) );
          wavelength[0] = start;
          numWavelengths = 1;
        }
        else
        {
          fprintf(stderr, "Error: -wavelength option not correct!\n" );
          return -1;
        }
        energy = 0;
      }
    }
    else if( strcmp(argv[i],"--energy_meV")==0 || strcmp(argv[i],"-e")==0 )
    {
      double start = 0.0;
      double step = 0.0;
      double stop = 0.0;
      tokptr = NULL;
      tokptr = strtok(argv[++i], ":");
      if( tokptr!=NULL )
      {
        start = atof(tokptr);
        tokptr = strtok(NULL, ":");
        if( tokptr!=NULL )
        {
          step = atof(tokptr);
          tokptr = strtok(NULL, ":");
          if( tokptr!=NULL )
          {
            stop = atof(tokptr);
          }
        }
        if( start>1E-6 && step>1E-6 && stop>1E-6)
        {
          numWavelengths = (int)(fabs(stop-start)/step)+1;
          wavelength = (double*)realloc( wavelength, numWavelengths*sizeof(double) );
          if( start>stop )
            step = -step;
          for( j=0; j<numWavelengths; j++ )
            wavelength[j] = 9.04457036579 / sqrt(step*j+start);
        }
        else if ( start>1E-6 )
        {
          wavelength = (double*)realloc( wavelength, sizeof(double) );
          wavelength[0] = 9.04457036579 / sqrt( start );
          numWavelengths = 1;
        }
        else
        {
          fprintf(stderr, "Error: -energy option not correct!\n" );
          return -1;
        }
        energy = 1;
      }
    }
    else if( strcmp(argv[i],"--parameter_file")==0 || strcmp(argv[i],"-f")==0 )
    {
      numAtomInfos = nxs_readParameterFile( argv[++i], &uc, &atomInfoList );
    }
    else if( strcmp(argv[i],"--space_group")==0 || strcmp(argv[i],"-s")==0 )
    {
      printf("%s\n",argv[++i]);
    }
    else if( strcmp(argv[i],"--lattice_a")==0 || strcmp(argv[i],"-a")==0 )
    {
      lattice_a = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--lattice_b")==0 || strcmp(argv[i],"-b")==0 )
    {
      lattice_b = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--lattice_c")==0 || strcmp(argv[i],"-c")==0 )
    {
      lattice_c = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--lattice_alpha")==0 || strcmp(argv[i],"-A")==0 )
    {
      lattice_alpha = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--lattice_beta")==0 || strcmp(argv[i],"-B")==0 )
    {
      lattice_beta = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--lattice_gamma")==0 || strcmp(argv[i],"-C")==0 )
    {
      lattice_gamma = atof( argv[++i] );
    }
    else if( strcmp(argv[i],"--debye_temp")==0 || strcmp(argv[i],"-D")==0 )
    {
      printf("%s\n",argv[++i]);
    }
    else if( strcmp(argv[i],"--temperature")==0 || strcmp(argv[i],"-t")==0 )
    {
      unitcell_temperature = atof( argv[++i] );
      if( unitcell_temperature<1E-6 )
        unitcell_temperature = 293.0;
    }
    else if( strcmp(argv[i],"--output")==0 || strcmp(argv[i],"-o")==0 )
    {
      tokptr = NULL;
      tokptr = strtok(argv[++i], ",");
      numOutputs = 0;
      while( tokptr!=NULL )
      {
        output = (int*)realloc( output, (++numOutputs)*sizeof(enum NXSCalc) );
        if( strcmp(tokptr,"attenuation")==0 ) output[numOutputs-1] = NXS_Atten;
        else if( strcmp(tokptr,"xsect_cohel")==0 ) output[numOutputs-1] = NXS_Cohel;
        else if( strcmp(tokptr,"xsect_cohinel")==0 ) output[numOutputs-1] = NXS_Cohinel;
        else if( strcmp(tokptr,"xsect_inel")==0 ) output[numOutputs-1] = NXS_Inel;
        else if( strcmp(tokptr,"xsect_inelBinder")==0 ) output[numOutputs-1] = NXS_InelBinder;
        else if( strcmp(tokptr,"xsect_sph")==0 ) output[numOutputs-1] = NXS_SPh;
        else if( strcmp(tokptr,"xsect_mph")==0 ) output[numOutputs-1] = NXS_MPh;
        else if( strcmp(tokptr,"xsect_mphCassels")==0 ) output[numOutputs-1] = NXS_MPhCassels;
        else if( strcmp(tokptr,"xsect_mphFreund")==0 ) output[numOutputs-1] = NXS_MPhFreund;
        else if( strcmp(tokptr,"xsect_incel")==0 ) output[numOutputs-1] = NXS_Incel;
        else if( strcmp(tokptr,"xsect_incinel")==0 ) output[numOutputs-1] = NXS_Incinel;
        else if( strcmp(tokptr,"xsect_abs")==0 ) output[numOutputs-1] = NXS_Abs;
        else if( strcmp(tokptr,"xsect_total")==0 ) output[numOutputs-1] = NXS_Total;
        else if( strcmp(tokptr,"transmission")==0 ) output[numOutputs-1] = NXS_Trans;
        else
        {
          fprintf(stderr, "Error: wrong output argument(s)!\n" );
          return -1;
        }
        tokptr = strtok(NULL, ",");
      }
    }
    else if( strcmp(argv[i],"--verbose")==0 || strcmp(argv[i],"-v")==0 )
    {
      verbose = 1;
    }
    else if( strcmp(argv[i],"--help")==0 || strcmp(argv[i],"-h")==0 )
    {
      usage();
      return 1;
    }
    else if( strcmp(argv[i],"--density")==0 || strcmp(argv[i],"-d")==0 )
    {
      unitcell_density = atof( argv[++i] );
      if( unitcell_density < 1E-6 )
      {
        fprintf(stderr, "Warning: wrong argument for --density, using auto calculation instead!\n" );
        unitcell_density = 0.0;
      }
    }
    else if( strcmp(argv[i],"--max_hkl_index")==0 || strcmp(argv[i],"-m")==0 )
    {
      maxHKL_index = atoi( argv[++i] );
      if( maxHKL_index < 1 )
      {
        fprintf(stderr, "Warning: wrong argument for --max_hkl_index, using --max_hkl_index = 8 instead!\n" );
        maxHKL_index = 8;
      }
    }
    else if( strcmp(argv[i],"--mph_c2")==0 || strcmp(argv[i],"-c2")==0 )
    {
      mph_c2 = atof( argv[++i] );
      if( mph_c2 < 1E-6 )
      {
        fprintf(stderr, "Warning: wrong argument for --mph_c2, using auto calculation instead!\n" );
        mph_c2 = 0.0;
      }
    }
    else if( strcmp(argv[i],"--length_cm")==0 || strcmp(argv[i],"-l")==0 )
    {
      length_cm = atof( argv[++i] );
      if( length_cm < 1E-6 )
      {
        fprintf(stderr, "Warning: wrong argument for --length_cm, using 1.0cm instead!\n" );
        length_cm = 1.0;
      }
    }
  }

  if( numAtomInfos>=0 && verbose )
  {
    printf("NXS unit cell definition is:\n"
          "space group number = %s\n"
          "a = %.5f \t\t alpha = %.5f\n"
          "b = %.5f \t\t beta  = %.5f\n"
          "c = %.5f \t\t gamma = %.5f\n"
          "# label  b_coherent  sigma_inc  sigma_abs  molar_mass  x  y  z\n",
          uc.spaceGroup, uc.a, uc.alpha, uc.b, uc.beta, uc.c, uc.gamma);
  }

  if( NXS_ERROR_OK != nxs_initUnitCell(&uc) )
  {
    fprintf(stderr, "ERRROR: nxs_initUnitCell() failed\n");
    return -1;
  }

  for(j=0; j<numAtomInfos; j++)
  {
    nxs_addAtomInfo( &uc, atomInfoList[j] );
    if( verbose )
      printf("%s  %f  %f  %f  %f  %f  %f  %f\n",
         atomInfoList[j].label, atomInfoList[j].b_coherent, atomInfoList[j].sigmaIncoherent,
         atomInfoList[j].sigmaAbsorption, atomInfoList[j].molarMass,
         atomInfoList[j].x[0], atomInfoList[j].y[0], atomInfoList[j].z[0] );
  }

  if( lattice_a>1E-6 ) uc.a = lattice_a;
  if( lattice_b>1E-6 ) uc.b = lattice_b;
  if( lattice_c>1E-6 ) uc.c = lattice_c;
  if( lattice_alpha>1E-6 ) uc.alpha = lattice_alpha;
  if( lattice_beta>1E-6 ) uc.beta = lattice_beta;
  if( lattice_gamma>1E-6 ) uc.gamma = lattice_gamma;

  uc.temperature = unitcell_temperature;

  if( unitcell_density>1E-6 )
    uc.density = unitcell_density;

  uc.maxHKL_index = maxHKL_index;

  if( mph_c2>1E-6 )
    uc.mph_c2 = mph_c2;

  nxs_initHKL( &uc );

  if( verbose )
  {
    printf( "Debye temp. = %f K\n", uc.debyeTemp );
    printf( "Temperature = %f K\n", uc.temperature );
    printf( "Unit cell density = %f g/cm^3\n", uc.density );
    printf( "Unit cell volume = %f A^3\n", uc.volume );
    printf( "Unit cell mass = %f\n", uc.mass );
    printf( "Max. hkl = %i\n", uc.maxHKL_index );
    printf( "C2 (mult. phonon) = %f\n", uc.mph_c2 );
    printf( "Transmission length = %f cm\n", length_cm );
    printf( "\n" );
  }

  for( i=0; i<numWavelengths; i++ )
  {
    if( energy )
      printf("%f ", 81.80425310173 / wavelength[i] / wavelength[i] );
    else
      printf("%f ", wavelength[i] );

    for( j=0; j<numOutputs; j++)
    {
      switch (output[j])
      {
        case NXS_Cohel:        printf("%f ", nxs_CoherentElastic(wavelength[i],&uc) ); break;
        case NXS_Abs:          printf("%f ", nxs_Absorption(wavelength[i],&uc) ); break;
        case NXS_Cohinel:      printf("%f ", nxs_CoherentInelastic(wavelength[i],&uc) ); break;
        case NXS_Inel:         printf("%f ", nxs_TotalInelastic(wavelength[i],&uc) ); break;
        case NXS_InelBinder:   printf("%f ", nxs_TotalInelastic_BINDER(wavelength[i],&uc) ); break;
        case NXS_SPh:          printf("%f ", nxs_SinglePhonon(wavelength[i],&uc) ); break;
        case NXS_MPh:          printf("%f ", nxs_MultiPhonon(wavelength[i],&uc) ); break;
        case NXS_MPhCassels:   printf("%f ", nxs_MultiPhonon_CASSELS(wavelength[i],&uc) ); break;
        case NXS_MPhFreund:    printf("%f ", nxs_MultiPhonon_FREUND(wavelength[i],&uc) ); break;
        case NXS_Incel:        printf("%f ", nxs_IncoherentElastic(wavelength[i],&uc) ); break;
        case NXS_Incinel:      printf("%f ", nxs_IncoherentInelastic(wavelength[i],&uc) ); break;
        case NXS_Total:        printf("%f ", totalcrossection(wavelength[i],&uc) ); break;
        case NXS_Trans:        printf("%f ", transmission(wavelength[i],&uc, length_cm) ); break;
        case NXS_Atten:
        default:               printf("%f ", attenuation(wavelength[i],&uc) ); break;
      }
    }
    printf("\n");
  }

  return 0;
}