#! /usr/bin/perl -w
#
#./build_DNA_base.pl nuc_base_r_theta_phi DNA_backbone_base_plane.pdb SEQ.dat > DNA.pdb
# 
#  nuc_base_r_theta_phi (input) : (r, theta, phi) coordinates for
#                                  for the base atoms of ade, thy, gua, cyt
#  DNA_backbone_base_plane.pdb (input) : a .pdb file containing backbonde 
#                             DNA strands and the plane definition (3 atoms) and
#                             NN (N1 for C, T; N3 for G, A) for each basis 
#                             plus any other protein molecules
#  SEQ.dat (input) : Sequence file
#                             SEQ A >ATATATATATAT ....
#                             SEQ B >TATATATATATA ....
#  DNA.pdb (output)         : The DNA with the threaded sequence
#
# Copyright (c) 2007 Peter Minary

#~~~~~~~~~~~~~~~~~~~~~~~~~~ define files ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  $nuc_base_file = $ARGV[0];
  $pdb_file      = $ARGV[1];
  $seq_file      = $ARGV[2];

  $pi = 3.1415926;  
  $inf = 1.e5;

#~~~~~~~~~~~~~~~~~~~~~~ read in nuc_base_file ~~~~~~~~~~~~~~~~~~~~~~~
  
  open(NUC_BASE,"$nuc_base_file");
  $nline = 0;
  while (<NUC_BASE>) {
         $nline++;
         chomp;
         $entry_base[$nline] = $_;
  }#endwhile#

  $iline = 0;
  $ntype = 4;
  for ($itype=1;$itype<=$ntype;$itype++) {
       $iline++;
#DEBUG
#    print "$entry_base[$iline]\n"; <STDIN>;
#DEBUG
       @list_entry_base = split(/ /,$entry_base[$iline]);
       $base_type = $list_entry_base[0];
       $ngroup    = $list_entry_base[1];
       for ($igroup=1;$igroup<=$ngroup;$igroup++) {
            $iline++;
            @list_entry_base = split(/ /,$entry_base[$iline]);
            $natm = $list_entry_base[3];      
            for ($iatm=1;$iatm<=$natm;$iatm++) {
                 $iline++;
                 @list_entry_base = split(/ /,$entry_base[$iline]);
                 $r{$base_type}[$iatm][$igroup]     = $list_entry_base[0];
                 $theta{$base_type}[$iatm][$igroup] = $list_entry_base[1];
                 $phi{$base_type}[$iatm][$igroup]   = $list_entry_base[2];
            }#endfor#
       }#endfor#
  }#endfor#

#DEBUG
#                 $base_type = 'A'; 
#                 $iatm   = 2;
#                 $igroup = 3;
#		 $r     =  $r{$base_type}[$iatm][$igroup];
#                 $theta =  $theta{$base_type}[$iatm][$igroup];
#                 $phi   =  $phi{$base_type}[$iatm][$igroup];
#                 print "$r, $theta, $phi\n";
#                <STDIN>;
#DEBUG

#~~~~~~~~~~~~~~~~~~~~~~ nuc_base definitions ~~~~~~~~~~~~~~~~~~~~~~~

# GUA
  @null       = ("0", "0");
  @gua_group1 = ("0", "N9");
  @gua_group2 = ("0", "C8", "H8", "N7", "C5", "C6", "O6", "N1", "H1", "C2", "N2", "N3", "C4");
  @gua_group3 = ("0", "H21", "H22");
  @gua_group  = (\@null, \@gua_group1, \@gua_group2, \@gua_group3);

  @gua_prim1 = ("0", "C3'", "C2'", "C1'", "N9");
  @gua_prim2 = ("0", "C2'", "C1'", "N9", "C8");
  @gua_prim3 = ("0", "N3", "C2", "N2", "H21");
  @gua_prim  = (\@null, \@gua_prim1, \@gua_prim2, \@gua_prim3);
  
  @gua_natom  = (0, 1, 12, 2 );
  $gua_ngroup =   3;

