Skip to content
Snippets Groups Projects
Commit 0166305c authored by slevy's avatar slevy
Browse files

Allow multiline input. 

Add precession.
Add determinant.  Make sure that &basis() returns a
positively-oriented matrix!
More precision in dms output.
Fix bug: qrotbtwn() had the wrong sign!
parent 3472ce48
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 140 additions and 18 deletions
scripts/tfm.pl
+ 140
− 18
View file @ 0166305c
#! /usr/bin/perl5
# Perl geometry calculator, with special assistance for astronomical coordinates.
# By Stuart Levy, NCSA, University of Illinois, 1997-2001. slevy@ncsa.uiuc.edu.
$pi = 3.14159265358979323846;
$choplimit = 2e-14;
&init_eq2gal;
......@@ -17,7 +20,7 @@ Here "T" is a 4x4 matrix as list of 16 numbers
tfm(ax,ay,ax, angle, cx,cy,cz) 4x4 rot about axis, fixing center cx,cy,cz
transpose( T ) NxN matrix transpose
tmul( T1, T2 ) => T1*T2 4x4 (or 3x3) matrix product
eucinv( T ) => Tinverse 4x4 inverse (assuming T Euclidean rot/trans/scale)
eucinv( T ) => Tinverse 4x4 inverse (assuming T Euclidean rot/trans/uniform scale)
hls2rgb(h,l,s) => (r,g,b) color conversion
svmul( s, v ) => s*v scalar * vector
vmmul( v4, T ) => v' 4-vector * 4x4 matrix => 4-vector
......@@ -32,17 +35,20 @@ Here "T" is a 4x4 matrix as list of 16 numbers
quat2t( q ) => T quaternion to 4x4 matrix T
quatmul(qa, qb) => qa*qb quaternion multiplication
qrotbtwn(v3a, v3b) => q quaternion which rotates 3-vector va into vb
lookat(from3,to3,up3,roll) construct w2c matrix.
aer2t(Ry,Rx,Rz) => T and t2aer(T) => Ry,Rx,Rz
vd2tfm(x,y,z,Rx,Ry,Rz) and tfm2vd(T) 4x4 matrix <=> tx ty tz rx ry rz
eq2dms(v3) => "hh:mm.m +dd:mm:ss dist"
lookat(from3,to3,up3,roll) construct 4x4 world->camera w2c matrix (pw * w2c = pc).
aer2t(Ry,Rx,Rz) => T and t2aer(T) => Ry,Rx,Rz Euler angle conversions
vd2tfm(x,y,z,Rx,Ry,Rz) and tfm2vd(T) 4x4 matrix <=> VirDir tx ty tz rx ry rz
eq2dms(v3) => text "hh:mm.m +dd:mm:ss dist"
radec2eq(ra,dec,dist) (ra h:m:s, dec d:m:s, dist) => J2000 3-vector
radec2eqbasis(ra,dec) => 3x3matrix (ra,dec DEGREES -> XY=sky-plane, +Ynorth)
Tprecess(fromyr,toyr) => 3x3matrix: Pfrom * Tprecess(fromyear,toyear) = Pto
\@Tab = 3x3matrix; a,b are: g(alactic) s(upergalactic) e(qJ2000) z(ecliptic) c(C-gal)
list("string") converts blank/comma/brace-separated string to list
put( list ) print N-vector, or 3x3 or 4x4 matrix
pt( list ) print list on one line (for copy/pasting)
Each line is a perl "eval", e.g.: \@a = (1,2,3); print vdot(\@a,\@a); sub me {...}
put( list ) print N-vector, or 2x2 or 3x3 or 4x4 matrix
pt( list ) print list on one line, full precision (for copy/pasting)
Each inputline is a perl "eval", e.g.: \@a = (1,2,3); print dot(\@a,\@a); sub me {...}
Previous line's answer saved in "\@_"; first scalar saved in \"\$_\".
For multiline input, have a not-yet-closed "{", or use "\\" at end of line.
EOF
}
......@@ -352,7 +358,7 @@ sub dms2rad {
sub dms2d {
local($d,$m,$s) = @_;
if($m eq "" && $d =~ /:/) {
if($d =~ /:/) {
($d,$m,$s) = split(/:/, $d);
}
local($sign) = ($d =~ s/^\s*-//) ? -1 : 1;
......@@ -369,7 +375,9 @@ sub rad2dms {
local($sign) = $_[1] || " ";
$sign = "-", $d = -$d if $d<0;
$d *= 180/$pi;
return sprintf("%s%02d:%04.1f", $sign, int($d), 60*($d - int($d)));
local($m) = 60*($d - int($d));
local($s) = 60*($m - int($m));
return sprintf("%s%02d:%02d:%04.1f", $sign, int($d), int($m), $s);
}
# &eqd2vec(ra(decimaldegrees), dec(decimaldegrees), dist) => J2000 3-vector
......@@ -390,16 +398,20 @@ sub radec2eq {
&eqd2vec($ra, $dec, $r);
}
# Supergalactic (sgx, sgy, sgz) to J2000 equatorial (RA, Dec, Dist)
sub sg2eq {
&eq2dms( &vm3mul( @_, @Tse ) );
}
# Equatorial (x, y, z) to supergalactic (sgx, sgy, sgz)
sub eq2sg {
&vm3mul( @_, @Tes );
}
sub radec2eqbasis {
local($ra, $dec) = @_;
$ra = &hms2d($ra) if $ra =~ /:/;
$dec = &dms2d($dec) if $dec =~ /:/;
$ra *= $pi/180;
$dec *= $pi/180;
# vector to galaxy (equatorial coords)
......@@ -422,6 +434,18 @@ sub eq2dms {
return(&rad2dms($ra/15), &rad2dms($dec,"+"), sqrt(&dot(@eq,@eq)));
}
# &gal2dms( Galactic 3-vector ) => ( londd:mm:ss.s, latdd:mm:ss.s, dist )
sub gal2dms {
local($r, $rxy);
local(@eq) = @_;
local($lon, $lat);
$rxy = sqrt($eq[0]*$eq[0] + $eq[1]*$eq[1]);
$lat = atan2($eq[2], $rxy);
$lon = atan2($eq[1], $eq[0]);
$lon += 2*$pi if($lon < 0);
return(&rad2dms($lon), &rad2dms($lat,"+"), sqrt(&dot(@eq,@eq)));
}
sub init_eq2gal {
$pi = 3.14159265358979;
......@@ -457,11 +481,75 @@ sub init_eq2gal {
@Tcg = @Tgc = (-1,0,0, 0,-1,0, 0,0,1);
@Tcs = &m3mmul( @Tcg, @Tgs );
@Tsc = &m3mmul( @Tsg, @Tgc );
@Tce = &m3mmul( @Tcg, @Tge );
@Tec = &transpose( @Tce );
@Tze = &m3( &tfm('x', 23.4393) ); # Ecliptic ("zodiac") to J2000 (3x3 matrix)
@Tez = &transpose(@Tze);
@Tzc = &m3mmul( @Tze, @Tec );
@Tcz = &transpose(@Tzc);
@Tzg = &m3mmul( @Tze, @Teg );
@Tgz = transpose(@Tzg);
@Tzs = &m3mmul( @Tze, @Tes );
@Tsz = transpose(@Tzs);
}
# &Tprecess(fromyear, toyear)
# precession of equatorial coordinates from one epoch to another.
# returns 3x3 matrix: vm3mul( p_fromyear, &Tprecess(fromyear,toyear) ) = p_toyear
# Adapted from P.T. Wallace's SLALIB routine PREC.
# Applies to limited range of years. From the notes:
# 1) The epochs are TDB (loosely ET) Julian epochs.
#
# 2) The matrix is in the sense V(EP1) = RMATP * V(EP0)
#
# 3) Though the matrix method itself is rigorous, the precession
# angles are expressed through canonical polynomials which are
# valid only for a limited time span. There are also known
# errors in the IAU precession rate. The absolute accuracy
# of the present formulation is better than 0.1 arcsec from
# 1960AD to 2040AD, better than 1 arcsec from 1640AD to 2360AD,
# and remains below 3 arcsec for the whole of the period
# 500BC to 3000AD. The errors exceed 10 arcsec outside the
# range 1200BC to 3900AD, exceed 100 arcsec outside 4200BC to
# 5600AD and exceed 1000 arcsec outside 6800BC to 8200AD.
# The SLALIB routine sla_PRECL implements a more elaborate
# model which is suitable for problems spanning several
# thousand years.
# References:
# Lieske,J.H., 1979. Astron.Astrophys.,73,282. equations (6) & (7), p283.
# Kaplan,G.H., 1981. USNO circular no. 163, pA2.
sub Tprecess {
local($from, $to) = @_;
$to = 2000 if $to == 0 ;
$from = 2000 if $from == 0;
local($t0) = ($from - 2000) / 100; # "from" in centuries-from-2000
local($dt) = ($to - $from) / 100; # delta time in centuries
# delta-time, scaled so $dtdeg * arc-sec-per-century = degrees
local($dtdeg) = $dt / 3600;
local($w) = 2306.2181+(1.39656-0.000139*$t0)*$t0;
local($zetadeg) = ($w+((0.30188-0.000344*$t0)+0.017998*$dt)*$dt)*$dtdeg;
local($zdeg) = ($w+((1.09468+0.000066*$t0)+0.018203*$dt)*$dt)*$dtdeg;
local($thetadeg) = ((2004.3109+(-0.85330-0.000217*$t0)*$t0)
+((-0.42665-0.000217*$t0)-0.041833*$dt)*$dt)*$dtdeg;
local(@q) = &quatmul( &quatmul(
&ax2quat( 'z', $zetadeg ),
&ax2quat( 'y', -$thetadeg ) ),
&ax2quat( 'z', $zdeg ) );
&m3( &quat2t( @q ) );
}
sub t2quat {
local(@t) = (0+@_ == 16) ? @_ : &tfm(@_);
local(@t) = (0+@_ == 9) ? &m4(@_) : (0+@_ == 16) ? @_ : &tfm(@_);
local($s) = &mag( @t[0..2] );
......@@ -741,7 +829,7 @@ sub qrotbtwn {
# Direction is (x2,y2,z2) cross (x1,y1,z1)
local(@ijk) = &normalize(&cross(@_));
local($cost) = &dot(@_) / sqrt(&dot(@_[0..2,0..2]) * &dot(@_[3..5,3..5]));
local($sinhalf) = -sqrt((1 - $cost) / 2);
local($sinhalf) = sqrt((1 - $cost) / 2);
# magnitude of ijk is sin(angle/2)
return ( sqrt((1+$cost)/2), $ijk[0]*$sinhalf, $ijk[1]*$sinhalf,
$ijk[2]*$sinhalf );
......@@ -782,6 +870,7 @@ sub pt {
printf "%.11g%s", shift(@_), (@_>0?" ":"\n");
}
}
"";
}
......@@ -856,7 +945,7 @@ sub round {
sub basis {
local($d) = shift;
local(@done) = (0) x $d;
local(@M, @T, @v, $j);
local(@M, @T, @v, $i, $j, $row);
# @M is the list of vectors assigned so far.
# @T is the final output array, with members of @M
# arranged in appropriate rows.
......@@ -894,9 +983,27 @@ sub basis {
@T[$row*$d .. $row*$d+$d-1] = @v;
$done[$row] = 1;
}
local($det) = &det(@T);
if($det < 0) {
@T[$row*$d..$row*$d+$d-1] = &svmul( -1, @T[$row*$d..$row*$d+$d-1] );
}
return @T;
}
sub det {
if(@_ == 1) {
return $_[0];
} elsif(@_ == 4) {
return $_[0]*$_[3] - $_[1]*$_[2];
} elsif(@_ == 9) {
return &dot( &cross( @_[0..2, 3..5] ), @_[6..8] );
} else {
# Oh well, maybe later.
print STDERR "Oops, ignoring determinant of ", sqrt(0+@_), "-rank matrix.\n";
return 1;
}
}
# Transpose a square matrix.
sub transpose {
local($d) = int(sqrt(@_));
......@@ -1012,16 +1119,31 @@ sub tfm_interact {
# If we were invoked as a shell command, act as a perl calculator.
local($tty) = (-t STDIN);
print STDERR "Type \"help\" for help\n> " if $tty;
@prefix = ();
while(($lastinput = <>) ne "" && $lastinput !~ /^(q|quit|exit|bye)$/i) {
$lastinput =~ s/^[\s>]*//;
$lastinput = "&help " if $lastinput =~ /^\s*\?\s*$/;
chop $lastinput;
push(@HIST, $lastinput);
@_ = eval($lastinput);
if($lastinput !~ /;$/ && $lastinput ne "" && $lastinput !~ /^&*help$/) {
$lastinput = "&help ", @prefix = () if $lastinput =~ /^\s*\?\s*$/;
chomp $lastinput;
if($lastinput =~ /(.*)\\$/) {
push(@prefix, $1);
print STDERR "\t" if $tty;
next;
}
$wholeinput = join("\n\t", @prefix, $lastinput);
if(($wholeinput =~ tr/{/{/) > ($wholeinput =~ tr/}/}/)) {
push(@prefix, $lastinput);
print STDERR "\t" if $tty;
next;
}
push(@HIST, $wholeinput);
@_ = eval($wholeinput);
if($@ != "") {
print $@, "\n";
} elsif($wholeinput !~ /;$/ && $wholeinput ne "" && $wholeinput !~ /^&*help$/) {
&put(@_);
$_ = $_[0];
}
@prefix = ();
print STDERR "> " if $tty;
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Contact help@ncsa.illinois.edu with questions regarding this site. | All rights reserved. ©2024 Board of Trustees of the University of Illinois. | Web Privacy Notice