#!perl
use strict;
use MaterialsScript qw(:all);
use Math::Trig;

my $doc = $Documents{"Layerhalf.xtd"};
my $setA = "Si";
my $setB = "wat";

# get the trajectory data
my $trj = $doc->Trajectory;
my $numFrames = $trj->NumFrames;

# create a study table for the results

my $std = Documents->New($doc->Name.".std");
$std->ColumnHeading(0) = "Angle (deg)";
$std->ColumnHeading(1) = "Distance (angstrom)";

my $setMolecules = $doc->UnitCell->Sets($setB)->Molecules;

# run over all frames in the trajectory

my $irow = 0;
my $setAtoms = $doc->UnitCell->Sets($setA)->Atoms;
my $plane = $doc->UnitCell->CreateBestFitPlane($setAtoms);
my $bestFitPlane = $plane->BestFitPlane;

for(my $frame = 201; $frame <= 300; $frame++)

{
   $trj->CurrentFrame = $frame;

   # get the normal vector of the plane
   my $normal = $plane->BestFitPlane->UnitNormalVector;

   foreach my $setMolecule(@$setMolecules)
   {  
     # z position of the molecule
     my $distance = $setMolecule->Center->Z;

     # get dipole moment
     my $dipoleMoment = $setMolecule->DipoleMoment;
     my $direction=Normalize($dipoleMoment);

     # calculate the angle between the dipole and the surface normal  
     my $cosAngle = InProduct($direction,$normal);
     my $angle =180/pi*acos($cosAngle);

     # store the data in the study table

     $std->Cell($irow,0) = $angle;       
     $std->Cell($irow++,1) = $distance;
   }
}
$doc->Discard; # discard any changes made

sub InProduct{
my($a, $b) = @_;
return $a->X*$b->X+$a->Y*$b->Y+$a->Z*$b->Z
}