[Date Prev][Date Next][Thread Prev][Thread Next] [Search] [Date Index] [Thread Index]

[MacPerl] Still More Perl Starmapping



I wasn't expecting to have this ready so soon!

You'll need a copy of the "Preliminary Version of the Third Catalog of Nearby Stars"

Preliminary Version of the Third Catalogue of Nearby Stars
     GLIESE W., JAHREISS H.
    <Astron. Rechen-Institut, Heidelberg (1991)>

which is available from

http://adc.gsfc.nasa.gov/

it's catalog 5070A.  Save this script as a droplet, and drag the catalog onto it.

Now, I am going to be adding some more MacPerl-specific things to this eventually, but I feel it's more important to get the core code working. This is now nearly identical in function to the C code I had as a model (although I didn't really borrow any techniques from there (except for the math, which is simple high-school trig)).  It's missing the TERSE output mode.  It does however flag potential locations for habitable planets, and can exclude based on habitability (which the original didn't do).

Richard Grubb said he'd work on the "nearest neighbor" routines, which produce the "trade routes" between stars in the catalog.  I've been thinking about doing red/green stereograms as a possible output.  At nearly 500 lines, this is my longest program to date, and it will probably only get longer.

Comments welcome.

#!usr/bin/perl -w

# USAGE: catalog.plx File_Name
#          where File_Name is a star catalog.
# Mac USAGE: save as droplet.  Drag & Drop target input file
# or just double-click on the droplet and let StandardFile do
# it for you.
#
# ver 0.5a.r1 
#    Brian McNett <webmaster@mycoinfo.com>
#    Richard L. Grubb, <richard.l.grubb@boeing.com>
#
# to do: Convert the input to a hash
#        This may require a complete
#        rewrite.

use strict;
use Math::Trig qw(deg2rad);
use vars qw(
    $file $input_file $output_file

    $i $strin_len

    $ident $Comp $DistRel
    $RAh $RAm $RAs
    $DE_sign $DEd $DEm
    $pm $u_pm $pmPA
    $RV $n_RV $Sp $r_Sp
    $Vmag $r_Vmag $n_Vmag
    $B_V $r_B_V $n_B_V
    $U_B $r_U_B $n_U_B
    $R_I $r_R_I $n_R_I 
    $trplx $e_trplx
    $plx $e_plx $n_plx 
    $Mv $n_Mv $q_Mv
    $U $V  $W 
    $HD $DM $Giclas $LHS $Other
    $Remarks

    $alpha $delta
    $phi $theta $rho
    $rVect $Xe $Ye $Ze
    $Xg1950 $Yg1950 $Zg1950
    $Xg2000 $Yg2000 $Zg2000
    
    $luminosity $habitable

    $ref_dist

    %settings

);

# for Mac users who forget to drag&drop the input file
if ($^O eq 'MacOS') {
    if ($#ARGV < 0) {
        eval "use Mac::StandardFile";
        $file = StandardGetFile('', 'TEXT');
        if ($file->sfGood()) {
            push(@ARGV, $file->sfFile());
        } else {
          exit(1);
        }
    }
}

open (STDERR, '>Dev:Console:Messages');

init_settings();

usr_console();

run_catalog();

close (IN_FILE);
close (OUT_FILE);

### SUBROUTINES ###

