#!/usr/bin/perl
#
# program za analizo spremenljivih zvezd
# pozeni tako:
#
#   stars.pl  *.dat
#
#
# september 2003,  B.Dintinjana & skupina studentov
#    izboljsana verzija: ce ni lokalnega kataloga vzame imena zvezd iz dat
#

$i=0;  # index za stevilo slik
#
# vska slika ima:
#
# $time             cas meritve
# $totalflux        sestevek primerjalnih zvezd
# $zvezde[0..j]
#
# primerjalne zvezde
$comp1 = "GSC_2004-01176";  # izberem po obcutku!
$comp2 = "TYC2_2004-01627-1";
$comp3 = "TYC2_2004-01230-1";

$stcomp=3;



# okolica zvezde
$eps=0.4/3600;  #/


# za M3 nimam:
#
# Imam katalog m67.2000.lst z imeni zvezd in koordinatami
# Katalog preberem
#read_local_katalog();
#
# Imam datoteko m67.rem.txt z imeni zvezd in referenacami iz literature
# preberem
#read_local_ref();

open (MAG,">magnitude.dat");
print MAG "      STAR_NAME              JD                 MAG          ERR \n";
# hocem izpis naprimer takole:
#           FMF_J085152.9+114636  2452704.51292824        4.805        0.032

for $a (@ARGV)
{
    printf "#Datoteka $a\n";
    $totalflux[$i] = 0,
    $j=0; # index po zvezdah v tej sliki
    $stc=0; # stevilo primerjalnih zvezd na sliki
    undef @ime;
    open (STA, $a) ;
    while (<STA>){

        if(/^Input/){
            # prva vrstica v datoteki, preberem cas
            ($n1,$n2,$n3,$n4,$n5,$n6) = split;
            $n5 =~ s/time://g;
            #printf "$n5\n\n";
            $time[$i] = $n5;
            next;
        }
        ($n,$x,$y,$flux[$j],$ra,$dec,$cat,$name) = split;

        $ime[$j] = find_star($ra,$dec,$name);

        if($ime[$j] eq $comp1 || $ime[$j] eq $comp2 || $ime[$j] eq $comp3 ) {
            $totalflux[$i] += $flux[$j];
            # stejem primerjalne zvezde
            $stc++;
            #print "Total flux: ", $totalflux[$i],"\n";
            print "Ime slike: $a Ime primerjalne zvezde: ", $ime[$j],"\n";
        }
        #printf " %20s %12.3f \n",$ime[$j],$flux[$j];
        $j++;
    }
    close (STA);

    #
    # izpis izmerjene magnitude za vse zvezde na tej sliki
    #
    if($stc != $stcomp ){
        die "Na sliki $a ni primerjalne zvezde!\n";
    }

    for($j=0;$j<$#ime;$j++){
        $mag = -2.5 * log( $flux[$j] / $totalflux[$i] )/log(10.0);
        $err = sqrt($flux[$j])/$flux[$j];
        printf MAG " %-20s %17.8f %12.3f %12.3f \n",$ime[$j],$time[$i],$mag,$err;

    }
    # naslednja slika /
    $i++;
}
close MAG;
#
# na koncu zapisem se seznam zvezd in dodam se navzkrizne reference
# iz baze simbad
#
open (REF,">reference.dat");
printf REF "# COORDINATE LIST\n";
printf REF "# Numbering system for M3.\n";
printf REF "#\n";
printf REF "#\n";
printf REF "#    STAR NAME          RA(J2000)   DEC(J2000)     N   REF.\n";

for ($k=0;$k<=$#seznam_ime;$k++){

    # samo zvezde z vec kot .. meritev
    if($seznem_n[$k] > 5){
        # zakomentiraj ce gre prepocasi
        # poizvedba v internetu iz baze simbad, hekersko
        #$znak=($seznam_dec[$k]>0)?"%2B":"";
        #$query=`lynx -dump \"http://simbad.u-strasbg.fr/sim-id.pl?protocol=html&Ident=$seznam_ra[$k]+$znak$seznam_dec[$k]&Radius=1&Radius.unit=arcsec\" | grep \"|\"| grep \"^\\[\"  |sed -e \"s/   //g\" `;
        #$query=~ s/\n//g;

        printf REF "%-20s %12s %12s %5d  %-30s %-50s\n",$seznam_ime[$k],dectohms1($seznam_ra[$k]/15,":"),dectodms1($seznam_dec[$k],":"),$seznem_n[$k],$seznam_rem[$k],$query;
        printf     "%-20s %12s %12s %5d  %-30s %-50s\n",$seznam_ime[$k],dectohms1($seznam_ra[$k]/15,":"),dectodms1($seznam_dec[$k],":"),$seznem_n[$k],$seznam_rem[$k],$query;
    }
}
close REF;








