#!/usr/bin/perl -w
####################################################################
# feff_tables.pl: generate tables of central atom phase shifts
#                 and scattering amplitudes and phase shifts
#     copyright (c) 2002, Bruce Ravel <ravel@phys.washington.edu>
#     Time-stamp: <2002/01/29 14:43:00 bruce>
#
# This script does a big loop in which it generates a feff.inp file
# and runs feff for each atom and for each edge from K to M5.  The
# central atoms phase shift column from the feff0001.dat file is read
# and saved in two forms: flat ascii, ifeffit pad, and perl
# Storable.  Edges below 100eV are not computed.
#
# Based on or developed using Distribution: GENFMT 2.0
# GENFMT 2.0 Copyright (c) [2002] University of Washington
####################################################################

use strict;
use Chemistry::Elements qw(get_name get_symbol);
use Xray::Absorption;
Xray::Absorption -> load("elam");
use Ifeffit qw(ifeffit put_array);
use File::Path;

use Getopt::Long;

## =================================================================
##    configuration variables
my $feffcmd   = 'feff828';	# how feff is run
my $feffid    = 'Feff 8.20x8';	# how feff identifies itself

my $capsout   = 'caps.dat';	# ascii data file
my $capsgroup = 'caps';		# ifeffit group name for CAPS
my $capsdir   = 'caps';		# subdir for PAD files

my $scatout   = 'scat.dat';	# ascii data file
my $scatgroup = 'scat';		# ifeffit group name for scattering
my $scatdir   = 'scat';		# subdir for PAD files

#my $storout   = 'caps.db';	# Storable data file
## =================================================================

my $version = 0.01;
my $date = "January 25, 2002";

my ($do_ascii, $do_pad) = (1,1);
&GetOptions("ascii!"  => \$do_ascii,
	    "pad!"    => \$do_pad,
	    "h"       => \&usage,
	    "help"    => \&usage,
	    "v"       => \&version,
	    "version" => \&version,
	   );

tee();

sub tee {
  my $pid;
  return if $pid = open(STDOUT, "| tee feff_tables.log");
  die "cannot fork: $!" unless defined $pid;
  $| = 1;
  while (<STDIN>) { print };
};


unless ($do_pad or $do_ascii) {
 print "You turned off all output types!  feff_tables.pl quitting.\n";
 exit;
};





## =================================================================
## Snarf up input file template
my @feffinp;
while (<DATA>) {
  push @feffinp, $_
};



## =================================================================
## write output file header
if ($do_ascii) {
  open CAPS, ">$capsout" or die "could not open $capsout for writing\n";
  open SCAT, ">$scatout" or die "could not open $scatout for writing\n";

  open H, "headers/top" or die "could not open headers/top for reading\n";
  while (<H>) {  print CAPS '# ', $_; print SCAT '# ', $_ }
  close H;

  foreach (@feffinp) { print CAPS '# ', $_; print SCAT '# ', $_ }

  open H, "headers/license" or die "could not open headers/license for reading\n";
  while (<H>) {  print CAPS '# ', $_; print SCAT '# ', $_ }
  close H;
};


