#!/usr/bin/perl -ws

use List::Util qw(max min sum reduce);
#use Data::Dump qw(pp);

our ($h, $help, $terminal, $range, $linewidth, $BB, $nuc,
     $label, $join, $fontsize,
     $jcamp);

if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print "
Usage: $program [options] files.xyzppm
Зависимости: Perl, gnuplot

Рисует спектры ЯМР 13С по расчитанным хим. сдвигам с помощью gnuplot.
Исходный файл д.б. в формате xyz с 5-й колонкой (в ней хим.сдвиги).
Может быть и 6-я колонка, в ней обозначения атомов, если они есть,
то будут использоваться вместо порядковых номеров атомом.

Опции:
-h              help

-terminal=x11   Терминал для gnuplot. По умолчания - x11. Если терминал - 
                файловый формат (f.e. postscript), то перенаправить в файл.
                Список терминалов: echo 'help set terminal' | gnuplot

-nuc=C          Спектр на углероде, -nuc=H - на протонах (без расщеплений)

-BB             Обычный спектр (все сигналы вверх). По умолчанию - jmod (для С)

-range=min,max  Диапазон спектра (в м.д.). Если нет, то помещаются все сигналы

-linewidth=0.05 Ширина линии (в м.д.). 
                Можно давать в виде Hz/MHz, например 5/100.

-label=0        Не надписывать пики. Это по умолчанию.
-label          Надписывать пики обозначениями атомов.
-label=1        То же самое
-label=2        Надписывать пики их хим. сдвигами.
-label=3        Надписывать пики их хим. сдвигами и обозначениями атомов.

-join=N         Объединять надписи пиков, между которыми менее N ширин линий.
                По умолчанию N=1.

-fontsize=10    Размер шрифта для надписей пиков (используется шрифт cour).

-jcamp         На каждый спектр делать JCAMP-DX (файл с расширением .dx).
               -terminal=null, если не нужен рисунок от gnuplot.
  ";
  exit;
}

#############################################################################

# Терминал для gnuplot если не задан, то x11
$terminal ||= 'x11';

# Кол-во точек в спектре будет определяться автоматически
# $SI = 1000; # можно задать фиксированное значение

$nuc ||= 'C';

$BB = 1 if $nuc ne 'C';

unless ($linewidth) {
  $linewidth = 5/100; # Ширина линии (Гц/МГц)
  $linewidth = 1/200 if $nuc eq 'H';
}
if ($linewidth =~ /^(\d+(?:\.\d+)?)(?:\/(\d+(?:\.\d+)?))?$/) {
  die "Invalid linewidth\n" if $2 && $2==0;
  $linewidth = $1/$2 if $2;
} else {
  die "Invalid linewidth\n";
}

# Для линий, отстоящих менее чем на $join полуширин, объединять метки
$join ||= 1;

## Умолчаемое значение для $label
#$label = 1 unless defined $label;
#$label = 0 if $nolabel;

# размер шрифта для надписей пиков (используется шрифт cour)
$fontsize ||= 10;

# Относительная интенсивность синглета, дублета, триплета и квартета в JMOD
@intensivity = (-0.3, 0.7, -0.7, 0.7);
#@intensivity = (-0.2, 0.5, -0.5, 0.5);
# @intensivity = (0.7, 0.7, 0.7, 0.7); # все сигналы одинаковой интенсивности

# gnuplot (программа для рисования графиков)
$gnuplot = `which gnuplot`;
chomp $gnuplot;

##############################################################################

$minppm = 1e5;
$maxppm = -1e5;

my $label_ppm_format = '%.1f';
$label_ppm_format = '%.2f' if $nuc eq 'H';

# Берем из опции диапазон, в котором рисовать спектр
if (defined $range) {
  ($minppm, $maxppm) = split /[^\d.-]+/, $range;
  m/^-?\d+(\.\d+)?$/ || die "Invalid range\n" for ($minppm, $maxppm);
  die "Invalid range\n" if $minppm == $maxppm;
  ($minppm, $maxppm) = ($maxppm, $minppm) if $minppm > $maxppm;
}

