#!/usr/bin/perl
# db_use: calculate radioactive decay chain. Loads ENSDF data converted
# by db_create to Storable files. Prints decay tree in text or
# graphviz input files for generating images.
# $Id: db_use 57 2009-12-28 10:17:54Z hatta $

use strict;
use Storable;
use Memoize;

my $decays_table = retrieve('ensdf-decays.dat');
my $gamma_transitions = retrieve('ensdf-gamma-transitions.dat');

# parse command line
my $debug = 0;
my $dot = 0;
my $show_table = 0;
my $show_reactions = 0;
my $show_intensity = 0;
my $rankdir_lr = 0;
my $show_gammas = 0;
my $show_isomers = 0;
unless (@ARGV > 0) {
  die "Usage: $0 <filename> {flags} or db_use <a> <z> {flags}\n";
}
my $mode;
my @input_nuclei;
my $flags_begin = 0;
my $input_text;
if (-r $ARGV[0]) {
  $mode = "file";
  $flags_begin = 1;
  open(INPUT, "< $ARGV[0]");
  while(<INPUT>) {
    my ($a, $z, $yield, $abundance, @fields) = split;
    next unless ($a > 0 && $z > 0);
    push @input_nuclei, {a => $a, z => $z, yield => $yield, abundance => $abundance, text => \@fields};
    $input_text .= $_;
  }
} elsif ($ARGV[0] =~ /^\d+$/ && $ARGV[1] =~ /^\d+$/) {
  $mode = "number";
  $flags_begin = 2;
  push @input_nuclei, {a => $ARGV[0], z => $ARGV[1], text => []};
}