# ADE 
  @null       = ("0", "0");
  @ade_group1 = ("0", "N9");
  @ade_group2 = ("0", "C8", "H8", "N7", "C5", "C6", "N6", "N1", "C2", "H2", "N3", "C4");
  @ade_group3 = ("0", "H61", "H62");
  @ade_group  = (\@null, \@ade_group1, \@ade_group2, \@ade_group3);

  @ade_prim1 = ("0", "C3'", "C2'", "C1'", "N9");
  @ade_prim2 = ("0", "C2'", "C1'", "N9", "C8");
  @ade_prim3 = ("0", "C5", "C6", "N6", "H61");
  @ade_prim  = (\@null, \@ade_prim1, \@ade_prim2, \@ade_prim3);
  
  @ade_natom  = (0, 1, 11, 2 );
  $ade_ngroup =   3;

# CYT 
  @null       = ("0", "0");
  @cyt_group1 = ("0", "N1");
  @cyt_group2 = ("0", "C2", "O2", "N3", "C4", "N4", "C5", "H5", "C6", "H6");
  @cyt_group3 = ("0", "H42", "H41");
  @cyt_group  = (\@null, \@cyt_group1, \@cyt_group2, \@cyt_group3);

  @cyt_prim1 = ("0", "C3'", "C2'", "C1'", "N1");
  @cyt_prim2 = ("0", "C2'", "C1'", "N1", "C2");
  @cyt_prim3 = ("0", "N3", "C4", "N4", "H41");
  @cyt_prim  = (\@null, \@cyt_prim1, \@cyt_prim2, \@cyt_prim3);
  
  @cyt_natom  = (0, 1, 9, 2 );
  $cyt_ngroup =   3;

# THY 
  @null       = ("0", "0");
  @thy_group1 = ("0", "N1");
  @thy_group2 = ("0", "C2", "O2", "N3", "H3", "C4", "O4", "C5", "C5M", "C6", "H6");
  @thy_group3 = ("0", "H51", "H52", "H53");
  @thy_group  = (\@null, \@thy_group1, \@thy_group2, \@thy_group3);

  @thy_prim1 = ("0", "C3'", "C2'", "C1'", "N1");
  @thy_prim2 = ("0", "C2'", "C1'", "N1", "C2");
  @thy_prim3 = ("0", "C6", "C5", "C5M", "H51");
  @thy_prim  = (\@null, \@thy_prim1, \@thy_prim2, \@thy_prim3);
  
  @thy_natom  = (0, 1, 10, 3 );
  $thy_ngroup =   3;

# Base Plane

  @label_plane = ("0", "C2", "C4", "C6");
  $natom_plane = 3;

#~~~~~~~~~~~~~~~~~~~~~~~~~~ read in SEQ.dat ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  open(SEQ,"$seq_file");
  $ndseq  = 0;
  while (<SEQ>) {
    chomp;
    $entry = $_;
    if (substr($entry,0,3) eq 'SEQ') {
	 $ndseq++;
       @list_entry = split;
       $dseq_chain[$ndseq]  = $list_entry[1];
       $dseq[$ndseq]        = $list_entry[2];
       $dseq_length[$ndseq] = length($dseq[$ndseq]) - 1;
    }#endif#
  }#endwhile#

#~~~~~~~~~~~~~~~~~~~~~~~~~~ read in .pdb ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  open(PDB,"$pdb_file");
  $natom = 0;
  while (<PDB>) {
      chomp;
      $entry = $_;
      if (substr($entry,0,4) eq 'ATOM') {
         @list_entry = split;
         $atom = $list_entry[2];
         $altLoc = substr($entry,16,1);
         if ($altLoc eq ' ' || $altLoc eq 'A') {
             $natom++;
             $atom[$natom]       = substr($entry,12,4); $atom[$natom]    =~ tr/ //d;
             $resName[$natom]    = substr($entry,17,3); $resName[$natom] =~ tr/ //d;
             $chainID[$natom]    = substr($entry,21,1); $chainID[$natom] =~ tr/ //d;
             $resSeq[$natom]     = substr($entry,22,4); $resSeq[$natom]  =~ tr/ //d;
# NOUSE
#             $char_first[$natom] = substr($atom,0,1);
# NOUSE
             $x[$natom]     = substr($entry,30,8); $x[$natom] =~ tr/ //d;
             $y[$natom]     = substr($entry,38,8); $y[$natom] =~ tr/ //d;
             $z[$natom]     = substr($entry,46,8); $z[$natom] =~ tr/ //d;
             $entry[$natom] = $entry;
         }
#DEBUG
# print "$entry\n";
#DEBUG
      }
#NOUSE
#      if (substr($entry,0,3) eq 'END' || substr($entry,0,3) eq 'TER') {
#          last;
#      }
#NOUSE
  }

