/* *************************************************** */
/* Distribute for free!! and feel free to modify       */
/* as long as this header remains unchanged            */
/* Original coding and design                          */
/* Sxren Koch;skoch@geo.aau.dk                         */
/* *************************************************** */
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <stdarg.h>
#include "vector_math.c"

#define TRUE 1
#define FALSE 0

float dspacing(int hkl[],float celleparametre[]);
void datainput(float celleparametre[]);
float volumen(float celleparametre[]);
void base(float x[],float y[],float z[],float celleparametre[]);
void vector(float v[],int hkl[],float a[],float b[], float c[]);
void hklinput(int hkl[]);
void convert(float a[],int b[]);
int my_printf(char *str, ...);

FILE *debug_file;
char *debug_file_name = "out.txt";

FILE *input;
static void input_construct (void) __attribute__((constructor));
static void input_construct (void) { input = stdin;};


main()
{
  float celleparametre[6],vinkel;
  int hkl[3],zonea[3],fortsaet,cont,valg,max,antal_vektorer;
  char input_buffer;
  float d,vol,x[3],y[3],z[3],ar[3],br[3],cr[3],lar,lbr,lcr,var,vbr,vcr;
  float hkla[3],hklb[3];
  int hklra[3],hklrb[3],hklx[3],hkly[3];
  float angle,vektor1,vektor2,dspacex,dspacey,vektor1d,vektor2d,angled;
  int hkl1[3],hkl2[3],overflowkontrol;
  debug_file = fopen(debug_file_name, "w");
  printf("\n\n           Crystallography\n\n");
  printf("               Made by:\n    Soren Koch: skoch@geo.aau.dk\n");
  printf("\n\n\n  Input all distances in Aangstroem\n");
  printf("\n Output is generated in the file: OUT.TXT\n\n");
  datainput(celleparametre);
  base(x,y,z,celleparametre);
  my_printf("\n\n Cartesian kordinate-system placed with b parallel to y og a in the x-y-plane");
  my_printf("\n     Cartesian celleparameters:\n\n");
  my_printf("     a: %f  %f  %f\n     b: %f  %f  %f\n     c: %f  %f  %f\n",x[0],x[1],x[2],y[0],y[1],y[2],z[0],z[1],z[2]);
  vol = volumen(celleparametre);
  my_printf("\n       cellvolume: %f A",vol);
  crossp(ar,y,z);
  ar[0]=ar[0]/vol;
  ar[1]=ar[1]/vol;
  ar[2]=ar[2]/vol;
  lar = sqrt(ar[0]*ar[0] + ar[1]*ar[1] + ar[2]*ar[2]);
  crossp(br,z,x);
  br[0]=br[0]/vol;
  br[1]=br[1]/vol;
  br[2]=br[2]/vol;
  lbr = sqrt(br[0]*br[0] + br[1]*br[1] + br[2]*br[2]);
  crossp(cr,x,y);
  cr[0]=cr[0]/vol;
  cr[1]=cr[1]/vol;
  cr[2]=cr[2]/vol;
  lcr = sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
  my_printf("\n       resiprocal lattice-parameters - cartesian basis\n");
  my_printf("     a* %f  %f  %f\n",ar[0],ar[1],ar[2]);
  my_printf("     b* %f  %f  %f\n",br[0],br[1],br[2]);
  my_printf("     c* %f  %f  %f\n\n",cr[0],cr[1],cr[2]);
  var = vektorvinkel(br,cr);
  vbr = vektorvinkel(cr,ar);
  vcr = vektorvinkel(ar,br);
  my_printf("\n       Reciprocal lattice-parameters\n");
  my_printf("       a*: %f   b*: %f   c*: %f\n",lar,lbr,lcr);
  my_printf("   angles: %f      %f      %f\n",var,vbr,vcr);
  overflowkontrol = 0;
  for(;;)
  {
    overflowkontrol++;
    if (overflowkontrol == 100)
      {
	return;
      }	
    /*  ----------MENU-------- */
    my_printf("\n         Select subprograme:\n");
    printf("     1: Calculate d-spacings\n");
    printf("     2: Calculate angles in reciprocal space\n");
    printf("     3: Index TEM-negative\n");
    printf("     4: Calculate angles between hkl vectors\n");
    printf("     5: Calculate reciprocal vector form hkl-indices\n");
    printf("\n    99: EXIT\n\n Select: ");
    fflush(stdout);
    fscanf(input,"%i",&valg);
    /*  ------- END MENU ----- */
    if (valg == 99)
    {
      return;
    }
    if (valg == 1)
    {
      /* calculation of d-spacings from manually given hkl-indices */
      fortsaet = TRUE;
      while(fortsaet)
      {
        hklinput(hkl);
        d = dspacing(hkl,celleparametre);
        my_printf("\nhkl-vector: %i %i %i\n",hkl[0],hkl[1],hkl[2]);
        my_printf("d-spacing:  %f\n",d);
        vector(hkla,hkl,ar,br,cr);
        /* my_printf("\nReciprocal vector:  %f %f %f\n\n",hkla[0],hkla[1],hkla[2]); */
        printf("more d-spacing calculations (1 for Yes)?");
        fflush(stdout);
        fscanf(input, "%i",&cont);
	fortsaet = FALSE;
	if(cont == 1)
        {
          fortsaet = TRUE;
        }
	else
	{  
          fortsaet = FALSE;
	  valg = 0;
        }
	/* end d-spacing calculations */
      }
    }
    if (valg == 2)
    {
      /* CALCULATION OF ANGLES IN RECIPROCAL SPACE FROM 2 HKL-VECTORS */
      printf("\n Input first hkl:\n");
      hklinput(hklra);
      vector(hkla,hklra,ar,br,cr);
      printf("\n Input second hkl:\n");
      hklinput(hklrb);
      vector(hklb,hklrb,ar,br,cr);
      vinkel=vektorvinkel(hkla,hklb);
      my_printf("\n h,k,l: %i %i %i\n h,k,l: %i %i %i\n",hklra[0],hklra[1],hklra[2],hklrb[0],hklrb[1],hklrb[2]);
      my_printf(" Reciprocal angle: %f\n",vinkel);
      valg = 0;
      /* END RECIPROCAL CALCULATIONS */
    }
    if (valg == 3)
    {
      /* INDICING TEM-NEGATIVES USING 2 SHORTEST VECTORS AND INTERVECTORANGLE */
      max = 1;
      /* DATAINPUT BY USER */
      printf("\n\n  Input length of first d(longest) ");
      fscanf(input,"%f",&vektor1);
      printf("\n  Input estimated error of first d (in percent) ");
      fscanf(input,"%f",&vektor1d);
      printf("\n  Input length of second d ");
      fscanf(input,"%f",&vektor2);
      printf("\n  Input estimated error of second d (in percent) ");
      fscanf(input,"%f",&vektor2d);
      printf("\n  Input angle ");
      fscanf(input,"%f",&angle);
      if (angle>90)
      {
        angle=180-angle;
        printf("calculating for acute angle; %f",angle);
      }
      printf("\n  Input input estimated error of angle (in degrees) ");
      fscanf(input,"%f",&angled);
      printf("\n input maximal indices:");
      fscanf(input,"%i",&max);
      my_printf("\nCalculating for vectors\n");
      my_printf("d:     %f  error: %f percent\n",vektor1,vektor1d);
      my_printf("d:     %f  error: %f percent\n",vektor2,vektor2d);
      my_printf("angle: %f  error: %f degrees\n",angle,angled);
      antal_vektorer = 0;
      /* END DATAINPUT */
      /* GENERATING HKL-INDICES */
      for (hklx[0]=-max;hklx[0]<max+1;++hklx[0])
      {
        for (hklx[1]=-max;hklx[1]<max+1;++hklx[1])
        {
          for (hklx[2]=-max;hklx[2]<max+1;++hklx[2])
          {
            /* COMPARING FIRST D-SPACING */
            dspacex=dspacing(hklx,celleparametre);
            if (dspacex<(vektor1+vektor1*vektor1d/100))
            {
              if (dspacex>(vektor1-vektor1*vektor1d/100))
              {
                for (hkly[0]=-max;hkly[0]<max+1;++hkly[0])
                {
                  for (hkly[1]=-max;hkly[1]<max+1;++hkly[1])
                  {
                    for (hkly[2]=-max;hkly[2]<max+1;++hkly[2])
                    {
                      /* COMPARING 2. D-SPACING AND INTERVECTOR ANGLE */
                      dspacey=dspacing(hkly,celleparametre);
                      if (dspacey<(vektor2+vektor2*vektor2d/100))
                      {
                        if (dspacey>(vektor2-vektor2*vektor2d/100))
                        {
                          vector(hkla,hklx,ar,br,cr);
                          vector(hklb,hkly,ar,br,cr);
                          vinkel=vektorvinkel(hkla,hklb);
                          printf("",vinkel);
                          if (vinkel >90)
                          {
                            vinkel=180-vinkel;
                          }
                          if (vinkel < (angle+angled))
                          {
                            if (vinkel > (angle-angled))
                            {
                              /* OUTPUT FITTING VECTORS */
                              my_printf("vector fit:\n");
                              my_printf("h,k,l: %i %i %i  d: %f\n",hklx[0],hklx[1],hklx[2],dspacex);
                              my_printf("h,k,l: %i %i %i  d: %f\n",hkly[0],hkly[1],hkly[2],dspacey);
                              /* GENERATING ZONEAXIS FOR FITTING VECTORS*/
                              crosspi(zonea,hklx,hkly);
                              my_printf("angle: %f  Zoneaxis: %i %i %i\n",vinkel,zonea[0],zonea[1],zonea[2]);
                              antal_vektorer++;
                              /* END OUTPUT */
                            }
                          }
                        }
                      }     
                    }
                  }
                }
              }
            }
          }
        }
      }
      my_printf("found %i fitting vectors",antal_vektorer);
      /* END GENERATING INDICES AND END COMPARING */
      valg = 0;
    }
    if (valg == 4)
    {
      /* CALCULATION OF ANGLES BETWEEN 2 HKL-VECTORS */
      printf("\n Input first hkl:\n");
      hklinput(hkl1);
      vector(hkla,hkl1,x,y,z);
      printf("\n Input second hkl:\n");
      hklinput(hkl2);
      vector(hklb,hkl2,x,y,z);
      vinkel=vektorvinkel(hkla,hklb);
      my_printf("\n\n   Vectors:\n");
      my_printf("   h,k,l: %i %i %i\n",hkl1[0],hkl1[1],hkl1[2]);
      my_printf("   h,k,l: %i %i %i\n",hkl2[0],hkl2[1],hkl2[2]);
      my_printf("   Angle: %f\n",vinkel);
      valg = 0;
      /* END CALCULATIONS */
    }
    if (valg == 5)
    {
      /* CALCULATION OF RECIPROCAL VECTOR USING CARTESIAN BASIS */
      fortsaet = TRUE;
      while(fortsaet)
      {
        hklinput(hkl);
        vector(hkla,hkl,ar,br,cr);
        my_printf("\nhkl-vector:        %i %i %i\n",hkl[0],hkl[1],hkl[2]);
        my_printf("Reciprocal vector: %f %f %f\n\n",hkla[0],hkla[1],hkla[2]);
        printf("more calculations (1 for Yes)?");
        fflush(stdout);
        fscanf(input, "%i",&cont);
        fortsaet = FALSE;
        if(cont == 1)
        {
          fortsaet = TRUE;
        }
        else
        {
          fortsaet = FALSE;
          valg = 0;
        }
        /* end reciprocal calculations */
      }
      /* END CALCULATIONS */
      valg = 0;
    }
  }
  return;
}