# Зачитываем xyz-файлы с хим. сдвигами, каждый в свой датасет
$dataset = 0;
while (&read_molden) {
  $title[$dataset] = $ARGV; # Название спектра - название xyzppm-файла
  #$title[$dataset] =~ s/\.xyz(ppm)?$//;
  
  # Атомы углерода, у которых есть хим.сдвиги и их мультиплетность
  for ($i=1; $i<=$N; $i++) {
    if (uc($atom[$i]) eq $nuc && $ppm[$i]=~/^-?\d+\.\d+/) {
      my $mult = 1;
      if ($nuc eq 'C') {
        $mult = 0;
        for ($j=1; $j<=$N; $j++) {
          if ( uc($atom[$j]) eq 'H' &&
               ($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2+($z[$i]-$z[$j])**2 < 1.44
              ) {
            $mult++;
          }
        }
      }
      @{$PPM[$dataset]{$atom[$i].$i}} = (0+$ppm[$i], $mult, $user_label[$i]);
      
      unless (defined $range) {
        $minppm = $ppm[$i] if $ppm[$i] < $minppm;
        $maxppm = $ppm[$i] if $ppm[$i] > $maxppm;
      }
    }
  }
  
  $dataset++;
}

exit unless @PPM;

# Если не задан диапазон хим.сдвигов
unless (defined $range) {
  my $r = $maxppm-$minppm;
  $r = $linewidth*50 if $r == 0;
  $minppm = $minppm-$r/10;
  $maxppm = $maxppm+$r/10;
}

# Если не задано количество точек в спектре
unless ($SI) {
  my $p = int(($maxppm-$minppm)/$linewidth/2);
  $SI = (1,2,4,4,5,8,8,8,10,10)[substr $p, 0, 1] . '0'x(length($p)-1);
  $SI = 500   if $SI < 500;
  $SI = 10000 if $SI > 10000;
}

# ppm.out - графики спектров в табличной форме для gnuplot
unlink 'ppm.out' if -f 'ppm.out';
open OUT, ">>", "ppm.out" or die "Can't write to ppm.out: $!\n";

#pp @PPM;
my @spc;
foreach $dataset (0..$#PPM) {
  next unless $PPM[$dataset];
  my ($spc,$peak) = specrtum(
                      $PPM[$dataset], 
                      maxppm=>$maxppm, 
                      minppm=>$minppm, 
                      linewidth=>$linewidth, 
                      SI=>$SI,
                      baseoffset=>$dataset,
                      BB=>$BB,
                      intensivity=>\@intensivity,
                    );
  # $spc - ссылка на массив точек (x,y) спектра 
  # $spc->[x_точки, y_точки]
  # $peak - ссылка на хэш вершин пиков 
  # $peak->{atom_label}[номер_точки, x_вершины, y_вершины, user_label]
  #pp $peak;
  foreach (@{$spc}) {
    print OUT ${$_}[0], " \t", ${$_}[1], "\n";
  }
  print OUT "\n";
  #push @peak, $peak if $label != 0;
  push @peak, $peak;
  push @spc, $spc;
}
close OUT;

unless (-s 'ppm.out') {
  unlink 'ppm.out';
  exit;
}

if ($jcamp) {
  my @files = unique_array(@title);
  for (my $i=0; $i<@spc; $i++) {
    write_jcamp($spc[$i],$files[$i]);
#    my $file = $files[$i];
#    print "$file\n";
  }
}

# Если есть gnuplot, то будем рисовать спектры, а если нет,
# то рисовать не будем, но и удалять ppm.out тоже не будем
if ($gnuplot && $terminal ne 'null') {
  open GP, "|-", "$gnuplot -persist" or die "Can't pipe to gnuplot: $!\n";
  print GP "
  set xrange [$maxppm:$minppm] reverse
  set xlabel 'ppm'
  set xtics mirror
  set tics out
  set noytics
  set border 1
  set terminal $terminal
  ";
  my $maxpeak = ['','',-1000];
  my $minpeak = ['','',1000];
  my $str_max = '';
  my $str_min = '';
  foreach $dataset (0..$#PPM) {
    $ytitle = $dataset+0.1;
    (my $title = $title[$dataset]) =~ s/\.xyz(ppm)?$//;
    print GP "set label '$title, $nuc' at $maxppm,$ytitle rotate\n";
#    if ($label) {
      my @labels; # Массив меток [text,x,y]
      # Формируем массив меток
      foreach (keys %{$peak[$dataset]}) {
        #my $xlabel = $peak[$dataset]{$_}[1]-$linewidth; 
        my $xlabel = $peak[$dataset]{$_}[1]; 
        my $ylabel = $peak[$dataset]{$_}[2]; 
        push @labels, [($peak[$dataset]{$_}[3] || $_),$xlabel,$ylabel];
      }
      next unless @labels;
      
      # Объединяем метки по хим.сдвигам (у полученного массива другой формат)
      my @l = combine(\@labels, 1, $join*$linewidth);
      
      my @join_labels;
      foreach (@l) {
        my @a;
        foreach my $i (0..$#{$_->[0]}) {
          push @a, [ $_->[0][$i],$_->[1][$i],$_->[2][$i] ];
        }
        
        push @join_labels, combine(\@a, 2, 0.2);
      }
      
      # Формируем новые названия и расположения меток
      #pp @join_labels;
      foreach (@join_labels) {
        #pp $_;
        #next unless defined $_->[0][0];
        $_->[0] = join ',', sort {($A=$a)=~s/(\d+)/sprintf '%05d',$1/eg;
                                  ($B=$b)=~s/(\d+)/sprintf '%05d',$1/eg;
                                  $A cmp $B} 
                  @{$_->[0]};
        $_->[1] = sum(@{$_->[1]})/@{$_->[1]};
        $_->[2] = max(@{$_->[2]});
      }
      
      # C1,C2,C3 => C1,2,3
      foreach (@join_labels) {
        1 while $_->[0]=~s/((?:^|,)$nuc(?:\d+,)+)$nuc(\d+(?:,|$))/$1$2/g;
      }
      my $max_peak = reduce {$a->[2]>$b->[2] ? $a : $b} @join_labels;
      $maxpeak = $max_peak if $maxpeak->[2] < $max_peak->[2];
      my $min_peak = reduce {$a->[2]<$b->[2] ? $a : $b} @join_labels;
      $minpeak = $min_peak if $minpeak->[2] > $min_peak->[2];
      #pp $max_peak;
      if ($label) {
        if      ($label == 3) {
          $str_max = sprintf "%.2f (%s)", $max_peak->[1],$max_peak->[0];
          $str_min = sprintf "%.2f (%s)", $min_peak->[1],$min_peak->[0];
        } elsif ($label == 1) {
          $str_max = sprintf "%s", $max_peak->[0];
          $str_min = sprintf "%s", $min_peak->[0];
        } elsif ($label == 2) {
          $str_max = sprintf "%.2f", $max_peak->[1];
          $str_min = sprintf "%.2f", $min_peak->[1];
       }
        # Печатаем метки в gnuplot
        foreach (@join_labels) {
          my $direct = $_->[2]>$dataset ? 'left' : 'right';
          if      ($label == 3) {
            printf GP "set label '$label_ppm_format ($_->[0])' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n",$_->[1];
          } elsif ($label == 1) {
            print GP "set label '$_->[0]' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n";
          } elsif ($label == 2) {
            printf GP "set label '$label_ppm_format' at $_->[1],$_->[2] $direct rotate font 'cour,$fontsize'\n",$_->[1];
          }
        }
      }
    }
#  }
  
  my $symbolwigth = $fontsize/10*1.37; # mm
  my $scale = 100; # Высота y-шкалы в mm
  my $max_y = $#PPM;
  my $min_y = 0;
  $max_y = $intensivity[1] if $max_y < $intensivity[1] && $BB;
  $min_y = $intensivity[0] if $min_y > $intensivity[0] && !$BB;
  
  $max_y = $maxpeak->[2] if $maxpeak->[2] > $max_y;
  $min_y = $minpeak->[2] if $minpeak->[2] < $min_y;
  
  #print "$min_y:$max_y\n";
  if ($label) {
    #my $range = ($max_y - $min_y);
    my ($length_max,$length_min) = (0,0);
    $length_max += length($str_max) if $max_y > $#PPM;
    $length_min += length($str_min) if $min_y < 0;
    #print "length $length_max + $length_min, ", $symbolwigth*($length_max + $length_min)," mm\n";
    my $k = ($scale - $symbolwigth*($length_max + $length_min))/($max_y - $min_y);
    $k = 10 if $k < 10;
    #print "k = $k\n";
    $max_y += $symbolwigth*$length_max/$k if $max_y > $#PPM;
    $min_y -= $symbolwigth*$length_min/$k if $min_y < 0;
  }
  #print "$min_y:$max_y\n";
  my $offset = ($max_y - $min_y)/20;
  print GP "set yrange [$min_y:$max_y]\n";
#  print GP "set ytics\n";
  print GP "set offsets 0, 0, $offset, $offset\n";
  print GP "plot 'ppm.out' notitle with line\n";
  close GP;

#  unlink 'ppm.out' if -f 'ppm.out';
}

sub read_molden {
  return undef if eof();
  # 1-st string
  if (<>=~/^\s*(\d+)\s*$/) {
    ($N,$energy,@atom,@x,@y,@z,@ppm,@user_label) = undef;
    $N = $1;
  } else {
    return undef;
  }
  # 2-nd string
  ($energy) = <>=~/(-?\d+\.\d+)/;
  # Geometry
  for ($i=1;$i<=$N;$i++) {
    ($atom[$i],$x[$i],$y[$i],$z[$i],$ppm[$i],$user_label[$i]) = split ' ', <>;
    return undef unless $atom[$i]=~/^\w\w?/ &&
                        $x[$i]=~/^-?\d+\.\d+/ &&
                        $y[$i]=~/^-?\d+\.\d+/ &&
                        $z[$i]=~/^-?\d+\.\d+/;
    $user_label[$i] =~ s/(['";,])/\\$1/g if $user_label[$i];
  }
  return $N;
}

sub specrtum {
  my $ppm = shift;
  my %param = (
    linewidth => 0.02,
    maxppm => 250,
    minppm => 0,
    SI => '8000',
    height => 0.7,
    baseoffset => 0,
    level => 0.001,
    BB => 0,
    intensivity => [-0.3, 0.7, -0.7, 0.7],
    @_
  );
  my $debug = 0;
  if ($debug) {
    print "%param:";
    print "\t$_ => $param{$_}\n" foreach (sort keys %param);
  }
  my $SW = $param{maxppm} - $param{minppm};
  $SW = $param{linewidth} if $SW==0;
  my $d = $SW / $param{SI};
  print "SW = $SW ppm, ppm per point = $d\n" if $debug;
  my @SPC;
  foreach (0..$param{SI}) {
    $SPC[$_][0] = $param{minppm} + $_*$d;
    $SPC[$_][1] = $param{baseoffset};
  }
  
  my @halfline;
  my $b = $param{linewidth}/2;
  my $y = 1;
  my $dd = $d;
  while ($y > $param{level}) {
    push @halfline, $y;
    my $t = $dd/$b;
    $y = 1 / (1 + $t*$t);
    $dd += $d;    
  }
  print "halfline: @halfline\n" if $debug;
  
  my %peak;
  foreach (keys %$ppm) {
    next if $ppm->{$_}[0] < $param{minppm} || $ppm->{$_}[0] > $param{maxppm};
    my $x = $ppm->{$_}[0];
    
    my $mult = $param{intensivity}->[$ppm->{$_}[1]];
    $mult = abs($mult) if $param{BB};
    
    # Помещаем вершину линии на ближайшую точку
    my $I = sprintf "%d", ($x - $param{minppm}) / $d;
    foreach (-$#halfline..$#halfline) {
      my $j = $I + $_;
      next if $j<0 or $j>$SI;
      $SPC[$j][1] += $halfline[abs($_)]*$mult; 
    }
    $peak{$_}[0] = $I if $I>=0 && $I<=$SI;
    
    $peak{$_}[3] = $ppm->{$_}[2];
  }
  
  foreach (keys %peak) {
    $peak{$_}[1] = $SPC[$peak{$_}[0]][0];
    $peak{$_}[2] = $SPC[$peak{$_}[0]][1];
    print "$_: point $peak{$_}[0], x=$peak{$_}[1], y=$peak{$_}[2]\n" if $debug;
  };
  
  return (\@SPC, \%peak);
}

sub combine {
  my @a = @{$_[0]};
  my $field = $_[1];
  my $criter = $_[2];
  
  @a = sort {$a->[$field]<=>$b->[$field]} @a;
  
  my @b = [ [$a[0][0]],[$a[0][1]],[$a[0][2]] ];
  for (my $i=1; $i<@a; $i++) {
    if ($a[$i][$field]-$a[$i-1][$field] <= $criter) {
      push @{$b[-1][0]}, $a[$i][0];
      push @{$b[-1][1]}, $a[$i][1];
      push @{$b[-1][2]}, $a[$i][2];
    } else {
      push @b, [ [$a[$i][0]],[$a[$i][1]],[$a[$i][2]] ];
    }
  }
  return @b;
}

sub write_jcamp {
  my ($spc, $file) = @_;
  
  my ($sec,$min,$hour,$mday,$mon,$year) = 
     localtime(-f $file ? (stat $file)[9] : time());
  my $DATE = sprintf "%02d/%02d/%02d", $year % 100, $mon+1, $mday;
  my $TIME = sprintf "%02d:%02d:%02d", $hour,$min,$sec;
  
  open F, '>', "$file.dx" or die "Can't write $file.dx: $!\n";
  print F "##TITLE= Calculated $nuc NMR from $file
##JCAMP-DX= 4.24
##DATA TYPE= NMR SPECTRUM
##ORIGIN= $0
##OWNER= $ENV{USER}
##DATE= $DATE
##TIME= $TIME
##DELTAX= @{[$spc->[1][0]-$spc->[0][0]]}
##XUNITS= PPM
##YUNITS= ABSORBANCE
##XFACTOR= 1.00
##YFACTOR= 1.00
##FIRSTX= $spc->[0][0]
##LASTX= $spc->[-1][0]
##NPOINTS= @{[$#{$spc}+1]}
##FIRSTY= @{[$spc->[0][1]*1e8]}
##XYDATA= (X++(Y..Y))
";
##.OBSERVE FREQUENCY= 75.476050452
##.OBSERVE NUCLEUS= ^13C
##SPECTROMETER/DATA SYSTEM= Prirda, Dimitri Laikov

  foreach (@$spc) {
    printf F "  %11.6f  %12d\n", $_->[0], $_->[1]*1e8;
  }
  print F "##END= \n";
  close F;
}

sub unique_array {
  my @a = @_;
  my %h;
  
  $h{$_}++ foreach @a;
  
  for(my $i = $#a; $i>=0; $i--) {
    next if $h{$a[$i]}==1;
    my $j = --$h{$a[$i]};
    $j++ while exists $h{"$a[$i]-$j"};
    $a[$i] .= "-$j";
    $h{$a[$i]}++;
  }
  return @a;
}