#~~~~~~~~~~~~~~~~~~~~~~~~~~ separate atoms into residues ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  $nres   = 1;
  $resSeq = $resSeq[1];
  $natom_res[$nres] = 0;
  for ($iatom=1;$iatom<=$natom;$iatom++) {
       if ($resSeq[$iatom] ne $resSeq) {
           $resSeq = $resSeq[$iatom];
           $nres++;
           $natom_res[$nres]++;
           $ind_atom[$nres][$natom_res[$nres]] = $iatom;
       } else {
           $natom_res[$nres]++;
           $ind_atom[$nres][$natom_res[$nres]] = $iatom;  
       }
  }

#~~~~~~~~~~~~~~~~~~~~~~~~~ separate residues into molecules based on chainID ~~~~~~~~~~~~~~~~~~
  $nmol = 1;
  $chainID = $chainID[1];
  $nres_mol[$nmol] = 0;
  $ind_chainID[$nmol] = $chainID;
  for ($ires=1;$ires<=$nres;$ires++) {
       $atom1 = $ind_atom[$ires][1];
#DEBUG
#  print "$chainID $chainID[$atom1]\n";
#  print "$entry[$atom1]\n";
#DEBUG
       if ($chainID[$atom1] ne $chainID) {
           $chainID = $chainID[$atom1];   
           $nmol++;
           $nres_mol[$nmol]++;
           $ind_chainID[$nmol] = $chainID;
           $ind_res[$nmol][$nres_mol[$nmol]] = $ires;
       } else {
           $nres_mol[$nmol]++;
           $ind_res[$nmol][$nres_mol[$nmol]] = $ires;
       }
  }

#~~~~~~~~~~~~~~~~~~~~~~~~~ identify molecules as DNA or protein  ~~~~~~~~~~~~~~~~~~

   $nDNA = 0;
#   $nprot = 0;
   for ($imol=1;$imol<=$nmol;$imol++) {
        $ires    = $ind_res[$imol][1];
        $iatm    = $ind_atom[$ires][1];
        $resName = $resName[$iatm];
        $length  = length($resName);
        if ($length == 3) { $ityp_mol[$imol] = 'protein'; }
        if ($length  < 3) { $ityp_mol[$imol] = 'DNA'; }
#DEBUG
#        print "$ityp_mol[$imol] \n";
#DEBUG
        if ($ityp_mol[$imol] eq 'DNA') {
            $nDNA++;
            $ind_DNA[$nDNA] = $imol;
        }
#        if ($ityp_mol[$imol] eq 'protein') {
#            $nprot++;
#            $ind_prot[$nprot] = $imol;
#        }
   }
   if ($nDNA != 2)
      { die "The number of DNA chains is not 2";}
   if ($nres_mol[$ind_DNA[1]] != $nres_mol[$ind_DNA[2]])
      { die "DNAs do not have the same # of residues\n";}


