#!/usr/bin/perl -w

use Math::Trig;

# Извлекаем опции из аргументов
while (<@ARGV>) {
  if (/^-h$/ or /^--help$/) {$help=1; last}
  if (/^-(\S+)$/) {
    push (@opts, $1);
    shift;
  }
}

if ($help) {
  $0 =~ s/^.*[\/\\]//;
  print <<HELP;
Do some simplest transformations of cartesian coordinates 
given in molden format file(s)

Usage: $0 [-options] [filename] [filenames] [> result]

Зависимости: Perl, Math::Trig

filename - файл с декартовыми координатами в формате molden'а;
           если отсутствует, то stdin
Опции:
h            этот help
x|y|z        поворот на 180 град. вокруг оси x (y, z)
xN|yN|zN     поворот на N град. вокруг оси x (y, z)
xy|yz|xz     отражение относительно плоскости xy (yz, xz). Default -xy 
xyz          инверсия относительно центра координат
c            преобразовать координаты к геометрическому центру
sxN|syN|szN  сдвиг по оси x (y, z) на N ангстрем

HELP
  exit;
}

# Regexp for number
my $d = qr(-?\d+(?:\.\d+)?);

# Проверка на валидность
if (@opts) {
  foreach (@opts) {
    unless (/^s?[xyz](?:$d)?$/ || /^(xy|yz|xz)$/ || /^xyz$/ || /^c$/) {
      die "Invalid option: $_\n";
    }
    # Нужно ли преобразовать координаты к центру
    $centre = 1 if /^c$/;
  }
} else {
  $opts[0] = 'xy';
}

# Удаляем "c" из опций
@opts = grep {!/^c$/} @opts if $centre;

# Поехали
while (<>) {
  # 1-st string
  if (/^\s*(\d+)\s/) {
    $N = $1;
    print;
  }
  else {
    print "Not molden format";
    exit;
  }
  # 2-nd string
  $_ = <>;
  print;

  # Зачитываем атомы и кординаты
  (@atom,@x,@y,@z,@rest) = undef;
  for ($i=1;$i<=$N;$i++) {
    $_ = <>;
    ($atom[$i],$x[$i],$y[$i],$z[$i],$rest[$i]) = split;
  }
  
  # Перемещаем в центр, если нужно
  if ($centre) {
    $xc = 0;                 $yc = 0;                 $zc = 0;
    $xc += $_ for @x[1..$N]; $yc += $_ for @y[1..$N]; $zc += $_ for @z[1..$N];
    $xc /= $N;               $yc /= $N;               $zc /= $N;
    $_ -= $xc for @x[1..$N]; $_ -= $yc for @y[1..$N]; $_ -= $zc for @z[1..$N];
  }
  
  for ($i=1;$i<=$N;$i++) {
    foreach (@opts) {
      if    ($_ eq 'x'  ) {             $y[$i]*=-1; $z[$i]*=-1; }
      elsif ($_ eq 'y'  ) { $x[$i]*=-1;             $z[$i]*=-1; }
      elsif ($_ eq 'z'  ) { $x[$i]*=-1; $y[$i]*=-1;             }
      elsif ($_ eq 'xy' ) {                         $z[$i]*=-1; }
      elsif ($_ eq 'yz' ) { $x[$i]*=-1;                         }
      elsif ($_ eq 'xz' ) {             $y[$i]*=-1;             }
      elsif ($_ eq 'xyz') { $x[$i]*=-1; $y[$i]*=-1; $z[$i]*=-1; }
      elsif (m/^([xyz])($d)$/) {
        # warn "Rotate around $1 on $2 degrees\n";
        $cos = cos(deg2rad($2)); 
        $sin = sin(deg2rad($2));
        $x = $x[$i]; $y = $y[$i]; $z = $z[$i];
        # Rotate around x
        if ($1 eq 'x') {
          $z[$i] = $cos*$z-$sin*$y;
          $y[$i] = $sin*$z+$cos*$y;
        }
        # Rotate around y
        if ($1 eq 'y') {
          $x[$i] = $cos*$x-$sin*$z;
          $z[$i] = $sin*$x+$cos*$z;
        }
        # Rotate around z
        if ($1 eq 'z') {
          $x[$i] = $cos*$x-$sin*$y;
          $y[$i] = $sin*$x+$cos*$y;
        }
      }
      elsif (m/^s([xyz])($d)$/) {
        $x[$i] += $2 if $1 eq 'x';
        $y[$i] += $2 if $1 eq 'y';
        $z[$i] += $2 if $1 eq 'z';
      }
    }
    printf " %-2s %12.8f %12.8f %12.8f",$atom[$i],$x[$i],$y[$i],$z[$i];
    print "   $rest[$i]" if $rest[$i];
    print "\n";
  }
}