sub usr_console {
my ($input, $done);

until (defined($done)) {

        usage();

        $input = <STDIN>;
        chomp($input);
        if ($input =~ /^1$/)  {
            print "Enter reference point X coordinate: ";
            $settings{ref_x} = <STDIN>;
            chomp($settings{ref_x});
                if ($settings{ref_x} =~ /[A-Za-z]+|\s+/) {
                print "X coordinate must be numeric.\n";
                $settings{ref_x} = 0;
                exit (1);
                }
            print "Enter reference point Y coordinate: ";
            $settings{ref_y} = <STDIN>;
            chomp($settings{ref_y});
                if ($settings{ref_y} =~ /[A-Za-z]+|\s+/) {
                print "Y coordinate must be numeric.\n";
                $settings{ref_y} = 0;
                exit (1);
                                }
            print "Enter reference point Z coordinate: ";
            $settings{ref_z} = <STDIN>;
            chomp($settings{ref_z});
            if ($settings{ref_z} =~ /[A-Za-z]+|\s+/) {
                print "Z coordinate must be numeric.\n";
                $settings{ref_z} = 0;
                exit (1);
            }
        }
        elsif ($input =~ /^2$/) {
            print "Minimum distance: ";
            $settings{mindist} = <STDIN>;
            chomp($settings{mindist});
            if ($settings{mindist} =~ /[A-Za-z]+|\s+/) {
                print "Distance needs to be numeric.\n";
                $settings{mindist} = 0;
                exit (1);
            }
        }
        elsif ($input =~ /^3$/) {
            print "Maximum distance: ";
            $settings{maxdist} = <STDIN>;
            chomp($settings{maxdist});
            if ($settings{maxdist} =~ /[A-Za-z]+|\s+/) {
                print "Distance must be numeric.\n";
                $settings{maxdist} = 5;
                exit (1);
            }
        }
        elsif ($input =~ /^4$/) {
            print "Minimum Luminosity: ";
            $settings{minlum} = <STDIN>;
            chomp($settings{minlum});
            if ($settings{minlum} =~ /[A-Za-z]+|\s+/) {
                print "Luminosity must be numeric.\n";
                $settings{minlum} = 0;
                exit (1);
            }
        }
        elsif ($input =~ /^5$/) {
            print "Maximum Luminosity: ";
            $settings{maxlum} = <STDIN>;
            chomp($settings{maxlum});
            if ($settings{maxlum} =~ /[A-Za-z]+|\s+/) {
                print "Luminosity must be numeric.\n";
                $settings{maxlum} = 200.0;
                exit (1);
            }
        }
        elsif ($input =~ /^6$/) {
            if ($settings{allow_binaries} eq "YES") {
                $settings{allow_binaries} = "NO";
            } else { $settings{allow_binaries} = "YES" }
        }
        elsif ($input =~ /^7$/) {
            if ($settings{just_habitables} eq "NO") {
                $settings{just_habitables} = "YES";
            } else { $settings{just_habitables} = "NO" }
        }
        elsif ($input =~ /^8$/) {
            if ($settings{output_style} eq "VERBOSE") {
                $settings{output_style} = "TERSE";
            } else { $settings{output_style} = "VERBOSE" }
        }
        elsif ($input =~ /^9$/) {
            print "Enter new filename: ";
            $settings{output_filename} = <STDIN>;
            chomp($settings{output_filename});
        }
        elsif ($input =~ /^0$/) {
            $done = "1";
        }
        else
        {
        print "I haven't the foggiest notion what you mean by $input\n";
        }
    }
}

sub run_catalog {
    print "Opening input data...\n";

    $input_file = $ARGV[0];
    open( IN_FILE, $input_file) || die "Can't open file $input_file $!\n" ;

    print "parsing\n";

    $output_file = $settings{output_filename};
    open (OUT_FILE, ">$output_file") || die "Can't open output $output_file $!\n" ;

    while (<IN_FILE>) {

        $i++;

        parse_gliese();

        calc_cartesian();    # Convert Location to Cartesian System

        calc_ref_distance(); # calculate distance to user-supplied reference point

# Calculate Neighbors

# Determine Habitable Systems
# ---------------------------
        is_habitable();

# Insert Planets, Moons

# Output logfile and write out resource compiler sources

        output_log();

        print  "[$i] $ident $Remarks\n"

    }

print "DONE!  Yeeehah!";
}


sub calc_cartesian {
if ($RAh eq '  ') {

    $rho    = 0.000;

    $Xe     = 0.000;
    $Ye     = 0.000;
    $Ze     = 0.000;

    $Xg1950 = 0.000;
    $Yg1950 = 0.000;
    $Zg1950 = 0.000;

    $alpha = 0;
    $delta = 0;

#   $Xg2000 = 0.000;
#   $Yg2000 = 0.000;
#   $Zg2000 = 0.000;

} else {

    $rho = 1000/$plx;

    $alpha = (($RAh * 15)+($RAm * 0.25)+($RAs * 0.0041666));
    $delta = join('', $DE_sign, ($DEd + $DEm/60));

    # calculate Phi, Theta (in radians)
    # ---------------------------------
    $phi     = deg2rad($alpha);
    $theta   = deg2rad($delta);

    # Earth-equatorial cartesian coords
    # ---------------------------------
    $rVect = $rho   * cos($theta);
    $Xe    = $rVect * cos($phi);
    $Ye    = $rVect * sin($phi);
    $Ze    = $rho   * sin($theta);

    # Galactic-equatorial cartesian coords (B1950)
    # --------------------------------------------
    $Xg1950 = -(0.0672 * $Xe) - (0.8727 * $Ye) - (0.4835 * $Ze);
    $Yg1950 =  (0.4927 * $Xe) - (0.4504 * $Ye) + (0.7445 * $Ze);
    $Zg1950 = -(0.8676 * $Xe) - (0.1884 * $Ye) + (0.4602 * $Ze);

    # Galactic-equatorial cartesian coords (J2000)
    # commented out until implementation
    # --------------------------------------------
    # $Xg2000 = -(0.0550 * $Xe) - (0.8734 * $Ye) - (0.4839 * $Ze);
    # $Yg2000 =  (0.4940 * $Xe) - (0.4449 * $Ye) + (0.7470 * $Ze);
    # $Zg2000 = -(0.8677 * $Xe) - (0.1979 * $Ye) + (0.4560 * $Ze);
    }
}

