//
// layers.C
//
// Call 'cone' to divide a cluster into surface and several sub-surface layers.
//
// Written by: Yanting Wang          01/21/2003
//

#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <stdio.h>
#include <libgen.h>
#include <stdlib.h>

int main( int argc, char *argv[] )
{
  const char *NAME[10]={"H","He","Li","Be","B","C","N","O","F","Ne"};

  if( argc != 2)
  {
    cerr << "Usage: " << basename( argv[0] ) << " #_of_layers(<=10)" << endl;
    return 1;
  }

  const int nlayers = atoi( argv[1] );
  
  int num;
  char s[1024];

  cin >> num;
  cin.getline( s, 1024 );
  cin.getline( s, 1024 );

  //
  // Read in the input configuration
  //
  char atom[num][3];   // naming array
  double x[num], y[num], z[num];

  for( int i=0; i<num; i++ )
  {
    cin >> s >> x[i] >> y[i] >> z[i];
    strcpy( atom[i], "Au" );
  }

  //
  // Do surface discrimination recursively
  //
  int n = num;
  double xx[n], yy[n], zz[n];

  for( int i=0; i<num; i++ )
  {
    xx[i] = x[i];   yy[i] = y[i];   zz[i] = z[i];
  }

  for( int i=0; i<nlayers; i++ )
  {
    ofstream fo( "input.xyz" );

    // write "Au" atoms into "input.xyz"
    fo << n << endl << endl;
    for( int j=0; j<n; j++ )
    {
      fo << "Au " << xx[j] << " " << yy[j] << " " << zz[j] << endl;
    }
    fo.close();

    // find surface atoms
    system( "cone 5.0 0.5 < input.xyz > output.xyz" );

    ifstream fi( "output.xyz" );
    int n2;
    fi >> n2;
    fi.getline( s, 1024 );
    fi.getline( s, 1024 );

    n = 0;
    int ptr = 0;   // pointer to the naming array atom[]
    for( int j=0; j<n2; j++ )
    {
      double x, y, z;
      fi >> s >> x >> y >> z;

      while( strcmp( atom[ptr], "Au" ) ) ptr++;
      
      if( !strcmp( s, "Au" ) )
      {
	xx[n] = x;   yy[n] = y;   zz[n] = z;
	n++;
      }
      else
      {
	strcpy( atom[ptr], NAME[i] );
      }

      ptr++;
    }

    fi.close();

    if( n<1 )
    { 
      cerr << "Run out of atoms." << endl;
      break;
    }
  }

  system( "rm -f input.xyz" );
  system( "rm -f output.xyz" );

  cout << num << endl << endl;
  for( int i=0; i<num; i++ )
    cout << atom[i] << " " << x[i] << " " << y[i] << " " << z[i] << endl;

  return 0;
}