#!/bin/bash

echo "parsing POSCAR file ..."

theta=$1

trg=POSCAR.Xtrafo$theta

sicnt=`sed -n 6p POSCAR | awk '{ print $1 }'`
ccnt=`sed -n 6p POSCAR | awk '{ print $2 }'`

lc=`sed -n 2p POSCAR | awk '{ print $1 }'`

x1=`sed -n 3p POSCAR | awk '{ print $1 }'`
x2=`sed -n 3p POSCAR | awk '{ print $2 }'`
x3=`sed -n 3p POSCAR | awk '{ print $3 }'`

y1=`sed -n 4p POSCAR | awk '{ print $1 }'`
y2=`sed -n 4p POSCAR | awk '{ print $2 }'`
y3=`sed -n 4p POSCAR | awk '{ print $3 }'`

z1=`sed -n 5p POSCAR | awk '{ print $1 }'`
z2=`sed -n 5p POSCAR | awk '{ print $2 }'`
z3=`sed -n 5p POSCAR | awk '{ print $3 }'`

((total=sicnt+ccnt))

echo "  Si: $sicnt, C: $ccnt, total: $total"
echo "  Lattice constant: $lc A"
echo "  Basis:"
echo "      $x1      $y1      $z1"
echo "  x = $x2  y = $y2  z = $z2"
echo "      $x3      $y3      $z3"
echo
echo "  -----"
echo "  Trafo: x axis, $theta degree"
echo "  -----"
echo

# determine trafo of cartesian coordinates to this one
normx=`echo $x1 $x2 $x3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
normy=`echo $y1 $y2 $y3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
normz=`echo $z1 $z2 $z3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
echo $normx $normz $normy
t11=`echo $x1 $normx | awk '{ print $1/$2 }'`
t12=`echo $x2 $normx | awk '{ print $1/$2 }'`
t13=`echo $x3 $normx | awk '{ print $1/$2 }'`
t21=`echo $y1 $normy | awk '{ print $1/$2 }'`
t22=`echo $y2 $normy | awk '{ print $1/$2 }'`
t23=`echo $y3 $normy | awk '{ print $1/$2 }'`
t31=`echo $z1 $normz | awk '{ print $1/$2 }'`
t32=`echo $z2 $normz | awk '{ print $1/$2 }'`
t33=`echo $z3 $normz | awk '{ print $1/$2 }'`
echo "  Matrix from cartesian to used coordinates:"
echo "  | $t11 $t21 $t31 |"
echo "  | $t12 $t22 $t32 |"
echo "  | $t13 $t23 $t33 |"
echo
i11=$t11
i12=$t21
i13=$t31
i21=$t12
i22=$t22
i23=$t32
i31=$t13
i32=$t23
i33=$t33
echo "  Inverse matrix (i=t^T):"
echo "  | $i11 $i21 $i31 |"
echo "  | $i12 $i22 $i32 |"
echo "  | $i13 $i23 $i33 |"
echo
# i * basis
X1=`echo $i11 $i21 $i31 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
X2=`echo $i12 $i22 $i32 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
X3=`echo $i13 $i23 $i33 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`

Y1=`echo $i11 $i21 $i31 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Y2=`echo $i12 $i22 $i32 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Y3=`echo $i13 $i23 $i33 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`

Z1=`echo $i11 $i21 $i31 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Z2=`echo $i12 $i22 $i32 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Z3=`echo $i13 $i23 $i33 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
echo "  Basis (in rotation system):"
echo "      $X1      $Y1      $Z1"
echo "  x = $X2  y = $Y2  z = $Z2"
echo "      $X3      $Y3      $Z3"
echo

costheta=`echo $theta | awk '{ print cos($1*3.1415927/180.0) }'`
sintheta=`echo $theta | awk '{ print sin($1*3.1415927/180.0) }'`

x1=$X1
x2=$X2
x3=$X3
X1=$x1
X2=`echo $x2 $x3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
X3=`echo $x2 $x3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`

y1=$Y1
y2=$Y2
y3=$Y3
Y1=$y1
Y2=`echo $y2 $y3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
Y3=`echo $y2 $y3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`

z1=$Z1
z2=$Z2
z3=$Z3
Z1=$z1
Z2=`echo $z2 $z3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
Z3=`echo $z2 $z3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`

echo "  Transformed basis (in rotation system):"
echo "       $X1       $Y1       $Z1"
echo "  x' = $X2  y' = $Y2  z' = $Z2"
echo "       $X3       $Y3       $Z3"
echo

# t * basis
x1=`echo $t11 $t21 $t31 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
x2=`echo $t12 $t22 $t32 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
x3=`echo $t13 $t23 $t33 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
X1=$x1
X2=$x2
X3=$x3

y1=`echo $t11 $t21 $t31 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
y2=`echo $t12 $t22 $t32 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
y3=`echo $t13 $t23 $t33 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Y1=$y1
Y2=$y2
Y3=$y3

z1=`echo $t11 $t21 $t31 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
z2=`echo $t12 $t22 $t32 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
z3=`echo $t13 $t23 $t33 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
Z1=$z1
Z2=$z2
Z3=$z3


echo "  Transformed basis (cartesian coordinates):"
echo "       $X1       $Y1       $Z1"
echo "  x' = $X2  y' = $Y2  z' = $Z2"
echo "       $X3       $Y3       $Z3"
echo

cnt=0
tcnt=0

cx=1.0
cy=1.0
cz=0.8

# header
sed -n 1p POSCAR > $trg
echo " $lc" >> $trg
echo " $X1 $X2 $X3" >> $trg
echo " $Y1 $Y2 $Y3" >> $trg
echo " $Z1 $Z2 $Z3" >> $trg
echo " $sicnt $ccnt" >> $trg
echo "selective dynamics" >> $trg
echo "direct" >> $trg

tail -${total} POSCAR | \
	while read x y z fx fy fz; do
		Y=`echo $y $z $costheta $sintheta | \
			awk '{ print $3*$1+$4*$2 }'`
		Z=`echo $y $z $costheta $sintheta | \
			awk '{ print $3*$2-$4*$1 }'`
		echo " $x $Y $Z $fx $fy $fz" >> $trg
	done

echo "done"