#!/usr/bin/perl
use Math::Integral::Romberg 'integral';
binmode STDOUT, ":utf8";

#14.99(wavelength in meters)*2.62x10^-13(permitivity of free space,etc) * integral (5000001(electrons per cubic meter))*(10x10^-12(nanotesla magnetic field)) dx from x=0 to 149597871475.30 (meters from sun to earth)
# wolframalpha:14^2 * (2.62*10^-13) * integrate ((5*10^6)*(10*10^-12)) dx from x=0 to 149597871475.30
# == 4.40352*10^6

my $speedoflight = 299792458;
my $freespace = 2.62*10**-13;#electroncharge^3/(8*pi*permitivityoffreespace*massofelectron*speedoflight)
my $electrondensity = 5*10**6; # electrons/meter^3 , about 5 electrons /cm^3
my $bfield = 10*10**-12; # 10 nanotesla
my $distancetoearth = 149597871475.30; # meters

print "(8*pi*permitivityoffreespace*massofelectron*speedoflight): $freespace\nB-field: $bfield tesla\nElectron Density: $electrondensity electrons/meter^3\n";

sub electrondensity { my ($x) = @_; ((1*10**6)*((2.39*10**8/$x^6)+(1.67*10**6/$x^2.3))) }
sub f { my ($x) = @_; ($bfield + $electrondensity) }
#sub f { my ($x) = @_; ($lump) }

my $rotationmeasure = $freespace*(integral(\&f, 0.001, $distancetoearth));
#can't use Math::Integral::Romberg because it sucks fucking donkey balls and returns bogus shit
#my $rotationmeasure = $freespace*integrate { $lump } from => 0.001, to => $distancetoearth;
print "\nfrequency | \x{3BB} meters | radians faraday rotation\n";
for (0.001,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1..20) { #wavelength squared
	#print (integral(\&f, 0, $distancetoearth));
	printf("%1e :",$speedoflight/$_);
	printf("%2g\t", $_);
	printf("%1g", ($_**2)*$rotationmeasure);
	print "\n";	
}
print "\nTotal electrons in path: ";
printf("%1g", integral(\&electrondensity, 0.01, 1));

sub integrate(&@) { #prototypes can blow me.
    local $x;
    my $sum;
    my $f = shift;
    my %o = (from => 0, by => 0.01, @_);
    for ($o{from} / $o{by} .. $o{to} / $o{by}) {
        $x = $_ * $o{by};
        $sum += &$f * $o{by};
    }
    $sum;
}