/* selenomet.c++  Routines related to using selenomet atom positons.
  Need to use space group operators, etc.
*/

#include "drawwindow.h"
#include "math.h"
#include "trace.h"
#include "selenomet.h"

#define TRUE		1
#define FALSE		0
#define LINE_WIDTH	100	/* Maximum width of line read in	*/

/* ***************************************************************** */
/* neg_power -- Takes x ^ n, with n a negative number.               */
/* ***************************************************************** */
double neg_power(int x, int n)
{
        double ans;
        n = -n;
        for (ans = 1; n > 0; --n)
                ans = ans * x;
        return(1. / ans);
}

/* ***************************************************************** */
/* read_number -- Takes a line and the pointer to a starting point   *
 *              and reads a number up until a non - numeric digit.   *
 *              Returns the number.  *pos is the position of the     *
 *              last digit.                                          */
/* ***************************************************************** */
float read_number(char op_line[LINE_WIDTH], int *pos, int sign)
{
        float num;
        int power;
/*
        printf("Line col = %d, sign = %d\n",*pos, sign);
*/
        for(num = 0; op_line[*pos] >='0' && op_line[*pos] <='9'; (*pos)++)
                num = 10 * num + op_line[*pos] - '0';
        if (op_line[*pos] == '.') {
                (*pos)++;
                for(power=-1;op_line[*pos]>='0' && op_line[*pos]<='9';(*pos)++){
                        num = num + neg_power(10, power)*(op_line[*pos] - '0');
                        power--;
                }
        }
/*
        printf ("returning %f, line col = %d\n", sign*num, *pos);
*/
        /* Decrement *pos (as it will be incremented by calling function */
        (*pos)--;
        return (sign * num);
}

/* ***************************************************************** */
/* make_operators -- Takes an array of characters corresponding to   *
 *              the symmetry operators, parses the input, and forms  *
 *              the symmetry and translation matrices in direct      *
 *              space.                                               */
/* ***************************************************************** */
void make_operators(int sym[NSYMOPS][3][3], float tranp[NSYMOPS][3],
                                        char op_line[LINE_WIDTH], int op_num)
{
        char digit;
        int line_col, sign;
        float num, denom;
        int row = 0;
        sign = 1;
        num = 0; denom = 1;

        printf("Symmetry operator number %d: %s", op_num + 1, op_line);

        /* Parse symmetry line and generate matrix      */
        for (line_col = 0; line_col < LINE_WIDTH; line_col++) {

                digit = op_line[line_col];
                if (digit == NULL || digit == '\n') {
                        break;
                }
                else if(digit == ',') {
                        if (!denom) {
                                printf("\nError! Division by zero\nExiting.\n");                                exit(1);
                        }

                        tranp[op_num][row] = num / denom;
                        sign = 1;
                        num = 0; denom = 1;
                        row++;
                        continue;
                }

                /* Set sign if a + or - is found */
                else if(digit == '+' || digit == '-')
                        sign = (digit == '+') ? 1 : -1;

                else if(digit == ' ' || digit == '/')
                        continue;

                /* Read in x, y, or z   */
                else if(digit == 'x' || digit == 'X') {
                        sym[op_num][row][0] = sign * 1;
                        sign = +1;
                } else if(digit == 'y' || digit == 'Y') {
                        sym[op_num][row][1] = sign * 1;
                        sign = +1;
                } else if(digit == 'z' || digit == 'Z') {
                        sym[op_num][row][2] = sign * 1;
                        sign = +1;
                } else if((digit >= '0' && digit <= '9') || digit == '.') {
                        if (!num) {
                                num = read_number(op_line, &line_col, sign);
                                sign = +1;
                        } else if (denom == 1) {
                                denom = read_number(op_line, &line_col, sign);
                                sign = +1;
                        }

                } else
                        printf("\nUnexpected character: %c\n", digit);
        }
        if (row < 2) {
                printf("\nWarning:  Incomplete symmetry generator\n");
        }

        if (!denom) {
                printf("\nError!  Division by zero\nExiting.\n");
                exit(1);
        }

        tranp[op_num][row] = num / denom;

}


/* ***************************************************************** */
/* fget_operators -- Asks for a space group, locates that group in   *
 *		reference file "symop.lib", and reads in the symm    *
 *		operators in direct space.  Returns # of symm. ops.  */
/* ***************************************************************** */
int fget_operators(int intnum,int sym[NSYMOPS][3][3], float tranp[NSYMOPS][3])
{
	FILE *fp;
	char op_line[LINE_WIDTH];
	int found_flag, group_no, op_num, zero;

	found_flag = FALSE;
	op_num = 0;

#if 0
	printf("\nInput the number of the space group ");
	printf("that you would like to use: ");
	scanf("%d", &group_no);
	spacegroup = group_no;
	printf("\n");
#endif
        group_no= intnum;

	if ((fp = fopen("symop.lib","r")) == NULL) {
		printf("Cannot open symop.lib data file.\nExit code 1\n");
		exit(1);
	}
	while(fgets(op_line, LINE_WIDTH, fp) != NULL) {
		
		/* This is just a variable to make read_number happy!	*/
		zero = 0;

		if (read_number(op_line, &zero, 1) == group_no) {
#if 0
			system("clear");
#endif
			printf("Space group %d found.\n", group_no);
			printf("\nData on spacegroup %d:\n", group_no);
			printf("Number, # symm ops, # primitive ");
			printf("symm ops, Spacegroup name, Pointgroup name\n");
			printf("%s\n", op_line);
			found_flag = TRUE;
			break;
		}
	}

	if (found_flag == FALSE) {
		printf("Could not find space group %d!\n", group_no);
		exit(1);
	}

	while(fgets(op_line, LINE_WIDTH, fp) != NULL && op_line[0] == ' ') {
		make_operators(sym, tranp, op_line, op_num);
		op_num++;
	}

	return(op_num);
}

