Beefy Boxes and Bandwidth Generously Provided by pair Networks
Clear questions and runnable code
get the best and fastest answer
 
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 all is quiet...

How do I use this? | Other CB clients
Other Users?
Others meditating upon the Monastery: (2)
As of 2017-12-16 21:42 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?
    What programming language do you hate the most?




















    Results (459 votes). Check out past polls.

    Notices?