#!/usr/bin/perl
use strict;

# input file
my $input_file = "1msun.mdl";

# get user input-file from command line
my $i = 0;
while ($i <= $#ARGV) 
{
    if ($ARGV[$i] eq '-f' || 
	$ARGV[$i] eq '-file' || 
	$ARGV[$i] eq '-input' || 
	$ARGV[$i] eq '-input_file') {
	$input_file = $ARGV[++$i];
    }
    $i++;
}

# column number with...
my $col_r = 2; # radius r of mass shell
my $col_rho = 4; # density of the mass shell
my $col_u = 21; # internal energy of the mass shell
my $col_m = 1; # mass co-ordinate

# some constants (cgs)
use constant G => 6.67259E-8;
use constant Msun => 1.9891E33;
use constant Rsun => 6.9598E10;
use constant Lsun => 3.8515E33;
use constant stefan_boltzmann => 5.67051E-5;

# gravitational energy
my $E_grav = 0.0;

# internal energy
my $E_int = 0.0;

open(FP, "<$input_file") or die "Could not open $input_file: $!\n";
$i = 0;
my @x;
my $prevm;
while (<FP>)
{
    # ignore first line because it contains only information about the format of the file
    chomp;
    @x=split(/\s+/,$_);

    # lagrangian mass co-ordinate
    my $m=$x[$col_m]*Msun;

    if(defined($prevm))
    {
	# radius (cgs)
	
	my $radius = Rsun * $x[$col_r];

	# mass of the shell (cgs)
	my $dm = $m - $prevm;

	# internal energy of mass shell (erg/g is stored : so multiply by the shell mass)
	my $u = $dm * $x[$col_u]; 
	
	############################################################################
	# Insert the correct calculation of the internal and gravitational energy! #
	############################################################################
	# sum up the gravitational energy
	$E_grav -= G*$m*$dm/$radius;

	# sum up the internal energy
	$E_int += $u;
	############################################################################
	
	if($dm>0.0)
	{
	    printf "Line $i at M=%g Msun, R=%g Rsun, U=%g (dm=%g Msun) Eg=%g Ei=%g\n",$m/Msun,$radius/Rsun,$u,$dm/Msun,$E_grav,$E_int;
	}
    }

    $prevm=$m;
    $i++;
}
close(FP);

# Mass and radius of the star
my $M=$x[$col_m];
my $R=$x[$col_r];

# ratio of E_grav and E_int
my $ratio = -$E_grav/$E_int;

# factor GM^2/R
my $alpha = -$E_grav/(G*$M*$M/$R*Msun*Msun/Rsun);

print "Gravitational energy: $E_grav erg\n";
print "Internal energy: $E_int erg\n";
print "Ratio of energies: $ratio\n";
print "Alpha-factor of Egrav: $alpha\n";