#!/usr/bin/perl -w

use strict;
use File::Basename;
use Getopt::Std;

# parsing commandline options

my (%opts);
getopts('f:', \%opts);

my ($seqname,$opts, $seq, %seqs, $seqlen, $i, $as, $horiz, %count, $entropie,   );

readseq();

# $seqs{$seqname} = $seq;
#  $seqs{7s} = gaucgaucgaucgaucga;
#  $seqs{5s} = acgaucgaugcagc;
#  $seqs{meine} = aucg;
# 
for ($horiz = 0; $horiz < $seqlen; $horiz++)
{

    # initialisieren: alles 0 setzen
    
    foreach $seqname (keys %seqs)   #vertikal: 7s, 5s, meine
    {  
    
        # zeile: seqname, spalte: $horiz
        $as = substr($seqs{$seqname}, $horiz, 1);
         $count{$as} += 1;
        print "dies ist $as in zeile $seqname spalte $horiz\n";
        print "zaehle count{$as} +1 =$count{$as}\n";
    }
    print "fertig mit spalte $horiz: jetzt rechnen\n";
    
# subroutines
#sub log 10 {
 # my $n = shift; 
 # return log($n)/log(10); 
#}
      
# Entropie
#    foreach ($count {$as} / 0.05)
    foreach $as (keys %count);     #ZEILE:50 FEHLERMELDUNG!!!!!
                                   # für jede aminosäure
    {
        #my $entropie=$count{$as} * ((log($n)/log(2)) * ($count{$as} / 0.05));
        my $entropie = $count{$as} * log2($count{$as} / 0.05);
        print "In der Spalte betraegt der Wert: $entropie fr $as\n";
    }

}


# end main

###############################################################################
sub log2 {
##########
    my $n = shift; 
    return log($n)/log(2); 
}

###############################################################################
# read sequence from file or command line
sub readseq {
#############
   # from file
        open(FILE, "$opts{f}")
            or die "can't open file:!";
        while (<FILE>){
            if ($_ =~ /^>/) {
                chomp;
                ($seqname = $_) =~ s/^> *//;
            } elsif (/^[ACDEFGHIKLMNPQRSTVWYacdefghiklmnpqrstvwy-]+$/) {
                chomp;
                $seq = uc($_);
                $seqs{$seqname} = $seq;
        $seqlen=length($seq);
            } elsif (/^$/) {
                next;
            } else {
                print "Sequence \"$seqname\" is not nucleic acid!\n";
        exit -1;
            }
        }
        close FILE;
}
# end readseq