#!/usr/bin/perl
use strict;
use warnings;

# Perl solving a physics / electrodynamics problem involving
# symbolic mathematics, derivatives and complex numbers:

use Math::Symbolic qw/:all/;
use Math::Complex;

# Given the following simple circuit:
#
#  ----|||||-----/\/\/\----       (R = resistor,
# |      R          L      |       L = solenoid,
# |                        |       U = alternating voltage)
#  ---------O ~ O----------
#            U(t)
#
# Question: What's the current in this circuit?
#
# We'll need some physics before letting the computer do the
# math:
# Applying Kirchhoff's rules, one quickly ends up with the
# following differential equation for the current:
#     (L * dI/dt) + (R * I)  =  U

my $left  = parse_from_string('L * total_derivative(I(t), t) + R * I(t)');
my $right = parse_from_string('U(t)');

<readmore>
# If we understand current and voltage to be complex functions,
# we'll be able to derive. ("'" denoting complex here)
#    I'(t) = I'_max * e^(i*omega*t)
#    U'(t) = U_max  * e^(i*omega*t)
# (Please note that omega is the frequency of the alternating voltage.
# For example, the voltage from German outlets has a frequency of 50Hz.)

my $argument = parse_from_string('e^(i*omega*t)');
my $current  = parse_from_string('I_max') * $argument;
my $voltage  = parse_from_string('U_max') * $argument;

# Putting it into the equation:
$left->implement( I  => $current );
$right->implement( U => $voltage );

$left = $left->apply_derivatives()->simplify();

# Now, we can solve the equation to get a complex function for
# the current:

$left  /= $argument;
$right /= $argument;
my $quotient = parse_from_string('R + i*omega*L');
$left  /= $quotient;
$right /= $quotient;

# Now we have:
#    $left    = $right
#    I_max(t) = U_max / (R + i*omega*L)
# But I_max(t) is still complex and so is the right-hand-side of the
# equation!

# Making the symbolic i a "literal" Math::Complex i
$right->implement(
    e => Math::Symbolic::Constant->euler(),
    i => Math::Symbolic::Constant->new(i),    # Math::Complex magic
);

print <<'HERE';
Sample of complex maximum current with the following values:
  U_max => 100
  R     => 10
  L     => 10
  omega => 1
HERE

print "Computed to: "
  . $right->value(
    U_max => 100,
    R     => 10,
    L     => 10,
    omega => 1,
  ),
  "\n\n";

# Now, we're dealing with alternating current and voltage.
# So let's make a generator that generates nice current
# functions of time!
#   I(t) = Re(I_max(t)) * cos(omega*t - phase);

# Usage: generate_current(U_Max, R, L, omega, phase)
sub generate_current {
    my $current = $right->new();    # cloning

    $current *= parse_from_string('cos(omega*t - phase)');

    $current->implement(
        U_max => $_[0],
        R     => $_[1],
        L     => $_[2],
        omega => $_[3],
        phase => $_[4],
    );
    $current = $current->simplify();
    return sub { Re( $current->value( t => $_[0] ) ) };
}

print "Sample current function with: 230V, 2Ohms, 0.1H, 50Hz, PI/4\n";
my $current_of_time = generate_current( 230, 2, 0.1, 50, PI / 4 );

print "The current at 0 seconds:   " . $current_of_time->(0) . "\n";
print "The current at 0.1 seconds: " . $current_of_time->(0.1) . "\n";
print "The current at 1 second:    " . $current_of_time->(1) . "\n";

</readmore>