## =================================================================
## Big loop over edges and elements
my $first_time = 1;
my $ncaps = 0;
my @k = ();
if ($do_pad) {
  (-d $capsdir) or mkpath($capsdir, 1);
  (-d $scatdir) or mkpath($scatdir, 1);
};
foreach my $z (1 .. 94) {
  my %m	    = ();
  my %p	    = ();
  my $found = 0;
  my ($name, $sym) = (get_name($z), get_symbol($z));
  foreach my $edge (qw(K L1 L2 L3 M1 M2 M3 M4 M5)) {
    ## this is an arbitrary cutoff...
    next if (Xray::Absorption->get_energy($z, $edge) < 100);
    ++$found;

    my @mag   = ();
    my @phase = ();

    ## fix up the template for this calculation
    open FI, ">feff.inp" or die "could not open feff.inp for writing\n";
    foreach my $x (@feffinp) {
      (my $y = $x) =~ s/\%n/$name/;
      $y =~ s/\%s/$sym/;
      $y =~ s/\%z/$z/;
      $y =~ s/\%e/$edge/;
      print FI $y;
    };
    close FI;

    ## run feff
    print $/;
    system $feffcmd;

    ## parse the feff0001.dat file and save columns 2, 3, 4 to each of
    ## the output formats
    my @this = ();
    open FD, "feff0001.dat" or die "could not open feff0001.dat for reading\n";
    my $switch = 0;
    my $i = 0;
    foreach (<FD>) {
      ($switch = 1), next if (/\@\#$/);
      next unless $switch;
      my @line = split;
      ($first_time) and push @k, sprintf("%11.1f", $line[0]);
      #my $phase = - ( $line[1] + 5*$k[$i] );
      my $phase = $line[1];
      ++$i;
      push @this,  sprintf("%11.5f", $phase);
      push @mag,   sprintf("%11.5f", $line[2]);
      push @phase, sprintf("%11.5f", $line[3]);
    };
    close FD;
    if ($first_time and $do_ascii) {
      print CAPS "wavenumber", @k, $/;
      print SCAT "wavenumber\n     ", @k, $/;
    }
    if ($do_ascii) {
      print CAPS sprintf("%2s %2s     ", $sym, $edge), @this, $/;
      print SCAT sprintf("%2s %2s\n     mag   ", $sym, $edge), @mag,
	"\n     phase ", @phase, $/;
    };
    $do_pad and put_array("$capsgroup.${sym}_$edge", \@this);
    $m{$edge} = \@mag;		# store scattering arrays
    $p{$edge} = \@phase;
    $first_time = 0;
  };

  next unless $found;

  if ($do_pad) {
    ++$ncaps;
    ## save the ifeffit PAD files for this atom
    ## first the CAPS arrays
    put_array("$capsgroup.k", \@k);
    my $zz = sprintf("%2.2d", $z);
    ifeffit("save(file=$capsdir/$zz.dat, no_strings, no_scalars, no_sys, with_arrays)\n");
    ifeffit("erase \@group $capsgroup\n");

    ## then the scattering amplitudes and phase shifts
    foreach my $edge (qw(K L1 L2 L3 M1 M2 M3 M4 M5)) {
      next unless exists $m{$edge};
      put_array("$scatgroup.${sym}_${edge}_mag", $m{$edge});
      put_array("$scatgroup.${sym}_${edge}_pha", $p{$edge});
    };
    put_array("$scatgroup.k", \@k);
    ifeffit("save(file=$scatdir/$zz.dat, no_strings, no_scalars, no_sys, with_arrays)\n");
    ifeffit("erase \@group $scatgroup\n");
  };
};

## save the ascii files and the Storable file
$do_ascii and close CAPS;
$do_ascii and close SCAT;

print "\n\n\n";
($do_ascii) and
  print "Wrote $capsout and $scatout.\n";
($do_pad) and
  print "Wrote $ncaps files to $capsdir/ and $scatdir/.\n";
print "\n";

## finally, clean up feff's spew
foreach (qw(atoms.dat axafs.dat chi.dat feff.bin feff.inp feff0001.dat
	    files.dat fms.bin fort.3 fpf0.dat geom.dat global.dat
	    ldos.inp list.dat log.dat log1.dat log2.dat log3.dat
	    log4.dat log5.dat log6.dat logdos.dat misc.dat mod1.inp
	    mod2.inp mod3.inp mod4.inp mod5.inp mod6.inp paths.dat
	    phase.bin pot.bin xmu.dat xsect.bin)) {
  next unless (-e $_);
  unlink $_;
};


sub version {
  print <<EOH

feff_tables.pl $version ($date)
Bruce ravel <ravel\@phys.washington.edu>
http://feff.phys.washington.edu/software/feff_tables/

EOH
  ;
  exit;
};


sub usage {
  print <<EOH;

Generate tables of central atom phase shifts and of scattering
amplitudes and phase shifts by running feff8 on a token local
structure.  By default the data are written to flat ascii and
ifeffit PAD files.

  feff_tables.pl [OPTION]
    OPTIONS:
      --ascii, --noascii     write/suppress flat text database
      --pad, --nopad         write/suppress ifeffit PAD files
      -v, --version          print version information and exit
      -h, --help             print this help and exit

EOH
  ;
  exit;
};

__END__
 TITLE %n (%s %z) at the %e edge
 TITLE making central atoms phase shifts

 EDGE  %e   1.0
 CONTROL   1      1     1     1
 PRINT     1      0     0     3
 RMAX        2.51

 POTENTIALS
        0   %z   %s
        1   %z   %s

 ATOMS
    0.00000    0.00000    0.00000  0    %s              0.00000
    1.44338    1.44338    1.44338  1    %s_1            2.50000
   -1.44338    1.44338    1.44338  1    %s_1            2.50000
    1.44338   -1.44338    1.44338  1    %s_1            2.50000
   -1.44338   -1.44338    1.44338  1    %s_1            2.50000
    1.44338    1.44338   -1.44338  1    %s_1            2.50000
   -1.44338    1.44338   -1.44338  1    %s_1            2.50000
    1.44338   -1.44338   -1.44338  1    %s_1            2.50000
   -1.44338   -1.44338   -1.44338  1    %s_1            2.50000
    2.88675    0.00000    0.00000  1    %s_2            2.88675
   -2.88675    0.00000    0.00000  1    %s_2            2.88675
    0.00000    2.88675    0.00000  1    %s_2            2.88675
    0.00000   -2.88675    0.00000  1    %s_2            2.88675
    0.00000    0.00000    2.88675  1    %s_2            2.88675
    0.00000    0.00000   -2.88675  1    %s_2            2.88675
 END