float dspacing(int hkl[],float celleparametre[])
{
  /* GENERATING D-SPACING USING HKL-INDICES AND CELLPARAMETERS */
  /* USING TRICLINIC CELL */
  float d,e,f,g,h,i,j,k,l,pi;
  float a,b,c;
  float alfa,beta,gamma;
  pi = 3.141592;
  a = celleparametre[0];
  b = celleparametre[1];
  c = celleparametre[2];
  alfa = celleparametre[3]*2*pi/360;
  beta = celleparametre[4]*2*pi/360;
  gamma = celleparametre[5]*2*pi/360;
  e = power(hkl[0],2)*power(sin(alfa),2)/(a*a);
  f = power(hkl[1],2)*power(sin(beta),2)/(b*b);
  g = power(hkl[2],2)*power(sin(gamma),2)/(c*c);
  i = 2*hkl[0]*hkl[1]*(cos(alfa)*cos(beta) - cos(gamma))/(a*b);
  j = 2*hkl[1]*hkl[2]*(cos(beta)*cos(gamma) - cos(alfa))/(b*c);
  k = 2*hkl[2]*hkl[0]*(cos(gamma)*cos(alfa) - cos(beta))/(c*a);
  h = power(cos(alfa),2) + power(cos(beta),2) + power(cos(gamma),2) - 2*cos(alfa)*cos(beta)*cos(gamma);
  l = e + f + g + i + j + k;
  d = 1/sqrt(l/(1-h));
  /*  d = 1/sqrt(l/(1 - power(cos(alfa),2) - power(cos(beta),2) - power(cos(gamma),2) + 2*cos(alfa)*cos(beta)*cos(gamma)));*/
  return d;
}