sub parse_gliese {
    chomp;

    $strin_len = length($_);       # check how long the input string is

    $ident    = substr($_,  0, 8); #A8   #Identifier ; see remarks.
    $Comp     = substr($_,  8, 2); #A2   #Components (A,B,C,... )
    $DistRel  = substr($_, 10, 1); #A1   #[pqsx] Reliability of the distance
                                   #12 is unused
    $RAh      = substr($_, 12, 2); #I2   #? Right Ascension B1950 (hours)
                                   #15 is unused
    $RAm      = substr($_, 15, 2); #I2   #? Right Ascension B1950 (minutes)
                                   #18 is unused
    $RAs      = substr($_, 18, 2); #I2   #? Right Ascension B1950 (seconds)
                                   #21 is unused
    $DE_sign  = substr($_, 21, 1); #A1   #Declination B1950 (sign)
    $DEd      = substr($_, 22, 2); #I2   #? Declination B1950 (degrees)
                                   #25 is unused
    $DEm      = substr($_, 25, 4); #F4.1 #? Declination B1950 (minutes)
                                   #30 is unused
    $pm       = substr($_, 30, 6); #F6.3 #? Total proper motion
    $u_pm     = substr($_, 36, 1); #A1   #Uncertainty flag (:) on pm
    $pmPA     = substr($_, 37, 5); #F5.1 #? Direction angle of proper motion
                                   #43 is unused
    $RV       = substr($_, 43, 6); #F6.1 #? Radial velocity
                                   #50 is unused
    $n_RV     = substr($_, 50, 3); #A3   #Remark on RV
                                   #54 is unused
    $Sp       = substr($_, 54,12); #A12  #Spectral type or color class
    $r_Sp     = substr($_, 66, 1); #A1   #Selected sources
    $Vmag     = substr($_, 67, 6); #F6.2 #Apparent magnitude
    $r_Vmag   = substr($_, 73, 1); #A1   #Note on origin of magnitude
    $n_Vmag   = substr($_, 74, 1); #A1   #[J] joint magnitude
    $B_V      = substr($_, 75, 5); #F5.2 #? color
    $r_B_V    = substr($_, 80, 1); #A1   #Note on origin of magnitude
    $n_B_V    = substr($_, 81, 1); #A1   #Joint color
    $U_B      = substr($_, 82, 5); #F5.2 #? color
    $r_U_B    = substr($_, 87, 1); #A1   #Note on origin of magnitude
    $n_U_B    = substr($_, 88, 1); #A1   #Joint color
    $R_I      = substr($_, 89, 5); #F5.2 #? color
    $r_R_I    = substr($_, 94, 1); #A1   #Note on origin of magnitude
    $n_R_I    = substr($_, 95, 1); #A1   #Joint color
    $trplx    = substr($_, 96, 6); #F6.1 #? Trigonometric parallax
    $e_trplx  = substr($_,102, 5); #F5.1 #? Standard error of trig. parallax
                                   #108 is unused
    $plx      = substr($_,108, 6); #F6.1 #? Resulting parallax
    $e_plx    = substr($_,114, 5); #F5.1 #? Standard error of res.  parallax
    $n_plx    = substr($_,119, 1); #A1   #[rwsop] Code on plx
                                   #121 is unused  
    $Mv       = substr($_,121, 5); #F5.2 #Absolute visual magnitude
    $n_Mv     = substr($_,126, 2); #A2   #Note on Mv, copied from cols 74-75
    $q_Mv     = substr($_,128, 1); #A1   #[a-f] Quality of absolute magnitude

    # Warning! CHECK THE VALIDITY OF THESE FIELDS BEFORE USING THEM!

    $U        = substr($_,131, 4); #I4   #? U space velocity component in the galactic plane and directed to the galactic center
                                   #136 is unused
    $V        = substr($_,136, 4); #I4   #? V space velocity component in the galactic plane and in the direction of galactic rotation
                                   #141 is unused
    $W        = substr($_,141, 4); #I4   #? W space velocity component in the galactic plane and in the direction of the North Galactic Pole

    # We can't be certain the line doesn't just end at this point, so...

    if ($strin_len <= 145) {
    $HD = " ";
    } else {
    $HD       = (substr($_,146, 6)); #I6   #[15/352860]? designation
    }
    if ($strin_len <= 152) {
    $DM = " ";
    } else {
    $DM       = substr($_,153,12); #A12  #Durchmusterung number BD / CD / CP
    }
    if ($strin_len <= 165) {
    $Giclas = " ";
    } else {
    $Giclas   = substr($_,166, 9); #A9   #number        <...>
    }
    if ($strin_len <= 175) {
    $LHS = " ";
    } else {
    $LHS      = substr($_,176, 5); #A5   #number        <...>
    }
    if ($strin_len <= 181) {
    $Other = " ";
    } else {
    $Other    = substr($_,182, 5); #A5   #Other designations
    }
    if ($strin_len <= 187) {
    $Remarks = " ";
    } else {
    $Remarks  = substr($_,188,69); #A69  #Additional identifications (LTT, LFT, Wolf, Ross, etc.) and remarks
    }
}

