Beefy Boxes and Bandwidth Generously Provided by pair Networks
Your skill will accomplish
what the force of many cannot
 
PerlMonks  

mwah's scratchpad

by mwah (Hermit)
on Sep 19, 2007 at 16:10 UTC ( #639921=scratchpad: print w/ replies, xml ) 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 my @tmpl = ( [0,0,0,'HEAD'], [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'], [9,0,0,'HEAD'] ); 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, $increment, deg2rad($increment); 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__



Log In?
Username:
Password:

What's my password?
Create A New User
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others contemplating the Monastery: (6)
As of 2015-03-05 03:12 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    When putting a smiley right before a closing parenthesis, do you:









    Results (135 votes), past polls