foreach (@ARGV[$flags_begin .. $#ARGV]) {
  $dot = 1 if /dot/i;
  $debug = 1 if /debug/i;
  $show_table = 1 if /table/i;
  $show_intensity = 1 if /intensity/i;
  $show_reactions = 1 if /reactions/i;
  $rankdir_lr = 1 if /rankdir/i;
  $show_gammas = 1 if /gamma/i;
  $show_isomers = 1 if /isomer/i;
}

my $fontspec = "fontname=Helvetica,fontsize=10";
my $nodespec = "fixedsize=false,height=0.2,width=0.7,shape=box,$fontspec,labeljust=l";

my %ztonam;
my %namtoz;
open(ELEMENTS, "< elements.txt") || die "can't open elements file.\n";
while (<ELEMENTS>) {
  if (/(\d+)\s*(\w+)/) {
    $ztonam{$1} = lc($2);
    $namtoz{lc($2)} = $1;
  }
}
close(ELEMENTS);

sub find_levels {
  my $levels = shift;
  my $energy = shift;
  my @res = sort {abs($energy - $a) <=> abs($energy - $b)} keys %$levels;
#  print "find_levels: given $energy return $res[0]\n";
  return ($res[0] > 0 ? $res[0] : "0.0");
}

sub find_gamma_transitions {
  my ($table, $ilevel, $p) = @_;
  my $metastable = $$table{metastable};
  my $levels_table = $$table{levels};
  my $gamma_table = $$table{gammas};
  my $t12 = $$table{halflife};
  my $t12_sec = $$table{halflife_sec};
  my $normalization = $$table{normalization};
  my %states;
  my $ilevel1 = find_levels($levels_table, $ilevel);
  # Either: a) this level is metastable
  if ($$metastable{$ilevel1} =~ /M/
      || $$t12{$ilevel1} =~ /STABLE/i
      || $$t12_sec{$ilevel1} > 0.1) {
    return {$ilevel1 => 1};
  }
  # b) this level is g.s.
  unless ($ilevel1 > 0) {
    return {"0.0" => 1};
  }
  # c) this level has no exits and is not metastable
  my @g = @{$$gamma_table{$ilevel1}};
  if(scalar(@g) == 0) {
    return {};
  }
  # d) this level is not metastable and has exits
  my $sum_i = 0;
  my $ti_flag = 0;
  $sum_i += $$_{ri} foreach (@g);
  if ($sum_i == 0) {
    $sum_i += $$_{ti} foreach (@g);
    $ti_flag = 1;
  }
  if ($sum_i == 0 and scalar(@g) > 0) {
    die "error: sum_i is 0 and transitions from $ilevel1 exist.\n";
    $sum_i = 1;
  }
  if ($$normalization{NRBR} > 0) {
    print "NRBR=$$normalization{NRBR}.\n";
  }
  for my $g (@g) {
    my $flevel = find_levels($levels_table, $ilevel1 - $$g{energy});
    my $fri = ($ti_flag ? $$g{ti} : $$g{ri}) / $sum_i;
#    printf "$ilevel1 to $flevel RI = $$g{ri}, TI = $$g{ti}, NRBR = $$normalization{NRBR}\n" if $p == 0;
    $sum_i > 0 or print "sum_i is 0, $ilevel1 to $flevel\n";
    my $sss = find_gamma_transitions($table, $flevel, $p+2);
    for my $s (keys %$sss) {
      if ($s == 0) {$s = "0.0";}
      $states{$s} += $fri*$$sss{$s};
    }
  }
  # print "-----------------\n";
  # if ($p == 0) {
  #   print "$_, $states{$_}\n" foreach (keys %states);
  # }
  return \%states;
}

# find_isomeric_decay: search for missing IT data in ADOPTED LEVELS
sub find_isomeric_decay {
  my ($a, $z, $ilevel, $name) = @_;
  my $gt = $$gamma_transitions{"$a-$z"};
  my $dt = $$decays_table{"$a-$z"};
  $ilevel = find_levels($$gt{levels}, $ilevel);
  if ($gt->{metastable}->{$ilevel} =~ /M/ 
      || $gt->{halflife_sec}->{$ilevel} > 0.1 
      || $gt->{halflife}->{$ilevel} =~ /STABLE/i) {
#    print "$a-$z $ilevel is metastable\n";
  } else {
    return undef;
  }
  my %a;
  my %i;
  for my $g (@{$gt->{gammas}->{$ilevel}}) {
    my $l = find_levels($$gt{levels}, $ilevel - $$g{energy});
    my $t = find_gamma_transitions($gt, $l, 0);
    %i = %$t;
    $a{$_} = 1 for (keys %$t);
  }
  my @f = keys %a;
  my @g;
  for my $a (@f) {
    push @g, {energy => $ilevel - $a, intensity => $i{$a}};
  }
  foreach my $d (@$dt) {
    my $decay_parent_level = find_levels($$gt{levels}, $$d{P});
    if ($$d{type} eq "isomeric" && abs($decay_parent_level - $ilevel) < 0.1) {return undef};
  }
  return {file => $$gt{file},
	  parent_a => $a,
	  parent_z => $z,
	  parent_name => $name,
	  final_a => $a,
	  final_z => $z,
	  final_name => $name,
	  type => "isomeric2",
	  P => $ilevel,
	  halflife => $gt->{halflife}->{$ilevel},
	  halflife_sec => $gt->{halflife_sec}->{$ilevel},
	  G => \@g,
	  L => \@f,
	  I => \%i,
	  NR => 1,
	  NT => 1,
	  BR => 1};
}

my %dot;
my %dot_nodes;
my %initial_nuclei;
my %dot_edges;
my %all_gammas_lines;
my %decay_constant_table;
sub add_decay_constant {
  my ($parent, $final, $c) = @_;
  my $row;
  if(exists $decay_constant_table{$parent}) {
    $row = $decay_constant_table{$parent};
  } else {
    $row = {};
  }
#   print STDERR "<small><i><b>warning:</b> add decay constant $c of $parent to $final: ", 
#     "already exists $$row{$final}.</i></small><br>\n" 
#       if(exists $$row{$final} && abs($$row{$final} - $c) > 2*abs($$row{$final}));
  $$row{$final} = $c;
  $decay_constant_table{$parent} = $row;
}

sub add_node {
  my ($nucl, $level, $spin_parity, $text) = @_;
  my $node_name = "$nucl";
  if ($level > 0) {
    $node_name .= " (${level}m)";
  }
  my $node = {name => $node_name,
	      level => $level,
	      spin_parity => $spin_parity,
	      text => $text};
  if (exists $dot_nodes{$nucl}) {
    my $nodes = $dot_nodes{$nucl};
    foreach my $dn (@$nodes) {
      return if $$dn{name} eq $node_name;
    }
    push @$nodes, $node;
  } else {
    $dot_nodes{$nucl} = [$node];
  }
}

sub add_edge {
  my ($name, $opt) = @_;
  if (exists($dot_edges{$name})) {
    foreach (@{$dot_edges{$name}}) {
      return if $opt eq $_;
    }
    push @{$dot_edges{$name}}, $opt;
  } else {
    $dot_edges{$name} = [$opt];
  }
}

sub add_gamma_line {
  my ($e, $text) = @_;
  if (exists $all_gammas_lines{$e}) {
    my $good = 1;
    foreach (@{$all_gammas_lines{$e}}) {
      $good = 0 if $text eq $_;
    }
    push @{$all_gammas_lines{$e}}, $text if $good;
  } else {
    $all_gammas_lines{$e} = [$text];
  }
}

sub trace_decays {
  my ($a, $z, $ilevel, $ppp) = @_;
  # find all decay paths from here
  my $gt = $$gamma_transitions{"$a-$z"};
  my $dt = $$decays_table{"$a-$z"};
  $ilevel = find_levels($$gt{levels}, $ilevel);
  my $parent_initial_levels = find_gamma_transitions($gt, $ilevel, $ppp);
  my @parent_initial_levels = keys %$parent_initial_levels;
  # a) no known decays
  return unless defined($dt);
  my @decays = @$dt;
  my $isomeric_decay = find_isomeric_decay($a, $z, $ilevel, $decays[0]->{parent_name});
  if (defined($isomeric_decay)) {
    my $sum_br = 0;
    for my $d (@decays) {
      for my $s (@parent_initial_levels) {
	my $stable_c = $$parent_initial_levels{$s};
	if (abs($$d{P} - $s) <= 0.0001) {
	  $sum_br += $$d{BR} * $stable_c * $$d{I}->{$s};
	}
      }
    }
    $$isomeric_decay{BR} = 1 - $sum_br;
    push @decays, $isomeric_decay if defined($isomeric_decay);
  }
  # b) iterate over all decay paths of this nucleus
  for my $d (@decays) {
#    if ($$d{parent_name} =~ /109PD/) { print "$$d{final_name}, $$d{P}, $$d{file}, $$d{type}, $ilevel\n"; }
    # over all final stable levels (possible from initial level)
    for my $s (@parent_initial_levels) {
      my $stable_c = $$parent_initial_levels{$s};
      my $decay_parent_level = find_levels($$gt{levels}, $$d{P});
      if (abs($decay_parent_level - $s) <= 0.0001) {
	# $d goes through all decays of ($a, $z, $ilevel)
	my %final_final_levels;
	my $final_level_intensities = $$d{I};
	# iterate over all final levels of this decay ($d)
	for my $f (@{$$d{L}}) {
	  my $intensity = $$final_level_intensities{$f};
	  my $final_final_levels_tmp =
	    find_gamma_transitions($$gamma_transitions{"$$d{final_a}-$$d{final_z}"}, $f, 0);
# 	  print %{$$gamma_transitions{"$$d{final_a}-$$d{final_z}"}}->{file}, ": ";
# 	  print "$f ->...-> $_  ($$tmp2{$_}); " foreach keys %$tmp2;
# 	  print "\n";
	  my @final_final_levels_tmp = keys %$final_final_levels_tmp;
	  foreach my $t4 (@final_final_levels_tmp) {
	    $final_final_levels{$t4} += $$final_final_levels_tmp{$t4}*$stable_c*$intensity;
	    #print "Iteration: $$final_final_levels_tmp{$t4} x $stable_c x $intensity.\n";
	  }
	}
#	print "At last:\n";
#	print "$_ $final_final_levels{$_};\n" foreach keys %final_final_levels;
	my @final_final_levels = sort keys %final_final_levels;
	for my $ff (@final_final_levels) {
	  my $parent_name = $$d{parent_name};
	  $parent_name .= " (${s}m)" if $s > 0;
	  my $final_name = $$d{final_name};
	  $final_name .= " (${ff}m)" if $ff > 0;
	  unless ($$d{parent_name} eq $$d{final_name} && $ff eq $s) {
	    # determine spin-parity of parent and final states
	    my $parent_gamma_table = $$gamma_transitions{"$$d{parent_a}-$$d{parent_z}"};
	    my $parent_sp = $$parent_gamma_table{spin_parity}->{$s};
	    my $final_gamma_table = $$gamma_transitions{"$$d{final_a}-$$d{final_z}"};
	    my $final_sp = $$final_gamma_table{spin_parity}->{$ff};
	    (my $halflife = $$d{halflife}) =~ s/^\s*(.*?)\s*$/$1/;
	    my $intensity_str =
	      sprintf("%.6f%% x %.6f%%", $$d{BR}*100, $final_final_levels{$ff}*100);
	    # record a photon coming from this decay
	    my $gammas = $$d{G};
	    foreach my $gamma (@$gammas) {
	      my $ti = $$gamma{intensity};
	      my $nt = $$d{NT};
	      $ti *= $nt if $nt > 0;
	      add_gamma_line($$gamma{energy}, "$ti ($parent_name $$d{type})") if $ti > 0.1;
	    }
	    print " "x$ppp,
	      "$$d{type} decay of $parent_name to $final_name, ",
		"intensity $intensity_str file $$d{file}\n" unless $dot;
	    if ($dot) {
	      my $opt = "$$d{type}\\n$halflife";
	      $intensity_str = sprintf("%.2f%%", $$d{BR}*$final_final_levels{$ff}*100);
	      $opt .= "\\n$intensity_str" if $show_intensity;
	      add_edge(qq("$parent_name" -> "$final_name"),
		       qq($fontspec,label="$opt"));
	      add_node($$d{parent_name}, $s, $parent_sp, []);
	      add_node($$d{final_name}, $ff, $final_sp, []);
	    }
	    if ($$d{halflife_sec} > 0) {
	      my $decay_constant = 0.69314718056 / $$d{halflife_sec};
	      $decay_constant *= $$d{BR} * $final_final_levels{$ff};
	      add_decay_constant($parent_name,$final_name, $decay_constant);
	    }
	    trace_decays($$d{final_a}, $$d{final_z}, $ff, $ppp+2);
	  } else {
	    print " "x$ppp,
	      "skipping $$d{type} decay of $parent_name ",
		"to $final_name file $$d{file}\n" unless $dot;
	  }
	}
      }
    }
  }
}

memoize('trace_decays');
memoize('find_gamma_transitions');
memoize('find_levels');
memoize('find_isomeric_decay');

sub print_nodes {
  my $nodname = shift;
  my @nodes = @{$dot_nodes{$nodname}};
  if (scalar @nodes > 1 ) {
    print "subgraph cluster$nodname {\n";
    print "compound=true;\n";
    my %text;
    foreach my $dn (@{$dot_nodes{$nodname}}) {
      foreach my $t (@{$$dn{text}}) {
	$text{$t} = 1;
      }
    }
    my $text;
    my $initial;
    foreach my $t (sort keys %text) {
      $text .= qq(<br/><font point-size="8">$t</font>);
      if ($t =~ /initial/i) {$initial = 1;}
    }
    print qq(label=<<font point-size="14">$nodname</font>$text>;\n);
    print "style=filled;\nfillcolor=grey;\n" if $initial;
    foreach my $dn (@{$dot_nodes{$nodname}}) {
      my $style = "shape=plaintext";
      if ($$dn{level} > 0) {
	$style = "shape=box,style=filled,fillcolor=pink";
      }
      my $level;
      if ($$dn{level} > 0) {
	$level = "$$dn{level} keV";
      } else {
	$level = "g. s.";
      }
      print qq("$$dn{name}" [$style,label=<<font point-size="12">$level</font>$$dn{spin_parity}>];\n);
    }
    print "}\n";
  } else {
    foreach my $dn (@{$dot_nodes{$nodname}}) {
      my $label = qq(<font point-size="14">$$dn{name}</font>$$dn{spin_parity});
      my $initial;
      foreach my $t (@{$$dn{text}}) {
	$label .= qq(<br/><font point-size="8">$t</font>);
	if ($t =~ /initial/i) {$initial = 1;}
      }
      my $opt;
      $opt = ",style=filled,fillcolor=grey" if $initial;
      print qq("$$dn{name}" [label=<$label>$opt];\n);
    }
  }
}
my $input_count = scalar @input_nuclei;
my %yields;
my %abundances;
unless($debug) {
  my $dot_size = "15,20";
  my @nuclei_for_tracing;
  my $initial_count = 1;
  for my $input (@input_nuclei) {
    my $aaa = $$input{a};
    my $zzz = uc($ztonam{$$input{z}});
    my $yield = $$input{yield};
    my $abundance = $$input{abundance};
    $yields{"$aaa$zzz"} = $yield;
    $abundances{"$aaa$zzz"} = $abundance;
    my $gamma_table = $$gamma_transitions{"$$input{a}-$$input{z}"};
    my $spin_parity= $$gamma_table{spin_parity}->{"0.0"} || $$gamma_table{spin_parity}->{0};
    my $initial = 0;
    my $strictly_initial = 1;
    # find all isomers of this nucleus if show_isomers is true
    my @stable_levels;
    foreach my $level (keys %{$$gamma_table{levels}}) {
      next unless $level > 0;
      if ($show_isomers) {
	if ($$gamma_table{halflife}->{$level} =~ /STABLE/i
	    || $$gamma_table{halflife_sec}->{$level} > 0.1
	    || $$gamma_table{metastable}->{$level} =~ /M/) {
	  push @stable_levels, {a => $$input{a}, z => $$input{z}, level => $level, spin_parity => $$gamma_table{spin_parity}->{$level}};
	}
      }
    }
    my $label_text;
    my $reactions;
    my @text;
    foreach my $t (@{$$input{text}}) {
      if ($t =~ /^(\d+\w+)(\(.*\))(\d+\w+)$/) {
	push @text, $t unless $show_reactions;
	$strictly_initial = 0;
      } else {
	push @text, $t;
      }
      $initial = 1 if $t =~ /^initial$/i;
    }
    $initial_nuclei{"$aaa$zzz"} = 1 if $initial;
    push @stable_levels, {a => $$input{a}, z => $$input{z}, level => 0, spin_parity => $spin_parity};
    foreach my $s (@stable_levels) {
      next if ($strictly_initial && $initial && $$s{level} > 0);
      my $name = "$aaa$zzz";
      if ($$s{level} > 0) {
	$name .= " ($$s{level}m)";
      }
      if ($initial && !($$s{level} > 0)) {
	$initial_count++;
      }
      add_node("$aaa$zzz", $$s{level}, $$s{spin_parity}, \@text);
      push @nuclei_for_tracing, $s;
      foreach my $t (@{$$input{text}}) {
	if ($t =~ /^(\d+\w+)(\(.*\))(\d+\w+)$/) {
	  add_edge(sprintf(qq("%s" -> "%s"), uc($1), $name), qq($fontspec,label="$2")) if $show_reactions;
	}
      }
    }
  }
  trace_decays($$_{a}, $$_{z}, $$_{level}, 0) foreach @nuclei_for_tracing;
  $dot_size="0,0";
  if ($dot) {
    my $rankdir = $rankdir_lr ? "LR" : "TB";
    print <<EOF;
digraph decays {
size="$dot_size";
ratio=compress;
nodesep=0.1;
ranksep=0;
rankdir="$rankdir";
overlap=false;
node [$nodespec];
edge [$fontspec];
fontname=Helvetica;
fontsize=10;
EOF
#    print "{\n\trank=same;\n";
    foreach my $dnh (keys %dot_nodes) {
      if ($initial_nuclei{$dnh} > 0) {
	print_nodes($dnh);
      }
    }
#    print "}\n";
    foreach my $dnh (keys %dot_nodes) {
      unless ($initial_nuclei{$dnh}) {
	print_nodes($dnh);
      }
    }
    foreach my $edge (keys %dot_edges) {
      foreach my $opt (@{$dot_edges{$edge}}) {
	print "$edge [$opt];\n";
      }
    }
    print "}\n";
  }
  if ($show_gammas) {
    my $html = "<h3>Gamma spectral lines</h3><table><tr><th>Energy, keV</th><th>Comment</th></tr>\n";
    foreach my $energy (sort {$a <=> $b} keys %all_gammas_lines) {
      foreach my $c (@{$all_gammas_lines{$energy}}) {
	$html .= "<tr><td>$energy</td><td>$c</td></tr>\n";
      }
    }
    $html .= "</table>\n";
    print STDERR $html;
  }
}

if($debug) {
  my $gamma_table = $$gamma_transitions{"234-91"};
  my $final = find_gamma_transitions($gamma_table, 166.72, 0);
  print "...-> $_  ($$final{$_}); " foreach keys %$final;
  print "\n";
}

sub make_applet_html {
  my ($class, $param) = @_;
  my $applet_html =<<EOF
<applet codebase="applet" archive="applet.jar,mantissa-7.2.jar" code="Main" width="100%" height="90%">
$param
</applet>
EOF
    ;
  return $applet_html;
}

# print decay constants table
if($show_table) {
  my @parents = keys %decay_constant_table;
  my %finals;
  my @finals;
  my $pnum = scalar @parents;
  my $fnum = 0;
  foreach my $p (@parents) {
    my $fnum1 = scalar(keys(%{$decay_constant_table{$p}}));
    $fnum = $fnum1 if $fnum1 > $fnum;
    $finals{$_} = 1 foreach keys %{$decay_constant_table{$p}};
  }
  @finals = keys %finals;
  # now print the same but for java applet
  print STDERR qq{<br><br><br><center>\n};
  my %all_names;
  $all_names{$_} = 1 foreach @parents;
  $all_names{$_} = 1 foreach @finals;
  my @all_names = keys %all_names;
  my $param_html = qq!<param name=constant value="{constant: [!;
  my $param2_html = qq!<param name=y value="{y: [!;
  my $param3_html = qq!<param name=n0 value="{n0: [!;
  my $param4_html = qq!<param name=names value='{names: [!;
  my @l;
  my @y;
  my @a;
  my @names;
  foreach my $p (@all_names) {
    my @c;
    foreach my $f (@all_names) {
      my $c = ${$decay_constant_table{$p}}{$f} || 0;
      push @c, $c;
    }
    my $s = join ", ", @c;
    push @l, "[$s]";
    my $y = ($yields{$p} > 0) ? $yields{$p} : 0;
    push @y, $y;
    my $a = ($abundances{$p} > 0 && $initial_nuclei{$p} > 0) ? $abundances{$p} : 0;
    push @a, $a;
    push @names, qq("$p");
  }
  $param_html .= join ",\n", @l;
  $param_html .= qq!]}">\n!;
  $param2_html .= join ", ", @y;
  $param2_html .= qq!]}">\n!;
  $param3_html .= join ", ", @a;
  $param3_html .= qq!]}">\n!;
  $param4_html .= join ", ", @names;
  $param4_html .= qq!]}'>\n!;
  print STDERR make_applet_html("Main.class", "$param_html\n$param2_html\n$param3_html\n$param4_html");
  print STDERR "</object></center>";
  print STDERR "<!-- Program input dump:\n$input_text\n -->\n";
}