sub is_habitable {
    $luminosity = 1/(2.5**($Mv - 4.85));

    if (($luminosity >= 0.24)
         && ($luminosity <= 2)
         && ($Sp !~ /[OBAMm]|I/)
       ) {
        $habitable = "Habitable";
        } else {
        $habitable = "Uninhabitable";
    }
}

sub output_log {

unless ($settings{just_habitables} =~ /^NO$/) {
    if ($habitable eq "Uninhabitable") { return(1) }
    }
unless ($settings{allow_binaries} =~ /^YES$/) {
    if ($Comp !~ /\s/) { return(1) }
}

if (($ref_dist     >= $settings{mindist})
   && ($ref_dist   <= $settings{maxdist})
   && ($luminosity >= $settings{minlum} )
   && ($luminosity <= $settings{maxlum})) {
        
printf OUT_FILE (<<'RECORD'
===========================================================
%s %s  Spectral Class: %s %s Luminosity: %4.5f (%s)
Distance  (Sol)  : %2.3f parsecs
Dist. (Reference): %2.3f parsecs
-----------------------------------------------------------
Remarks: %s
-----------------------------------------------------------
Coordinates: (Earth)    [ %2.3f , %2.3f , %2.3f ]
             (Galactic) [ %2.3f , %2.3f , %2.3f ]
-----------------------------------------------------------
Alpha: %2.5f °       Delta: %2.5f °
RA: %s h %s m %s s   Dec:  %s%s ° %s '
HD#: %s              Giclas: %s
LHS: %s              DM#: %s
Other: %s
RECORD
  ,
  $ident, $Comp, $Sp, $r_Sp ,$luminosity, $habitable,
  $rho, $ref_dist,
  $Remarks,
  $Xe, $Ye, $Ze,
  $Xg1950, $Yg1950, $Zg1950,
  $alpha, $delta, $RAh, $RAm, $RAs, $DE_sign,
  $DEd, $DEm,
  $HD, $Giclas, $LHS, $DM, $Other
  );
  }
}

sub usage {
    print "Enter an item number as indicated below: \n\n";
    printf STDOUT (<<'    USAGE'
    Current output settings are
    ---------------------------
    1). Reference Point      : %2.3f  %2.3f  %2.3f
    2). Minimum Distance     : %2.3f
    3). Maximum Distance     : %2.3f
    4). Minimum Luminosity   : %3.5f
    5). Maximum Luminosity   : %3.1f
    6). Binaries Allowed?    : %s
    7). Only Habitable Stars?: %s
    8). Output Style         : %s
    9). Output File Name     : %s
    ----------------------------
    Enter a number to modify the output settings,
    or type "0" to run.
    USAGE
    ,
    $settings{ref_x}, $settings{ref_y}, $settings{ref_z},
    $settings{mindist},
    $settings{maxdist},
    $settings{minlum},
    $settings{maxlum},
    $settings{allow_binaries},
    $settings{just_habitables},
    $settings{output_style},
    $settings{output_filename});
}

sub init_settings {
    %settings = (
    'ref_x'           => 0.0,
    'ref_y'           => 0.0,
    'ref_z'           => 0.0,
    'mindist'         => 0.0,
    'maxdist'         => 5.0,
    'minlum'          => 0.0,
    'maxlum'          => 200.0,
    'allow_binaries'   => 'YES',
    'output_style'     => 'VERBOSE',
    'output_filename'  => "gliese3.log",
    'interactive'     => 'ON',
    'just_habitables' => 'NO'
    );
}

sub calc_ref_distance {

if (($settings{ref_x} == 0) &&
    ($settings{ref_y} == 0) &&
    ($settings{ref_z} == 0)) {
    $ref_dist = $rho
    } else {
    $ref_dist = sqrt(
    (($settings{ref_x} - $Xg1950)**2) +
    (($settings{ref_y} - $Yg1950)**2) +
    (($settings{ref_z} - $Zg1950)**2)
    );
    }
  }

__END__


# Fungal Parataxonomy                   Mycology Information (Mycoinfo)
# Webmaster, Staff Writer      **The World's First Mycology E-Journal**   
# <mailto:webmaster@mycoinfo.com>            <http://www.mycoinfo.com/> 
#
# First they ignore you. Then they laugh at you. Then they fight you.
# Then you win.                                     --Mohandas Gandhi


===== Want to unsubscribe from this list?
===== Send mail with body "unsubscribe" to macperl-request@macperl.org