#~~~~~~~~~~~~~~~~~~~~~~~~~ build on DNA backbone ~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ write out DNA ~~~~~~~~~~~~~~~~~~~~~~

   for ($imol=1;$imol<=$nmol;$imol++) {
       if ($ityp_mol[$imol] ne 'DNA') {next;}
        $design = 0;
	for ($idseq=1;$idseq<=$ndseq;$idseq++) {
#DEBUG
#  print "$dseq_chain[$idseq]\n"; <STDIN>;
#DEBUG
             if ($dseq_chain[$idseq] eq  $ind_chainID[$imol]) {$design = 1; $iidseq = $idseq; last;}
        }#endfor#
        if ($design == 1) {
 	  if ( $dseq_length[$iidseq] != $nres_mol[$imol] ) { die "design sequence length mismatch"; }
        }#endif#
 # $dseq[$iidseq];
        $nres = $nres_mol[$imol];
        for ($ires=1;$ires<=$nres;$ires++) {
	     $iires = $ind_res[$imol][$ires];
             $natm  = $natom_res[$iires];
             for ($iatm=1;$iatm<=$natm;$iatm++) {
  		  $iiatm = $ind_atom[$iires][$iatm];
                  $res_type = substr($dseq[$iidseq],$ires,1);
                  substr($entry[$iiatm],17,3) = $res_type."  ";
#DEBUG
#		  print "res_type $res_type\n";
#DEBUG
                  if (($atom[$iiatm] ne "PLN") && 
                      ($atom[$iiatm] ne "NN")) {
                      print "$entry[$iiatm]\n";
                  }#endif#
		  if ($atom[$iiatm] eq "NN") {$ind_NN = $iiatm;}
             }#endfor-iatm#

             if ($design == 1) {
                    $res_type = substr($dseq[$iidseq],$ires,1);

                    $natom_plane = 0;
                    for ($iatm=1;$iatm<=$natm;$iatm++) {
			 $iiatm = $ind_atom[$iires][$iatm];
                         if ($atom[$iiatm] eq "PLN") {
                             $natom_plane++;
                             $atom_plane[$natom_plane] = $iiatm;
                         }#endif# 
                    }#endfor#                     

#DEBUG
# print "$iires residie is being designed, res type $res_type\n"; <STDIN>;
#DEBUG
                if ($res_type eq "G"){@prim=@gua_prim;@group=@gua_group;$ngroup=$gua_ngroup;
                                      @natm_group=@gua_natom;}
                if ($res_type eq "A"){@prim=@ade_prim;@group=@ade_group;$ngroup=$ade_ngroup;
                                      @natm_group=@ade_natom;}
                if ($res_type eq "C"){@prim=@cyt_prim;@group=@cyt_group;$ngroup=$cyt_ngroup;
                                      @natm_group=@cyt_natom;}
                if ($res_type eq "T"){@prim=@thy_prim;@group=@thy_group;$ngroup=$thy_ngroup;
                                      @natm_group=@thy_natom;}
                for ($igroup=1;$igroup<=$ngroup;$igroup++) {
                     $ngatm = $natm_group[$igroup];
                     for ($iprim=1;$iprim<=3;$iprim++) {
                          $label = $prim[$igroup][$iprim];
                          $natm  = $natom_res[$iires];
                          $atom_found = 0;
                          for ($iatm=1;$iatm<=$natm;$iatm++) {
    		               $iiatm = $ind_atom[$iires][$iatm];
#DEBUG                        
# print "atom #$atom[$iiatm]# label #$label# iprim $iprim igroup $igroup\n"; <STDIN>;
#DEBUG
                               if ($atom[$iiatm] eq $label) {
#DEBUG
# print "atom #$atom[$iiatm]# label #$label#\n"; <STDIN>;
#DEBUG
                                   $rxp[$iprim] = $x[$iiatm];
                                   $ryp[$iprim] = $y[$iiatm];
                                   $rzp[$iprim] = $z[$iiatm];
#                                   $ind_atm_prim[$iprim][$igroup] = $iiatm; $atom_found = 1;
                               }#endif#
			       if ($atom_found==0) {
                                   for ($jgroup=1;$jgroup<$igroup;$jgroup++) {
				        for ($jatm=1;$jatm<=$natm_group[$jgroup];$jatm++) {
					     if ($label eq $group[$jgroup][$jatm]) {
                                                 $rxp[$iprim] = $x_group[$jgroup][$jatm];
                                                 $ryp[$iprim] = $y_group[$jgroup][$jatm];
                                                 $rzp[$iprim] = $z_group[$jgroup][$jatm];
                                             }#endif#
					}#endfor#
                                   }#endfor#
                               }#endif#
                          }#endfor#
                     }#endfor#
              
                     # define r1, r2, r3
                     $r1x = $rxp[1];  $r2x = $rxp[2];  $r3x = $rxp[3];
                     $r1y = $ryp[1];  $r2y = $ryp[2];  $r3y = $ryp[3];
		     $r1z = $rzp[1];  $r2z = $rzp[2];  $r3z = $rzp[3];

                     # find plane atoms in group
                     $ap_found = 0;
		     for ($iap=1;$iap<=$natom_plane;$iap++) {
                         for ($igatm=1;$igatm<=$ngatm;$igatm++) {
                          if ($label_plane[$iap] eq $group[$igroup][$igatm]) {
                              $ind_atom_plane[$iap] = $igatm;
                              $ap_found = 1;
                          }#endif#
                         }#endfor# 
                     }#endfor#
#DEBUG
#		      print "$ap_found \n"; <STDIN>;
#DEBUG      
                         

                     # determine phi_shift
                     if ($ap_found == 1) {
			 $phi_shift = 0;
                         $nshift = 100;
                         $error_min = $inf; 
                         for ($ishift=-$nshift;$ishift<=$nshift;$ishift++) {
                           $phi_shift = 2.0*$ishift/$nshift*$pi;  
                           $error_dist = 0;
			   for ($iap=1;$iap<=$natom_plane;$iap++) {
			  	  $igatm = $ind_atom_plane[$iap];
#DEBUG	#DEBUG
#			    print "$igatm \n"; <STDIN>;
#DEBUG	#DEBUG
			  	  $r     = $r{$res_type}[$igatm][$igroup]; 
			  	  $theta = $theta{$res_type}[$igatm][$igroup];
			  	  $phi   = $phi{$res_type}[$igatm][$igroup];
                                  $phi  += $phi_shift;  
                                
                                &build_r4; 
	
                                $xagp[$iap] = $r4x;
                                $yagp[$iap] = $r4y;
                                $zagp[$iap] = $r4z;
 	 
                                $xap[$iap] = $x[$atom_plane[$iap]]; 
                                $yap[$iap] = $y[$atom_plane[$iap]]; 
                                $zap[$iap] = $z[$atom_plane[$iap]];

                              $error_dist +=sqrt (
                                             ($xagp[$iap]-$xap[$iap])*($xagp[$iap]-$xap[$iap]) 
                                            +($yagp[$iap]-$yap[$iap])*($yagp[$iap]-$yap[$iap]) 
                                            +($zagp[$iap]-$zap[$iap])*($zagp[$iap]-$zap[$iap])
                                                  );
			   }#endfor#
			   $x21 = $xap[2] - $xap[1];
			   $y21 = $yap[2] - $yap[1];
			   $z21 = $zap[2] - $zap[1];
	
			   $x31 = $xap[3] - $xap[1];
			   $y31 = $yap[3] - $yap[1];
			   $z31 = $zap[3] - $zap[1];
                           
  	  		   $xn  = $y21*$z31 - $z21*$y31;
  	  		   $yn  = $z21*$x31 - $x21*$z31;
 	 		   $zn  = $x21*$y31 - $y21*$x31;
                           $n2  = $xn*$xn + $yn*$yn + $zn*$zn;                     
	
                           $error_plane = 0; 
                           for ($iap=1;$iap<=$natom_plane;$iap++) {
                                $scl_prod = $xn*($xap[1]-$xagp[$iap]) 
                                          + $yn*($yap[1]-$yagp[$iap])
                                          + $zn*($zap[1]-$zagp[$iap]);
                                $error_plane += sqrt($scl_prod*$scl_prod)/sqrt($n2);
                           }#endfor#
			   $error = 10.0*$error_plane + $error_dist;
                           if (($error) < $error_min) {
                               $error_min = $error;
                               $phi_shift_correct = $phi_shift;
                           }#endif# 
#DEBUG
#			   print "$error\n";
#			   <STDIN>;
#DEBUG
		       }#ishift#
#DEBUG
#			   print "error $error_min\n";
#			   <STDIN>;
#DEBUG
		     }#endif-ap_found#


                     for ($igatm=1;$igatm<=$ngatm;$igatm++) {
			  $r     = $r{$res_type}[$igatm][$igroup]; 
			  $theta = $theta{$res_type}[$igatm][$igroup];  
			  $phi   = $phi{$res_type}[$igatm][$igroup];  
                          
                          if ($ap_found) {
                              $phi  += $phi_shift_correct;
                          }#endif#        

                          &build_r4;

                          if ($igroup==1) {
                              $r4x = $x[$ind_NN];
                              $r4y = $y[$ind_NN];
                              $r4z = $z[$ind_NN];
                          }#endif#

                          $x_group[$igroup][$igatm] = $r4x;
                          $y_group[$igroup][$igatm] = $r4y;
                          $z_group[$igroup][$igatm] = $r4z;
#DEBUG
#   print "iiatm $iiatm\n"; <STDIN>;
#DEBUG
                          $entry_new = $entry[$iiatm];
                          $atom_new  = sprintf("%-4s",$group[$igroup][$igatm]);
#DEBUG
#   print "atom_new #$atom_new#\n"; <STDIN>;
#DEBUG
                          substr($entry_new,12,4) = $atom_new;
                          $x_field = sprintf("%8.5g",$r4x);
                          $y_field = sprintf("%8.5g",$r4y);
                          $z_field = sprintf("%8.5g",$r4z);
                          substr($entry_new,30,8) = $x_field;
                          substr($entry_new,38,8) = $y_field;
                          substr($entry_new,46,8) = $z_field;
#                          $length = length($entry_new);
#                          substr($entry_new,$length-6,6) = "     ";
#DEBUG
#   print "$entry[$iiatm]\n"; <STDIN>;
#DEBUG
                          print "$entry_new\n";
                     }#endfor-igatm#
                }#endfor-igroup#

             }#endif-design#  


        }#endfor- ires#
   }#endfor - imol#