sub find_star
{
    my ($ra,$dec,$name) = @_;
    #
    # najprej preverim ce je zvezda ze v mojem seznamu
    # indes k tece po seznamu od 0 .. k
    #
    for ($k=0;$k<=$#seznam_ime;$k++){
        if(abs($ra-$seznam_ra[$k]) < $eps && abs($dec-$seznam_dec[$k]) < $eps) {
            # nasel zvezdo v seznamu
            # stejem meritv
            $seznem_n[$k]++;
            return $seznam_ime[$k];
        }
    }
    #
    # nato pogledam ce je zvezda v lokalnem katalogu
    # index l tece po katalogu, index k pa kaze na prazen (nov) zapis v seznamu
    # ce katakoga ni, gre kar naprej
    #
    for ($l=0;$l<=$#katalog_ime;$l++){
        if(abs($ra-$katalog_ra[$l]) < $eps && abs($dec-$katalog_dec[$l]) < $eps) {
            #
            # nasel zvezdo v lokalnem katalogu
            # pogledam se ce je kaksna referenca
            for ($r=0;$r<=$#rem_ime;$r++){
                if($katalog_ime[$l] eq $rem_ime[$r]){
                    $seznam_rem[$k] = $rem_rem[$r];
                    last;
                }
            }
            # zapisem v moj seznam referenc
            # index k kaze na zadnji +1
            $seznam_ime[$k] = $katalog_ime[$l];
            $seznam_ra [$k] = $ra;
            $seznam_dec[$k] = $dec;
            $seznem_n  [$k] = 1;
            return $seznam_ime[$k];
        }
    }
    # zvezde ni nasel, torej nov vpis v seznam
    # index k kaze na zadnji+1
    if($name ne ""){
        # ce zvezda ze ima ime iz dat datoteke
        $seznam_ime[$k] = $name;
        $seznam_ra[$k]  = $ra;
        $seznam_dec[$k] = $dec;
        $seznem_n  [$k] = 1;
    }else{
        #sicer naredim novo ime
        $seznam_ime[$k] = new_star_name($ra,$dec);
        $seznam_ra[$k]  = $ra;
        $seznam_dec[$k] = $dec;
        $seznem_n  [$k] = 1;
    }
    return $seznam_ime[$k];
}

sub new_star_name
{
    # naredim ime po IAU priporocilih, standardu
	my ($ra,$dec) = @_;
	my $name = "FMF_J" . dectohms1($ra/15,"") . dectodms1($dec,"");
	return $name;
}

sub read_local_katalog
{
   open(LST,"local.lst") || die "cannot open file local.lst !\n";
   $i=0;
   while (<LST>){
      if(/^#/) {next}
      ($katalog_ime[$i],$ra_lst[$i],$dec_lst[$i])=split;
      $katalog_ra [$i] = strtora($ra_lst[$i]) /3600.0;
      $katalog_dec[$i] = strtodec($dec_lst[$i]) /3600.0;
      #printf "%12s\t%11s\t%11s\n", $katalog_ime[$i],$ra_lst[$i],$dec_lst[$i];
      $i++;
   }
   printf "Read $i stars from local catalog.\n";
   close (LST);
}

sub read_local_ref
{
    open(LST,"m67.rem.txt") || die "cannot open file m67.rem.txt !\n";
    $r=0;
    while (<LST>){
        if(/^#/) {next}
        ($rem_ime[$r],$rem_rem[$r])=split("\t");
        $rem_rem[$r] =~ s/\n//g;
        $r++;
   }
   printf "Read $r stars from local reference.\n";
   close (LST);
}


#  Transform a string to right ascension in arseconds
sub strtora

{
    my ($ras) = @_;
    my ($h, $m, $s);
    ($h,$m,$s) = split(":",$ras);
    return 15 * ($h * 3600.0 + $m * 60.0 + $s);
}

#  Transform a string to declination in arseconds
sub strtodec

{
    my ($decs) = @_;
    my ($d, $m, $s)=split(":",$decs);
    if(index($decs,"-") > -1){
      return -1.0*(abs($d) * 3600.0 + abs($m) * 60.0 + abs($s));
    }
    else {
       return (abs($d) * 3600.0 + abs($m) * 60.0 + abs($s));
    }

}

#  Spremeni stevilo iz decimalnega v format hhmmss.ss
#  Podano je se locilo
sub dectohms
{
    my ($x, $sep) = @_;
    my ($h, $m, $s);

    $h = int($x);
    $x -= $h;
    $x *= 60.0;
    $m = int($x + 1e-9);
    $s = ($x - $m + 1e-9) * 60;
    return sprintf("%02d%s%02d%s%05.2f", $h, $sep, $m, $sep, $s);
}



#  Spremeni stevilo iz decimalnega v format ddmmss.s
#  Podano je se locilo
sub dectodms
{
    my ($x, $sep) = @_;
    my ($h, $m, $s, $c);

    if ($x < 0) {
	$c = '-';
	$x = -$x;
    }
    else {
	$c = '+';
    }
    $h = int($x);
    $x -= $h;
    $x *= 60.0;
    $m = int($x + 1e-9);
    $s = ($x - $m + 1e-9) * 60;
    return sprintf("%s%02d%s%02d%s%04.1f", $c, $h, $sep, $m, $sep, $s);
}



sub dectodms1
{
    my ($x, $sep) = @_;
    my ($h, $m, $s, $c);

    if ($x < 0) {
	$c = '-';
	$x = -$x;
    }
    else {
	$c = '+';
    }
    $h = int($x);
    $x -= $h;
    $x *= 60.0;
    $m = int($x + 1e-9);
    $s = ($x - $m + 1e-9) * 60;
    return sprintf("%s%02d%s%02d%s%02.0f", $c, $h, $sep, $m, $sep, $s);
}

sub dectohms1
{
    my ($x, $sep) = @_;
    my ($h, $m, $s);

    $h = int($x);
    $x -= $h;
    $x *= 60.0;
    $m = int($x + 1e-9);
    $s = ($x - $m + 1e-9) * 60;
    return sprintf("%02d%s%02d%s%04.1f", $h, $sep, $m, $sep, $s);
}