void datainput(float celleparametre[])
{
  /* USERINPUT OF CELLPARAMETERS */
  /* FIRST 3 ELEMENTS IN ARRAY IS LENGTH */
  /* LAST 3 ELEMENTS IN ARRAY IS ANGLES */
  float a,b,c,alfa,beta,gamma;
  printf("\ninput cellparameters\n");
  printf("   a:   ");
  fscanf(input,"%f",&a);
  printf("   b:   ");
  fscanf(input,"%f",&b);
  printf("   c:   ");
  fscanf(input,"%f",&c);
  printf("  alfa: ");
  fscanf(input,"%f",&alfa);
  printf("  beta: ");
  fscanf(input,"%f",&beta);
  printf(" gamma: ");
  fscanf(input,"%f",&gamma);
  /* OUTPUT TO ARRAY */
  celleparametre[0]=a;
  celleparametre[1]=b;
  celleparametre[2]=c;
  celleparametre[3]=alfa;
  celleparametre[4]=beta;
  celleparametre[5]=gamma;
  return;
}

float volumen(float celleparametre[])
{
  /* CALCULATION OF CELLVOLUME USING CELLPARAMETERS */
  /* B-AXIS ALONG Y AND A-AXIS IN THE X-Y-PLANE */
  float volumen;
  float x[3],y[3],z[3],q[3];
  base(x,y,z,celleparametre);
  crossp(q,y,z);
  volumen = prikp(x,q);
  return volumen;
}