#~~~~~~~~~~~~~~~~~~~~~~ write everything but DNA ~~~~~~~~~~~~~~~~~~~
#   print "CBLC >CD\n";
   for ($imol=1;$imol<=$nmol;$imol++) {
        if ($ityp_mol[$imol] eq 'DNA') {next;}
        $nres = $nres_mol[$imol];
        for ($ires=1;$ires<=$nres;$ires++) {
	     $iires = $ind_res[$imol][$ires];
             $natm  = $natom_res[$iires];
             for ($iatm=1;$iatm<=$natm;$iatm++) {
  		  $iiatm = $ind_atom[$iires][$iatm];
                  print "$entry[$iiatm]\n";
	     }#endfor#
	}#endfor#
   }#endfor#


#~~~~~~~~~~~~~~~~~~~~~~~~ subroutines ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 sub build_r4 
 {
     # builds r4 based on (r1, r2, r3, r, theta, phi)
     # defined in the main body of the program 

     # define e1 = (r3-r2)/||r3-r2||, e2 = (r1-r2)/||r1-r2||
     $e1x = $r3x - $r2x; $e1y = $r3y - $r2y; $e1z = $r3z - $r2z; 
     $e2x = $r1x - $r2x; $e2y = $r1y - $r2y; $e2z = $r1z - $r2z;

     $e1_2 = $e1x*$e1x + $e1y*$e1y + $e1z*$e1z;
     $e1   = sqrt($e1_2);
     $e1x  /= $e1; $e1y /= $e1; $e1z /= $e1;                 

     $e2_2 = $e2x*$e2x + $e2y*$e2y + $e2z*$e2z;
     $e2   = sqrt($e2_2);
     $e2x  /= $e2; $e2y /= $e2; $e2z /= $e2;                 

     # define ez = e1
     $ez_x = $e1x; $ez_y = $e1y; $ez_z = $e1z;
     
     # define ey = (e1 x e2) / ||e1 x e2||
     $ey_x  = $e1y*$e2z - $e1z*$e2y;
     $ey_y  = $e1z*$e2x - $e1x*$e2z;
     $ey_z  = $e1x*$e2y - $e1y*$e2x;

     $ey_2 = $ey_x*$ey_x + $ey_y*$ey_y + $ey_z*$ey_z;
     $ey   = sqrt($ey_2);
     $ey_x  /= $ey; $ey_y /= $ey; $ey_z /= $ey;                 
     
     # define ex = (ey x ez) / ||ey x ez||
     $ex_x  = $ey_y*$ez_z - $ey_z*$ez_y;
     $ex_y  = $ey_z*$ez_x - $ey_x*$ez_z;
     $ex_z  = $ey_x*$ez_y - $ey_y*$ez_x;

     $ex_2 = $ex_x*$ex_x + $ex_y*$ex_y + $ex_z*$ex_z;
     $ex   = sqrt($ex_2);
     $ex_x  /= $ex; $ex_y /= $ex; $ex_z /= $ex;                 

     #reconstruct in reference frame (ex, ey, ez)

     $rx    = $r*sin($theta)*cos($phi);
     $ry    = $r*sin($theta)*sin($phi);
     $rz    = $r*cos($theta);
     
     #rotate from reference frame to real frame
     $rrx   = $ex_x*$rx + $ey_x*$ry + $ez_x*$rz;
     $rry   = $ex_y*$rx + $ey_y*$ry + $ez_y*$rz;
     $rrz   = $ex_z*$rx + $ey_z*$ry + $ez_z*$rz;
     
     #place vector by adding it to $r3
     $r4x = $r3x + $rrx; 
     $r4y = $r3y + $rry; 
     $r4z = $r3z + $rrz;

 }