deepsky.pl

来自「celestia源代码」· PL 代码 · 共 796 行 · 第 1/2 页

#!/usr/bin/perl -w
#
# Author: [e-mail]
#
#use Math::Libm ':all';
use Math::Trig;
use Math::Quaternion qw(slerp);
#
open(DSC, ">deepsky.dsc") || die "Can not create deepsky.dsc\n";
#
# Boiler plate
#-------------
#
($ver = "Revision: 1.20 ") =~ s/\$//g;
($me = $0) =~ s/.*\///;
($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = gmtime;
$year += 1900;
$mon += 1;
print DSC "# Revised NGC and IC Catalog, Wolfgang Steinicke, April 5, 2005\n";
print DSC "# http://www.ngcic.org/steinicke/default.htm\n";
print DSC "#\n";
print DSC "# Augmented by\n";
print DSC "# Revised 3rd Reference Catalog of Bright Galaxies (RC3,VII/155)\n";
print DSC "# http://cdsweb.u-strasbg.fr/viz-bin/Cat?VII/155\n";
print DSC "#\n";
print DSC "# Augmented by\n";
print DSC "# Mark III Catalog of Galaxy Peculiar Velocities (Willick+ 1997,VII/198)\n";
print DSC "# http://cdsweb.u-strasbg.fr/viz-bin/Cat?VII/198#sRM2.18\n";
print DSC "#\n";
print DSC "# Augmented by distances from\n";
print DSC "# The SBF Survey of Galaxy Distances. IV.\n";
print DSC "# SBF Magnitudes, Colors, and Distances, \n";
print DSC "# J.L. Tonry et al., Astrophys J 546, 681 (2001)\n";
print DSC "#\n";
print DSC "# Augmented by distances from\n";
print DSC "# Compilation of \"200 Brightest Galaxies\"\n";
print DSC "# http://www.anzwers.org/free/universe/galax200.html\n";
print DSC "#\n";
print DSC "# Using today's Hubble constant = 72 [km/sec/Mpc]  (WMAP 2005)\n";
print DSC "#\n";
print DSC "# Abreviations for various distance methods used:\n";
print DSC "# S= SBF (Surface Brightness Fluctuations), T-F= Tully-Fischer\n";
print DSC "# V = rad. velocity in CMB frame & Hubble law \n";
print DSC "# C=Cepheids, P=photometric, N(G)=planetary nebula (globular cluster) luminosity function\n";
print DSC "#\n";
print DSC "# Adapted for Celestia with Perl script: $me $ver\n";
print DSC "# Processed $year-$mon-$mday $wday $yday $isdst $hour:$min:$sec UTC\n";
print DSC "#\n";
print DSC "# by Dr. Fridger Schrempp, fridger.schrempp\@desy.de\n";
print DSC "# ------------------------------------------------------ \n";
print DSC "\n\n";
#
# Constants
#----------
#
$H0 = 72.0;      # today's Hubble constant [km/sec/Mpc], WMAP 2005
$pc2ly = 3.26167;       # parsec-to-lightyear conversion
$epsilon =23.4392911;       # ecliptic angle
$M_abs_ell = -20.133; # absolute (Blue) magnitude for elliptical galaxies
# from peak of SBF-survey distribution <=== Maple
$M_abs_spir = -19.07; # absolute (Blue) magnitude for spiral galaxies
# from peak of SBF-survey distribution <=== Maple
#
# Constant hashes and arrays
# --------------------------
#
# hash with distance methods used, key = $dist_flag
%method =  ("0" => "  # distance uncertain!\n",
            "1" => "  # method: S\n",
            "2" => "  # method: T-F\n",
            "3" => "  # method: V\n",
            "4" => "  # method: B200\n");

# Hubble types of nearby galaxies
%HType = ("LMC"      => "Irr",
          "SMC"      => "Irr",
          "Sculptor" => "E0",
          "Fornax"   => "E4",
          "WLM"      => "Irr",
          "Leo I"    => "E3");

# array with nearby galaxies
@local = ("LMC","SMC","Sculptor","Fornax","WLM","Leo I");

#
# read in revised SBF-survey data IV for fast lookup
#
open(SBF," < catalogs/sbf.txt")|| die "Can not read sbf.txt\n";
$n_sbf=0;
while (<SBF>) {
    next if (/^Pref_name/);
    next if (/^E\d+/ || /^U\d+/);
    $name = &ss(1,6);
    $name =~ s/^N0*(\w+)\s*/NGC\ $1/g;
    $name =~ s/^I0*(\w+)\s*/IC\ $1/g;
    #
    # SBF-distances [ly]
    # ------------------
    #
    $deltam   = &ss(77,81);
    $dm{$name} = $deltam;
    $dm{$name} =~ s/ //g;
    $sbf_distance{$name} = 10**(($deltam+5)/5)*$pc2ly;
    #$v_CMB_SBF{$name} = &ss(29,32);
    $n_sbf++;
}
close (SBF);
print "\n";
printf "%5d SBF galaxy distance data\n",$n_sbf;
#
# read in TF + Dn-sigma distance data for fast lookup
#
# read in first GINnumber ($gin), name and BT and B-V
open(TFE," < catalogs/TF/egalf1.dat") || die "Can not read egalf1.dat\n";
while (<TFE>) {
    $gin = &ss(1,5);
    $gin =~ s/ //g;
    $name = &ss(7,16);
    # get rid of extra spaces
    $name =~ s/^\s*//;
    $name =~ s/\s*$//;
    $name =~ s/\s+/ /g;
    $name =~ s/^N\s(\d+)\s*/NGC\ $1/g;
    $name =~ s/^I\s(\d+)\s*/IC\ $1/g;
    $name =~ s/^U\s(\d+)\s*/UGC\ $1/g;
    $name =~ s/^E0*(\d+)-G0*(\d+)\s*/ESO\ $1-$2/g;
    $name =~ s/^[ZC]0*(\d{1,3})_?0*(\w+)\s*/ZWG\ $1\.$2/g;
    $name =~ s/^A0*(\d+)-S*(\d+)/A\ $1-$2/g;
    $name =~ s/^DC(\d+)/DC\ $1/g;
    $name_f1{$gin} = $name;
}
close (TFE);
# next enter the corresponding TF distances
open(TF," < catalogs/TF/egalf2.dat") || die "Can not read egalf2.dat\n";
$n_tf_e=0;
while (<TF>) {
    $gin = &ss(1,4);
    $gin =~ s/ //g;
    $name = $name_f1{$gin};
    $distance_TFE = &ss(35,41);
    $distance_TFE =~ s/ //g;
    $distance_TFE{$name} = $distance_TFE/$H0 * 1.0e6*$pc2ly; # [ly]
    # $v_CMB_TFE{$name} = &ss(32,38) + &ss(40,46);
    #printf "%5d %15s   %10.4g\n",$gin,$name, $distance_TFE{$name};
    $n_tf_e++;
}
close(TF);
printf "%5d TF elliptical galaxy distance data\n",$n_tf_e;

open(TF," < catalogs/TF/tf-spiral.dat") || die "Can not read tf-spiral.dat\n";
$n_tf_s=0;
while (<TF>) {
    $pgc = &ss(6,10);
    $pgc =~ s/ //g;
    $name = &ss(12,20);
#printf "%15s",$name;
   $name =~ s/^NG?C?0*(\w+)\[?V?\]?\s*/NGC\ $1/g;
    $name =~ s/^IC?0*(\w+)\[?V?\]?\s*/IC\ $1/g;
    $name =~ s/^UG?C?0*(\w+)\s*/UGC\ $1/g;
    $name =~ s/^E?0*(\d+)-G?0*(\d+)\s*/ESO\ $1-$2/g;
    $name =~ s/^[ZC]0*(\d{1,3})_?0?(\w+)\s*/ZWG\ $1\.$2/g;
    $name =~ s/^A0*(\d+)-S*(\d+)/A\ $1-$2/g;
    $name_TFS{$pgc} = $name;
    $distance_TF = &ss(82,86);
    $distance_TF =~ s/ //g;
    $distance_TFS{$name} = $distance_TF/$H0 * 1.0e6*$pc2ly; # [ly]
    $logr = &ss(66,70);
    $logr =~ s/ //g;
    #printf "PGC %6d %15s   %10.4g\n",$pgc, $name, $distance_TFS{$name};
    $n_tf_s++;
}
close(TF);
printf "%5d TF spiral galaxy distance data\n",$n_tf_s;
#
# read in data on 200 brightest galaxies
#
open( B200," < catalogs/200.txt") || die "Can not read 200.txt\n";
while (<B200>) {
    next if (/^\s+/);
    $name = &ss(1,8);
    $name = &cmpname ($name) if ($name =~ /^N|^I/);;
    $name =~ s/\s*$//g;
    $messier = &ss(10,13);
    $messier =~ s/ //g;
    if ($messier){
    $messier =~ s/M(\d+)/M $1/;
    $name = $messier;
    }
    $distance_B200{$name} = &ss(74,79);
    $distance_B200{$name} =~ s/ //g;
    $distance_B200{$name} = 1e6*$distance_B200{$name};
    $method_B200{$name} = &ss(83,90);
    $method_B200{$name} =~ s/\s+$//g;
    #print "$name   $distance_B200{$name} $method_B200{$name}\n";
}
close(B200);
#
# read in revised RC3 data rc3.txt for fast lookup
#
open(RC3," < catalogs/rc3.txt")|| die "Can not read rc3.txt\n";
$n_rc3 = 0;
while (<RC3>) {
    $pgc = &ss(106,116);
    $pgc =~ s/^PGC\s*0*(\d+)\s*$/$1/g;

    #
    # Galaxy name in RC3 as tag
    #--------------------------
    #
    $n1 = &ss(63,74);
    $altname = &ss(75,89);
    # squeeze out all spaces
    $n1 =~ s/ //g;
    $altname =~ s/ //g;
    $altname =~ s/UGC(\d+)/UGC\ $1/g;
    $altname =~ s/ESO(\d+-\w+)/ESO\ $1/g;
    $altname =~ s/CGCG0*(\d{1,3})-0*(\w+)/ZWG\ $1\.$2/g;
    $altname =~ s/MCG-?(\w+)/MCG\ $1/g;
    $name = "";
    if ($n1 =~ /NGC|IC|^A/) {
    $name = $n1;
    # compactify name
    $name = &cmpname ($name);
    } elsif ($n1 =~ /LMC|SMC|Sculptor|Fornax|WLM/) {
    $name = $n1;
    $altname = "NGC 292" if ( $name =~ /SMC/);
    } elsif ($n1 =~ /Leo(\w+)\b/) {
    $name = "Leo $1";
    } elsif (!$n1) {
    # empty string?
    $name = $altname;
    }
    #print "$name    $altname\n";
    # drop galaxies without name
    next if (! $name);
    $name_RC3{$pgc} = $name;

    if ($name =~ /LMC|SMC|Sculptor|Fornax|WLM|Leo\ I\b/) {

    $altname1{$name} = $altname;

    #
    # J2000 coordinates
    #-------------------
    #
    # RA
    $rh = &ss(1,2);
    $rm = &ss(3,4);
    $rs = &ss(5,8);

    # convert to hours

    $alpha{$name} = &ra2h_rc3 ($rh,$rm,$rs);

    # DEC
    $dd = &ss(11,12);
    $dm = &ss(13,14);
    $ds = 0;
    $ds = &ss(15,16) if (&ss(15,16) =~ /[0-9]/);
    $dsign = &ss(10,10);

    # convert to degrees

    $delta{$name} = &dec2deg_rc3 ($dsign,$dd,$dm,$ds);
    #
    # maj, min size  [arcmins]
    #-------------------------
    #
    $log10D25 = &ss(152,155);
    $log10D25 =~ s/ //g;
    $D = (10**$log10D25)/10; # units -> [0.1 arcmins]!
    $d = "";
    $d = $D / (10 ** &ss(162,165)) if (&ss(162,165) =~ /[0-9]/);
    #

    # Radius [ly] and inclination [degrees]
    #-----------------------------------------
    #

    $radius{$name} = 0.5*deg2rad($D/60)*$distance_B200{$name};
    if ($d){
        $i{$name} = &inclination($d/$D);
    } else {
        # equate min size $d to to max size, if $d lacking
        $d = $D;
        $i{$name} = 3.0;
    }
    #
    # PA, degrees E of N
    #---------------
    #
    $PA{$name} = &ss(186,188);
    $PA{$name} =~ s/ //g;
    # assign to 0 if undefined
    $PA{$name} = 0 if (! $PA{$name});
        #
        # apparent RC3 magnitudes
        #--------------------------
        #
        $Bmag_RC3{$name} = &ss(190,194);
        $Bmag_RC3{$name} =~ s/ //g;

    }

    #
    # RC3 distances, [ly]
    #--------------------
    # from Hubble's law, d ~ v_recession/H0,
    # if velocity >~ 1000 km/sec (positive!)
    # use weighted mean recession velocity corrected to the reference frame
    # defined by the 3K microwave background radiation, units [km/sec]
    $vrad = &ss(359,363);
    if ($vrad =~ /[0-9]/) {
        $n_rc3++ if ($vrad > 80.0);
    $rc3_distance{$name} = $vrad/$H0 * 1.0e6*$pc2ly;
    $v_CMB_RC3{$name} = $vrad;
    }
}
close (RC3);
printf "%5d RC3 galaxy distance data\n",$n_rc3;

$ngcic_flag = 0;        # =1 for NGC =2 for IC
$count =0;
$c_sbf = 0;
$c_hubble = 0;
$c_tf = 0;
$c_B200 = 0;

open(NGCIC,"<catalogs/revngcic.txt")|| die "Can not read revngcic.txt\n";
while (<NGCIC>) {
    $pgc = &ss(89,94);
    $pgc =~ s/ //g;
    $pgc =~ s/0*(\d+)/$1/g;
  #
  # Galaxy name
  #------------
  #
  # skip header and identify NGC or IC catalogs
  if (/^NGC/) {
    $ngcic_flag = 1;
    next;
  } elsif (/^\ IC/) {
    $ngcic_flag = 2;
    next;
  } elsif (/^-/) {
    next;
  }
  # skip extension lines
  next if (&ss(9,10) =~ /2|10/);

  # NGC number including C(omponent)-column
  $name = &ss(1,6);
  # get rid of extra spaces
  $name =~ s/^\s*//;
  $name =~ s/\s*$//;
  $name =~ s/\s+/ /g;
  # connect the different components (C-column) of an object with a dash
  $name =~ s/\s(\d)$/-$1/g;
  if ($ngcic_flag == 1) {
    $name = "NGC ".$name;
  } elsif ($ngcic_flag == 2) {
    $name = "IC ".$name;
  } else {
    next;
  }

  # status ( reading out the S-column)
  $status = &ss(11,12);
  $status =~ s/^\s*//;
  # drop all but galaxies S=1 (drop also  parts of galaxies S=6)
  # print "$name\n" if ($status =~ /6/);
  next if ($status =~ /[234567890]/);

  # NGC 292 = SMC will be taken care of separately
  next if($name =~ /NGC 292/);

  #
  # J2000 coordinates
  #-------------------
  #
  $ra = &ss(21,30);
  $dec = &ss(32,40);

  # Convert coordinates into hours (ra) and degrees (dec):

  $alpha = &ra2h ($ra);
  $delta = &dec2deg ($dec);
  #
  #
  # Galaxy type, translation to simplified Hubble classes
  #------------------------------------------------------
  #
  $Type = &ss(79,87);
  $Type =~ s/^\s*//;
  $Type =~ s/\s*$//;
  $Type =~ s/\s+/ /g;
  if ($Type =~ /S[\+\?0abcdm]?|^[ES]$|SB[?0abcdm]?|d?E[\+\?0-7]?|Irr|IB?m?|C|RING|Ring|GxyP|P|D|^R/) {
    $HType = $Type;