void base(float x[],float y[],float z[],float celleparametre[])
{
  /* CALCULATION AF CARTESIAN CELLPARAMETERS USING TRICLINIC CELLPARAMETERS */
  /* B-AXIS ALONG Y AND A-AXIS IN X-Y-PLANE */
  /* CARTESIAN VECTORS:X,Y,Z */
  float a,b,c,alfa,beta,gamma;
  float pi;
  pi=3.141592;
  a=celleparametre[0];
  b=celleparametre[1];
  c=celleparametre[2];
  alfa=celleparametre[3];
  beta=celleparametre[4];
  gamma=celleparametre[5];
  y[0] = 0;
  y[1] = b;
  y[2] = 0;
  x[0] = a*sin(gamma*2*pi/360);
  x[1] = a*cos(gamma*2*pi/360);
  x[2] = 0;
  z[1] = c*cos(alfa*2*pi/360);
  z[0] = (c*cos(beta*2*pi/360) - cos(gamma*2*pi/360)*c*cos(alfa*2*pi/360))/sin(gamma*2*pi/360);
  z[2] = sqrt(c*c - z[0]*z[0] - z[1]*z[1]);
  return;
}

void vector(float v[], int hkl[], float a[], float b[], float c[]) 
{
  /* CALCULATION OF A VECTORS FOR A GIVEN HKL-INDICES AND CARTESIAN BASIS */
  /* RETURNVECTOR:V, INDICES:hkl AND CARTESIAN CELLPARAMETERS:a,b,c */
  v[0] = hkl[0]*a[0] + hkl[1]*b[0] + hkl[2]*c[0];
  v[1] = hkl[0]*a[1] + hkl[1]*b[1] + hkl[2]*c[1];
  v[2] = hkl[0]*a[2] + hkl[1]*b[2] + hkl[2]*c[2];
  return;
}

void hklinput(int hkl[])
{
  /* USERINPUT OF HKL-INDICES */
  printf("\nInput hkl-indices\n");
  printf("  h? ");
  fscanf(input,"%i",&hkl[0]);
  printf("  k? ");
  fscanf(input,"%i",&hkl[1]);
  printf("  l? ");
  fscanf(input,"%i",&hkl[2]);
  return;
}

int my_printf(char *str, ...) 
{
  va_list vargs;
  int res;
  va_start(vargs, str);

  res = vprintf(str, vargs);
  va_end(vargs);
  if (debug_file) {
    va_start(vargs, str);
    res = vfprintf(debug_file, str, vargs);
    va_end(vargs);
  }
  return res;
}