#include "pdb.h"

double real_margin = 8.0;
int chain_mapping[128] ;
int generation_index = 0;

file_records *the_file  ;
void initialize_chain_mapping ()
{
  char * arg = command_line_arg("chains") ;

  if (arg != NULL)
    {
      int len = strlen(arg) ;
      int ii ;
      
      for(ii=0; ii< len ; ii +=2)
	{
	  chain_mapping[arg[ii]] = arg[ii+1] ;
	  chain_mapping[arg[ii+1]] = arg[ii] ;
	}
    }
  else
    {
      int ii ;
      chain_mapping[' '] = '_' ;
      chain_mapping['_'] = ' ' ;

      for(ii = 'A' ; ii<= 'Z'; ii++)
	{
	  int inv = 'Z' - (ii -'A') ;
	  int low = ii - 'A' + 'a' ;
	  int invlow = 'z' - (low -'a') ;
	  chain_mapping[ii] = inv;
	  chain_mapping[inv] = ii;
	  chain_mapping[low] = invlow;
	  chain_mapping[invlow] = low;
	}
    }
}

#define PUT_IN_MINIMUM(x,max,min ) \
   if (x < min)                                \
    {                                          \
      while (x < min )                         \
	x += 1.0;                              \
      if (x > max ) return ;                   \
    }                                          \
  else                                         \
    {                                          \
      while (x > min )                         \
	x -= 1.0;                              \
	  x += 1.0 ;                           \
    }                                          
  


void atom_genearte_symmetry_1 (atom_record *atom, symmetry_operation *so, int index)
{
  double tx,ty,tz ;
  double tty,ttz ;
  atom_record tatom ;

  OPERATION_TRANSFORM(so,tx,ty,tz,atom->x,atom->y,atom->z) ;
  
  PUT_IN_MINIMUM(tx,maximum_x,minimum_x) ;
  PUT_IN_MINIMUM(ty,maximum_y,minimum_y) ;
  PUT_IN_MINIMUM(tz,maximum_z,minimum_z) ;

  bcopy(atom,&tatom,sizeof(atom_record)) ;

  while(tx < maximum_x)
    {
      for(tty = ty ; tty < maximum_y; tty += 1.0)
	{
	  for(ttz = tz ; ttz < maximum_z; ttz += 1.0)
	    {
	      tatom.x = tx;
	      tatom.y = tty;
	      tatom.z = ttz;
	      tatom.chain = chain_mapping[atom->chain] ;
	      tatom.number  = generation_index++;
	      un_fractionalize_atom(&tatom,the_file,0) ;
	      write_pdb_record(stdout,&tatom,BROOKHAVEN|IO_HYDROGEN) ;
	    }
	}
      tx += 1.0;
    }
}

static neibour_count = 0 ;
void *record_neibour_count(atom_record *atom, atom_record **neibours      ,int count)
{
  neibour_count = count ;
  return NULL;
}

void atom_genearte_symmetry (atom_record *atom, symmetry_operation *so, int index)
{
  double tx,ty,tz ;
  double tty,ttz ;
  atom_record tatom ;

  OPERATION_TRANSFORM(so,tx,ty,tz,atom->x,atom->y,atom->z) ;
  
  PUT_IN_MINIMUM(tx,maximum_x,minimum_x) ;
  PUT_IN_MINIMUM(ty,maximum_y,minimum_y) ;
  PUT_IN_MINIMUM(tz,maximum_z,minimum_z) ;

  bcopy(atom,&tatom,sizeof(atom_record)) ;

  while(tx < maximum_x)
    {
      for(tty = ty ; tty < maximum_y; tty += 1.0)
	{
	  for(ttz = tz ; ttz < maximum_z; ttz += 1.0)
	    {
	      if(fabs(ttz-atom->z) > 0.1 || 
		 fabs(tty-atom->y) > 0.1 || 
		 fabs(tx-atom->x) > 0.1 )
		{
		  tatom.x = tx;
		  tatom.y = tty;
		  tatom.z = ttz;

/*
void call_with_neibours(atom_record *atom,void (*function)() , int what );
*/
		  call_with_neibours(&tatom,record_neibour_count , CWN_ALL_ATOMS) ;
		  if(neibour_count)
		    {
		      tatom.chain = chain_mapping[atom->chain] ;
		      tatom.number  = generation_index++;
		      un_fractionalize_atom(&tatom,the_file,0) ;
		      write_pdb_record(stdout,&tatom,BROOKHAVEN) ;
		    }
		}
	    }
	}
      tx += 1.0;
    }
}

void generate_symmetry_atoms(file_records *file, symmetry_operation *so)
{
  traverse_file_atom_records(file, atom_genearte_symmetry,so) ;
}
int main(int argc,char ** argv)
{
  SAY_HI_FIRST ;
  CITATION ;
  do_initialization(argc,argv) ;
  initialize_chain_mapping() ;
  SET_DOUBLE_FROM_OPTION(real_margin, "margin") ;
  STDERR("Need to have the file symop.lib in the current directory or a link to it");
  {
    file_records *file = read_next_file_and_check(IO_HYDROGEN) ;
    space_group *sg = &file->sg ; 
    double smallest = (file->a > file->b) ? ((file->b > file->c) ?  file->c : file->b)
                                          :((file->a > file->c) ?  file->c : file->a);
    the_file = file ;
    distance_check = real_margin /smallest ;



    generation_index = file->atomnum +1 ;
    fractionalize(file) ;
    setup_neibouring_cube(file, distance_check) ;
    minimum_x -= real_margin / file->a ;
    maximum_x +=real_margin/file->a ;
    minimum_y -=  real_margin /file->b ;
    maximum_y += real_margin/file->b ;
    minimum_z -= real_margin/file->c ;
    maximum_z += real_margin/file->c ;


    if (! fill_space_group(sg, file->sg.s_name)) 
      serious_error(ARG_TYPE_STRING, "Failed tofind specae group",file->sg.s_name) ;
    LOOP_ALL_SG_OPERATIONS(sg ,
			    generate_symmetry_atoms(file,so) ;) ;
    
  }


}