SASA

#!/usr/bin/perl

# SASA using Shrake-Rupley algorithm
# calculate Solvent Accessible Surface Area
# SASA
# VR 2013

use Math::Trig;
use Data::Dumper;


#H2O radius:
$prad=1.4;

#tasku sk. aplink atomus (daugiau tiksliau: ~simtai):
$M=500;

#atoms example 0,1,...  (x,y,z,radius):
@N[0] = [1,2,3,4];
@N[1] = [0,1,2,3];


$i = 0;
$irad = 0;

while(@N[$i]) {
 print "skersmuo: ", $N[$i][3];
 $irad=$N[$i][3] + $prad;
 while ($k<$M) {
  $u=rand();
  $v=rand();
  $theta=2*pi*$u;
  $phi=acos(2*$v-1);
  $x=cos($theta)*sin($phi);
  $y=sin($theta)*sin($phi);
  $z=cos($phi);
  #sukuriamas taskas:
  $point[0]=$N[$i][0]+$x*$irad;
  $point[1]=$N[$i][1]+$y*$irad;
  $point[2]=$N[$i][2]+$z*$irad;
  push @{$pts[$i][$k]}, @point;
 $k++;
 }
$i++;
$k = 0;
}

#tasku visuma:
#print Dumper \@pts;

#eit per atomus ir istrinti taskus kurie kito atomo kelyje:
$i = 0;
$k = 0;

while(@N[$i]) {
 $irad=$N[$i][3] + $prad;
 $Mp=$M;
 while ($k<$M) {
  #pt(:)=pts(i,k,:)
  $fail=0;
  $j=$i+1;
  while(@N[$j]) {
   $jrad=$N[$j][3] + $prad;
   #jcoord=(/x,y,z/) of atom j
   $r=sqrt(( ($pts[$i][$k][0]-$N[$j][0])+($pts[$i][$k][1]-$N[$j][1])+($pts[$i][$k][2]-$N[$j][2]) )**2);
   if($r <= $jrad) {
    $fail=1;
   }
  $j++;
  }
  if($fail) {
   $Mp--;
  }
 $k++;
 }
 #SASA atomu i:
 $sasa[$i]=4*pi*$irad*$irad*$Mp/$M;
$i++;
$k = 0;
}

print "sasa: \n";
print Dumper \@sasa;