#if 0

void readselenomet(char selenometfile[],int *intnum,int *nsel,float selpos[MAXSEL][3])
/* Read the selenomet file
  Assumes that positions are in fractional format
*/
{
   int i;
   char name[200];
   FILE *fp;

   if( (fp=fopen(selenometfile,"r")) == NULL){
          printf("Selenomet file %s does note exit- setting nsel = 0\n",selenometfile);
          *nsel = 0;
   }
   else{ 
      	fgets(name,200,fp); /* read space group number */
	sscanf(name,"%d",intnum);
      	fgets(name,200,fp); /* read number selenomtg */
	sscanf(name,"%d",nsel);
        if(*nsel>MAXSEL-2){
           printf("ERROR - more than %d selenomet sites - truncating\n",MAXSEL);
           *nsel = MAXSEL-2;
        } 
        for(i=1;i<= *nsel;i++){
      	  fgets(name,200,fp); /* read number of title lines */
	  sscanf(name,"%f %f %f",&selpos[i][0],&selpos[i][1],&selpos[i][2]);
        }
   }
   fclose(fp);
}

void symop(int isym,int sym[NSYMOPS][3][3],float ftranp[NSYMOPS][3],
        float pos1[3],float pos2[3])
{
   int i;

   multvec(sym[isym],pos1,pos2); /* returns vector v = m*v1*/
   for(i=0;i<=1;i++)
       pos2[i] = pos2[i]+ftranp[isym][i];
}


void findselmap(struct griddata& den1,int intnum,int *nsel,float selpos[MAXSEL][3])
/* Finds all the selenomet postions that lie within the map region
  1) Apply all sym operators to all selenomet positions, and then
      and then keep increasing and decreasing fraction by 1 until pass
      through map.
  intnum is number of the space group
*/
{
   int i,j,k,nsym,tempnsel,npos[10][3],isel,sym[NSYMOPS][3][3];
   float ftranp[NSYMOPS][3],tempselpos[MAXSEL][3];
   float frmin[3],frmax[3],xx[10][3];

   orthtofr(den1,den1.ximin,frmin); // convert real space mapmin to fractional
   orthtofr(den1,den1.ximax,frmax); // convert real space mapmin to fractional
   for(i=0;i<NSYMOPS;i++)/* set all operators = 0 */
      for(j=0;j<=2;j++){
         ftranp[i][j]=0;
         for(k=0;k<=2;k++){
            sym[i][k][j]=0;
         }
      }          
    nsym=fget_operators(intnum,sym,ftranp);
    tempnsel =0; //counts number of selomet in map region
    for(isym=1;isym<=nsym;isym++){
     for(isel=1;isel<= *nsel;isel++){
       symop(isym,sym,ftranp,selpos[isel],pos2);
       for(k=0;k<=2;k++){ // find x,y,z positons in map
          npos[k]=0; // counts # of k pos in map
          if(pos2[k]<frmax[k]) // make all pos2 > frmax
              pos2[k] = pos2[k] + (int)(frmax[k]-pos2[k]) +1; // round up
       }
       for(k=0;k<=2;k++){ // find x,y,z positons in map
          while(pos2[k]>frmax[k]){ // find first pos < frmax
             pos2[k]= pos2[k]-1;
          }
          npos[k] = npos[k]+1;
          if(npos[k]>9){
             printf("ERROR: map seems to be too large - npos = %d\n",npos[k]);
             goto quitrout;
          }
          xx[ npos[k] ][k] = pos2[k];
          do{
             pos2[k] = pos2[k]-1;
             if(pos2[k]>=frmin[k])
               npos[k] = npos[k]+1;
               if(npos[k]>9){
                 printf("ERROR: map seems to be too large - npos = %d\n",npos[k]);
                 goto quitrout;
               }
               xx[ npos[k] ][k] = pos2[k];
             }
          while(pos2[k]>=frmin[k]);
        }
        for(i=1;i<=npos[0];i++)
           for(j=1;j<npos[1];j++)
              for(k=1;k<npos[2];k++){
                 tempnsel= tempnsel+1;
                 tempselpos[ tempnsel ][0] = xx[i][0];  
                 tempselpos[ tempnsel ][1] = xx[j][1];  
                 tempselpos[ tempnsel ][2] = xx[k][2]; 
                 if(tempnsel>MAXSEL-2){
                     printf("ERROR finselmap - more than %d selenomet sites - truncating\n",MAXSEL)
                     goto findpos;
                 }
              }
      }//isel 
    } // isym 
    findpos:;
    *nsel = tempnsel;
    for(i=1;i<=tempnsel;i++){
       for(k=0;k<=2;k++)
         selpos[i][k] = tempselpos[i][k];
    }                  
}

void findselpos(struct griddata& den1,char selenometfile[],int *intnum,
      int *nsel,float selpos[MAXSEL][3])
/* Reads selenomet file, 
  the applys sym operator to all selenomet positions, and translates
    into map region
  intnum is spacegroup number
*/
{
    readselenomet(selenometfile,intnum,nsel,selpos);
    findselmap(den1,*intnum,nsel,selpos);
}
#endif
