package CelCoords; # Author: Masaaki Sakano - SSC/LUX our $VERSION = '4.21'; use strict; use Exporter 'import'; # our @EXPORT = qw(CelCoords angulardistance to_s); # symbols to export in default # our @EXPORT_OK = qw(); # symbols to export on request use strict; use Coords; use Math::Trig; # pi, acos, rad2deg, deg2rad... use POSIX qw(fmod); # fmod: Float MODulo. use Data::Dumper; use SSCLib qw(isnumber hashmerge sscnote sscwarn); our $Taskname = 'CelCoords'; =head1 NAME CelCoords - Class for Celestial-Coordinates objects. =head1 SYNOPSIS use CelCoords; =head1 DESCRIPTION Class for the celestial coordinates. This is a subclass of Coords class. Many methods are defined in Coords class. The method to_s is overwritten. =item new ( ARGUMENTS ) Input is either a hash in the form of ("x%val"=>123.45, "y%val"=>-67.8, "x%unit"=>"deg", "y%unit"=>"deg", "all%sys"=>"FK5", "all%epoch"=>2000.0) or an array with 4 components (x[deg], y[deg], COORDSYS, EPOCH), such as, (123.45, -67.8, "FK5", 2000.0), in which case all must be there in this order. The hash key parameters can include "allarg" and "arrayhash". See the reference of Coords.pm for detail --- the interface is the same, except the aliases of "sys" and "epoch" (for "all%sys" and "all%epoch", respectively) are defined for the input in this class (nb., they nevertheless are stored with the key "all%..." in the object as usual). Note all%type=="spherical-celestial" is automatically assigned in default unless otherwise explicitly given by the user. =item angulardistance ( ARGUMENT ) If this is called as Class method, the arguments are either (CelCoords1, CelCoords2) or (RA1, DEC1, RA2, DEC2) [degree]. my ($angle) = &CelCoords::angulardistance($c1, $c2); my ($angle) = &CelCoords::angulardistance($ra1, $dec1, $ra2, $dec2); If this is called as normal method, the argument is (CelCoords2). my ($angle) = $c1->angulardistance($c2); Returns the angular distance between the two CelCoords objects in unit the same as that of the X-axis of the first object ($c1->{axis}->[0]->{x}). So far, only the valid units are degree and radian. If something is wrong, returns a negative value (-1). =item distance ( ARGUMENT ) Alias for angulardistance(). =item to_s ( ) Returns the string of the instance (to print). =cut ; # For Emacs ############################################################ our @ISA = qw(Coords); # Inheritance #--------------------------------------------------------------------- sub new { #--------------------------------------------------------------------- ### Input: ("x%val"=>123.45, "y%val"=>-67.8, "x%unit"=>"deg", "y%unit"=>"deg", "sys"=>"FK5", "epoch"=>2000.0) ### or (123.45, -67.8, "FK5", 2000.0) (All must be there in this order). ### If the second form is chosen, the units have to be in degree. ### Return: (Instance of CelCoords class) ### Description: ### Constraints: ### Global: ### Usage example: my $class = shift; # The first argument is CelCoords my @ararg = @_; # my $self = {}; my($self); my(@ar, $s1, $s2); my(%h) = ('all%type'=>"spherical-celestial"); ## Reading arguments if ( (4 == scalar(@ararg)) && isnumber($ararg[0]) # isnumber in SSCLib && isnumber($ararg[1]) && isnumber($ararg[3]) ) { # Array input %h = ('x%val' => $ararg[0], 'y%val' => $ararg[1], 'all%sys' => $ararg[2], 'all%epoch' => $ararg[3], 'x%unit' => 'deg', 'y%unit' => 'deg', ); } else { # Hash input my %hsarg = @_; while (scalar(@ar=each(%hsarg)) > 0){ if ($ar[0] =~ /^(allarg|arrayhash|(sub)?name)$/i) { $h{lc($&)} = $ar[1]; } elsif ($ar[0] =~ /^(sys|epoch)$/) { ## If another new top-level keyword is invented, then ## the above "IF" condition should be modified as such. $h{'all%'.$&} = $ar[1]; # Special form for (sys|epoch) } elsif ($ar[0] =~ /^([\w_\-]+)%([\w_\-]+)$/) { $h{$ar[0]} = $ar[1]; } else { $h{'all%'.$ar[0]} = $ar[1]; &sscwarn("The arguments in new(@ararg) does not look right (Pair==($ar[0])($ar[1])).") if isnumber($ar[0]); # isnumber in SSCLib } } } $self = Coords->new(%h); # $self = {}; bless $self, $class; # Inheritance is applied (with this 2-arguments bless), if there is any. # $self->_initialize(@ararg); return $self; } #--------------------------------------------------------------------- sub distance { #--------------------------------------------------------------------- ### Alias for angulardistance(), ### to be used in isNearlyEqual() in the base class, Coords. return(angulardistance(@_)); } #--------------------------------------------------------------------- sub angulardistance { #--------------------------------------------------------------------- ### Input: (CelCoords1, CelCoords2) or ### : (RA1, DEC1, RA2, DEC2) [degree] if called as Class method. ### : (CelCoords2) if called as normal method. ### Return: Angular distance between two CelCoords objects ### in unit the same as that of the X-axis of the first object ($c1->{axis}->[0]->{x}). ### So far, only the valid units are degree and radian. ### If something is wrong, returns a negative value (-1). ### Description: my $self=shift; my $halfPI = pi/2; # pi defined Math::Trig my $RETWRONG = -1; my $cosinegamma; my $isdeg = 0; # False if (ref($self)) { # Normal method. return( ref($self)->angulardistance($self,@_) ); # Practically calling the following. } else { # Class method. my($alpha1,$delta1, $alpha2,$delta2); # in radian. my $retangle; if ((scalar(@_) == 4) && (isnumber(@_)) ) { # isnumber in SSCLib $isdeg = 1; # True in default. $alpha1 = °2rad($_[0]); # °2rad() in Math::Trig $delta1 = °2rad($_[1]); $alpha2 = °2rad($_[2]); # -2pi < alpha < 2pi, guaranteed. $delta2 = °2rad($_[3]); $alpha1 += 2*pi if $alpha1 < 0; # pi in Math::Trig $alpha2 += 2*pi if $alpha2 < 0; } else { my $c1=shift; my $c2=shift; $isdeg = 1 if ($c1->{axis}->[0]->{unit} && ($c1->{axis}->[0]->{unit} =~ /^(arc)?deg(ree)?/i)); if ( (($c1->{all}->{sys} && $c2->{all}->{sys}) && ($c1->{all}->{sys} ne $c2->{all}->{sys})) || ((defined($c1->{all}->{epoch}) && defined($c2->{all}->{epoch})) && ($c1->{all}->{epoch} ne $c2->{all}->{epoch}) ) ) { return($RETWRONG); } if ( defined($c1->{axis}->[0]) && defined($c1->{axis}->[1]) && defined($c2->{axis}->[0]) && defined($c2->{axis}->[1]) ) { $alpha1 = CelCoords->get_longitude_radian_positive($c1->{axis}->[0]); $delta1 = CelCoords->get_latitude_radian($c1->{axis}->[1]); $alpha2 = CelCoords->get_longitude_radian_positive($c2->{axis}->[0]); $delta2 = CelCoords->get_latitude_radian($c2->{axis}->[1]); } else { return undef; } } $cosinegamma = cos($halfPI-$delta1) * cos($halfPI-$delta2) + sin($halfPI-$delta1) * sin($halfPI-$delta2) * cos($alpha1-$alpha2); # http://www.krysstal.com/sphertrig.html # http://spiff.rit.edu/classes/phys440/lectures/coords/coords.html $retangle = &acos($cosinegamma); # in Math::Trig # $retangle = atan2(sqrt(1-$cosinegamma**2), abs($cosinegamma)); if ($isdeg) { return( &rad2deg($retangle) ); } else { return( $retangle ); } } } #--------------------------------------------------------------------- sub get_longitude_radian_positive { #--------------------------------------------------------------------- ### Input: (longitude, [unit(Def:deg)]) ### or ($coord->{axis}->[0]) # Must have ->{val, unit}. ### Return: Modulo of the longitude (0<=RA<360[deg]) in radian. ### Description: Class method, but can be used as normal method. my $self=shift; # Ignore my $alpha = shift; # Longitude my $unit = "deg"; # Default $unit = shift if (scalar(@_) > 0); if ((ref $alpha) and (defined $alpha->{val})) { if ((defined $alpha->{unit}) and ($alpha->{unit} ne "")){ $unit = $alpha->{unit}; } $alpha = $alpha->{val}; } if ( $unit =~ /^(arc)?deg(ree)?/i) { $alpha = °2rad($alpha); # °2rad() in Math::Trig } elsif ($unit =~ /^rad(ian)?/i) { $alpha = &fmod($alpha, 360); # &fmod() in POSIX } else { ## Do nothing (namely, passed as it is). } ## Convert the potential negative into positive. $alpha += 2*pi if ($alpha < 0); # pi in Math::Trig return($alpha); } #--------------------------------------------------------------------- sub get_latitude_radian { #--------------------------------------------------------------------- ### Input: (latitude, [unit(Def:deg)]) ### or ($coord->{axis}->[1]}) # Must have ->{val, unit}. ### Return: The latitude in radian. ### Constraint: The input has to be in the range of -90<=delta<=90 [deg]. ### Description: Class method, but can be used as normal method. my $self=shift; # Ignore my $delta = shift; # Latitude my $unit = "deg"; # Default $unit = shift if (scalar(@_) > 0); if ((ref $delta) and (defined $delta->{val})) { if ((defined $delta->{unit}) and ($delta->{unit} ne "")){ $unit = $delta->{unit}; } $delta = $delta->{val}; } if ( $unit =~ /^(arc)?deg(ree)?/i) { $delta = °2rad($delta); # °2rad() in Math::Trig } elsif ($unit =~ /^rad(ian)?/i) { # Do nothing. } return($delta); } #--------------------------------------------------------------------- sub to_s { #--------------------------------------------------------------------- ### Input: (['indent'=>"", 'withname'=>1(Default)]) ### Return: String expression of the instance (to print). ### Description: ### If the optional indent, such as, "| " is given, it is added ### to the head of each line. ### Overwriting "to_s" method in the base class, Coords. my $self = shift; my %arghashpartdef = ('indent' => "", 'withname' => 1); my %hsarg = @_; %hsarg = hashmerge(\%hsarg, \%arghashpartdef); # in SSCLib my $indent = ($hsarg{'indent'} || ""); my(@arunit, $i, $s, $strundef, $fmt); $strundef = sprintf("%11s", 'UNDEFINED'); my(%hsnum) = ('val' => [$strundef, $strundef], 'CDELT'=> [$strundef, $strundef], 'CROTA'=> [$strundef, $strundef], ); my(%hssysepoch) = ('sys'=>'UNDEFINED', 'epoch'=>'UNDEFINED'); my(%hssysepochfmt) = ('epoch'=>'%f'); my(%flaghsnum) = (); foreach $s (keys(%hsnum)) { $flaghsnum{$s} = 0; # Flag whether the key is defined or not. } foreach $i (0..1) { $arunit[$i] = 'UNDEFINED'; if ((ref($self->{axis}) eq 'ARRAY') && ($#{$self->{axis}} >= $i)) { if ( ref($self->{axis}->[$i]) eq 'HASH' ) { foreach $s (keys(%hsnum)) { if ( defined($self->{axis}->[$i]{$s}) ) { #(exists( $self->{axis}->[$i]{val})) && # Unnecessary $flaghsnum{$s} = 1; # The keyword is defined in at least 1 axis. if ($s eq 'val') { $hsnum{$s}[$i] = sprintf("%11.6f", $self->{axis}->[$i]->{$s}); } else { $hsnum{$s}[$i] = sprintf("%e", $self->{axis}->[$i]->{$s}); } } } if ( defined($self->{axis}->[$i]{unit}) ) { $arunit[$i] = $self->{axis}->[$i]->{unit}; } } } } # foreach $i (0..1) foreach $s (keys(%hssysepoch)) { if (defined($self->{all}->{$s})) { $fmt = ($hssysepochfmt{$s} || '%s'); $hssysepoch{$s} = sprintf($fmt, $self->{all}->{$s}); } } my ($retstr) = ""; ### 1st line $retstr = sprintf("NAME (SUB): %s (%s)\n" , ($self->{name} ? $self->{name} : "NONE") , ($self->{subname} ? $self->{subname} : "NONE") ) if ($hsarg{'withname'}); ### 2nd line #return sprintf("(%11.6f, %11.6f) [%s, %s: (%s, epoch=%s)]" $retstr .= sprintf("(%s, %s) [%s, %s] (%s, epoch=%s)" , $hsnum{'val'}[0], $hsnum{'val'}[1], $arunit[0], $arunit[1] # , $arval[0], $arval[1], $arunit[0], $arunit[1] , $hssysepoch{'sys'}, $hssysepoch{'epoch'} ); foreach $s (qw(CDELT CROTA)) { if ($flaghsnum{$s}) { $retstr .= sprintf(" %s=(%s, %s)", $s, $hsnum{$s}[0], $hsnum{$s}[1]); } } $retstr .= "\n"; ### 3rd line ## Common things foreach $i (keys(%{$self->{all}})) { next if ($i =~ /^(sys|epoch)$/); if (defined $self->{all}->{$i}) { $retstr .= " $i=(" . $self->{all}->{$i} . "),"; } } $retstr = substr($retstr,0,length($retstr)-1); # Remove "," in the end $retstr .= "\n"; ### 4th and following lines foreach $i (keys(%{$self->{comment}})) { if (defined $self->{comment}->{$i}) { $retstr .= sprintf("# %s: %s\n", $i, $self->{comment}->{$i}); } } $retstr =~ s/^/$indent/gm; # Add an indent, if needed. return($retstr); } sub DESTROY { # Called upon the garbage collection of objects. } 1;