#!/usr/bin/perl -w

# calculate the bond angle alpha and dihedral angle tor
# of the backbone of a amino acids
# input is one pdbfile

if(!@ARGV){
	printf "ccAngle.pl: calculate the dihedral angles and bond angles of Ca-Ca seudo bonds.\n";
	printf "example: ccAngle.pl 1a32.pdb\n";
	exit(0);
}

sub acos {

   my($x) = @_;
   my $ret = atan2(sqrt(1 - $x**2), $x);
   return $ret;
}

sub distance {
    my ($p1,$p2) = @_;
    $x=@$p1[0]-@$p2[0];
    $y=@$p1[1]-@$p2[1];
    $z=@$p1[2]-@$p2[2];
    return sqrt($x*$x+$y*$y+$z*$z);
};

sub angle{
    my($p1,$p2,$p3) =@_;
    $d12 = distance(\@$p1,\@$p2);
    $d23 = distance(\@$p2,\@$p3);
    $x12=@$p1[0]-@$p2[0]; $y12=@$p1[1]-@$p2[1]; $z12=@$p1[2]-@$p2[2];
    $x23=@$p2[0]-@$p3[0]; $y23=@$p2[1]-@$p3[1]; $z23=@$p2[2]-@$p3[2];
    $dot=$x12*$x23 + $y12*$y23 + $z12*$z23;
    $answer = 180 * acos($dot/($d12*$d23))/3.1415926;
    return (180-$answer);
};

sub torsion{
    my($p1,$p2,$p3,$p4) =@_;
    $x1=@$p2[0]-@$p1[0]; $y1= @$p2[1]-@$p1[1]; $z1 = @$p2[2]-@$p1[2];
    $x2=@$p2[0]-@$p3[0]; $y2= @$p2[1]-@$p3[1]; $z2 = @$p2[2]-@$p3[2];
    $x3=@$p3[0]-@$p4[0]; $y3=@$p3[1]-@$p4[1]; $z3 =@$p3[2]-@$p4[2];
    $axbbxc = ($y1*$z2-$z1*$y2)*($y2*$z3-$z2*$y3)+($x2*$z1-$x1*$z2)*($x3*$z2-$x2*$z3)+($x1*$y2-$x2*$y1)*($x2*$y3-$x3*$y2);
    $vaxb = ($y1*$z2-$z1*$y2)*($y1*$z2-$z1*$y2)+($x2*$z1-$x1*$z2)*($x2*$z1-$x1*$z2)+($x1*$y2-$x2*$y1)*($x1*$y2-$x2*$y1);
$vbxc = ($y2*$z3-$z2*$y3)*($y2*$z3-$z2*$y3)+($x3*$z2-$x2*$z3)*($x3*$z2-$x2*$z3)+($x2*$y3-$x3*$y2)*($x2*$y3-$x3*$y2);
    $answer = $axbbxc/sqrt($vaxb*$vbxc);
    $answer = 180*acos($answer)/3.14;
    $sign = $x1*($y2*$z3-$y3*$z2)-$y1*($x2*$z3-$x3*$z2)+$z1*($x2*$y3-$x3*$y2);

    if($sign>=0){
    	if ($answer>0) {$answer=$answer;}
    	if ($answer<0) {$answer=-$answer;}
    }
    if($sign<0){
    	if ($answer>=0) {$answer=-$answer;}
    	if ($answer<0) {$answer=$answer;}
    }
    if($answer>=0) {
	return 180-$answer;
    }
    else { 
	return -180-$answer;
    }
};

if($ARGV[0] =~ /^(\S+).ab.pdb$/){
	$outname=$1;
}
elsif($ARGV[0] =~ /^(\S+).pdb$/){
	$outname=$1;
}
else {
	$outname =$ARGV[0];
}

$outname =~ tr/A-Z/a-z/;
if(-e "$outname.ab.pdb"){
    	open PDBFILE, "$outname.ab.pdb";
}
elsif(-e "$outname.sc.pdb") {
    	open PDBFILE, "$outname.sc.pdb";
}
elsif(-e "$outname.ca.pdb") {
    	open PDBFILE, "$outname.ca.pdb";
}
elsif(-e "$outname.pdb"){
    	open PDBFILE, "$outname.pdb";
}
else { 
	printf "cant open $outname.pdb.\n";
}

open OUTFILE,">$outname.ang";

$pos=0;
$res_co=[['0.0','0.0','0.0'],['0.0','0.0','0.0'],['0.0','0.0','0.0'],['0.0','0.0','0.0']];
@res_name={' ',' ',' ',' '};

$round=0;

for($i=0;$i<3;++$i) {
    	$line=<PDBFILE>;
	if($line !~ /^ATOM/) {
    		next;
	}
	if(substr($line,16,4) =~ /[BCD23][A-Z][A-Z][A-Z]/) {
		next;
	}
	$atomName= substr($line, 12,4); 
	$atomName=~ s/^\s+|\s+$//g;
	if($atomName ne "CA") {
		$i--;
		next;
	}

	$res_co->[$round][0]=substr($line,30,8);$res_co->[$round][0]=~ s/^\s+|\s+$//g;
	$res_co->[$round][1]=substr($line,38,8);$res_co->[$round][1]=~ s/^\s+|\s+$//g;
	$res_co->[$round][2]=substr($line,46,8);$res_co->[$round][2]=~ s/^\s+|\s+$//g;
	$res_name[$round] = substr($line,17,3); $res_name[$round] =~ s/^\s+|\s+$//g;
	$resNum[$round]  = substr($line, 22,4); $resNum[$round]  =~ s/^\s+|\s+$//g;
	++$pos;
	$round=($round+1)%4;

} # end of for($i..

while(defined($line=<PDBFILE>)){  
    	if($line !~ /^ATOM/) {
		next;
}

if(substr($line,16,4) =~ /[BCD23][A-Z][A-Z][A-Z]/) {
    	next;
}
$atomName= substr($line, 12,4); $atomName=~ s/^\s+|\s+$//g;
if($atomName ne "CA") {
	next;
}

$res_co->[$round][0]=substr($line,30,8);$res_co->[$round][0]=~ s/^\s+|\s+$//g;
$res_co->[$round][1]=substr($line,38,8);$res_co->[$round][1]=~ s/^\s+|\s+$//g;
$res_co->[$round][2]=substr($line,46,8);$res_co->[$round][2]=~ s/^\s+|\s+$//g;
$res_name[$round] = substr($line,17,3);$res_name[$round] =~ s/^\s+|\s+$//g;
$resNum[$round]  = substr($line, 22,4); $resNum[$round]  =~ s/^\s+|\s+$//g;
++$pos;
$round=($round+1)%4;
$ang= angle(\@{$res_co->[($round-3)%4]},\@{$res_co->[($round-2)%4]},\@{$res_co->[($round-1)%4]});
$dis = distance(\@{$res_co->[($round-3)%4]},\@{$res_co->[($round-2)%4]});
$tor = torsion(\@{$res_co->[($round-4)%4]},\@{$res_co->[($round-3)%4]},\@{$res_co->[($round-2)%4]},\@{$res_co->[($round-1)%4]});
$res_name[$round-2] =~ tr/[a-z]/[A-Z]/;
printf OUTFILE "%5s %4d %6.2f %6.2f %5.3f\n",$res_name[$round-2],$resNum[$round-2],$ang,$tor,$dis;

} # end of while
