Welcome to the Monastery PerlMonks

by mwah (Hermit)
 on Sep 19, 2007 at 16:10 UTC Need Help??

This code is related to a question (Doing rigid body rotations in Perl) I asked some time ago:

```#
# example generate a perfect spiral (helix) of (10-Atom)-rods
#
use Math::Trig;

my \$axis = [ 0, 0, 1 ]; # z axis

[1,0,0,'SEG'],
[2,0,0,'SEG'],
[3,0,0,'SEG'],
[4,0,0,'SEG'],
[5,0,0,'SEG'],
[6,0,0,'SEG'],
[7,0,0,'SEG'],
[8,0,0,'SEG'],
);

my (\$nMol, \$nPeriods) = (360, 12);  # 360 Molecules make 12 helix tur
+ns
my \$increment = \$nPeriods*360.*(1./\$nMol); # put that into output

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
my @ensemble;

for my \$id (0..\$nMol-1) {
my \$w = deg2rad      (\$nPeriods * 360. * (\$id/\$nMol));
my \$m = setrotmatrix ( \$w, \$axis );
my \$p = deepcopy_body( \@tmpl );
center_body   ( \$p );
translate_body( \$p => [0,0,\$id] );
rotate_body   ( \$p => \$m );
push @ensemble, \$p;
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

{ # do another rotation of the whole thing
my \$mat = setrotmatrix( deg2rad(36), normalize([ 1, 1, 0 ]) );
for my \$en (@ensemble) { rotate_body( \$en, \$mat ) }
}

my \$fn = sprintf "ihx-b%02d-n%04d-p%d.kar", 0+@tmpl, \$nMol, \$nPeriods
+;
open my \$fh, '>',\$fn or die "can't open \$!";

printf \$fh "%d\t # %f/%f deg/rad increment\n",
@tmpl * @ensemble,

my(\$molecule, \$idmol, \$boffs) = (0,1,1);

# coordinates
for \$molecule (@ensemble) {
for my \$coord (@\$molecule) { print \$fh "@\$coord \$idmol 1\n" }
\$idmol++;
}
# bonds
for \$molecule (@ensemble) {
for my \$n (0..@\$molecule-2) { print \$fh \$boffs+\$n, ' ', \$boffs+\$n+
+1, "\n" }
\$boffs += @tmpl;
}
print \$fh "0 0\n";
close \$fh;

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

# sub deg2rad { pi * \$_[0] / 180. } => already defined in Math::Trig
# sub rad2deg { 180. * \$_[0] / pi } => already defined in Math::Trig

sub setrotmatrix
{
my (\$theta, \$axis) = @_;
my (\$x, \$y, \$z) = @\$axis;
my \$c = cos \$theta;
my \$s = sin \$theta,
my \$t = 1. - \$c;
return [
[\$t * \$x * \$x + \$c,       \$t * \$x * \$y + \$s * \$z,  \$t * \$x * \$z - \$s
+ * \$y],
[\$t * \$x * \$y - \$s * \$z,  \$t * \$y * \$y + \$c,       \$t * \$y * \$z + \$s
+ * \$x],
[\$t * \$x * \$z + \$s * \$y,  \$t * \$y * \$z - \$s * \$x,  \$t * \$z * \$z + \$c
+     ]
]
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub rotate
{
my (\$v, \$m) = @_;
my (\$x, \$y, \$z) = @\$v;
\$v->[0] =  \$m->[0][0]*\$x + \$m->[0][1]*\$y + \$m->[0][2]*\$z ;
\$v->[1] =  \$m->[1][0]*\$x + \$m->[1][1]*\$y + \$m->[1][2]*\$z ;
\$v->[2] =  \$m->[2][0]*\$x + \$m->[2][1]*\$y + \$m->[2][2]*\$z ;
return \$v;
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub veclength2
{
my (\$v) = @_;
my (\$x, \$y, \$z) = @\$v;
return \$x**2 + \$y**2 + \$z**2
}

sub veclength {
my (\$v) = @_;
return sqrt( veclength2(\$v) )
}

sub normalize
{
my (\$v) = @_;
my \$len = veclength(\$v);
\$len = (\$len < 1e-10) ? 0. : 1./\$len;
\$v->[0] *= \$len, \$v->[1] *= \$len, \$v->[2] *= \$len;
return \$v;
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub deepcopy_body
{
my (\$body) = @_;
my \$copy = [];
for my \$p ( @\$body ) { push @\$copy, [ @\$p ] }
return \$copy
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub translate_body
{
my (\$body, \$disp) = @_;
for my \$p (@\$body) {
\$p->[0] += \$disp->[0], \$p->[1] += \$disp->[1], \$p->[2] += \$disp->[2
+];
}
return \$body
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub center_body
{
my (\$body) = @_;
my \$cnt = @\$body;
my (\$pcoord, @com) = 0, (0,0,0);

for \$pcoord (@\$body) {
\$com[0] += \$pcoord->[0],
\$com[1] += \$pcoord->[1],
\$com[2] += \$pcoord->[2]
}
\$com[0] /= \$cnt,
\$com[1] /= \$cnt,
\$com[2] /= \$cnt;

for \$pcoord (@\$body) {
\$pcoord->[0] -= \$com[0],
\$pcoord->[1] -= \$com[1],
\$pcoord->[2] -= \$com[2]
}
return \$body
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #

sub rotate_body
{
my (\$body, \$m) = @_;
for my \$v (@\$body) {
rotate(\$v, \$m);
}
return \$body
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+- - - - #
__END__

Create A New User
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others surveying the Monastery: (4)
As of 2021-04-18 09:21 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?

No recent polls found

Notices?