#!/usr/bin/perl -w
use strict;

################## INITIALIZATION ###################
# Reaction to be modeled is S --> P
my $S = 1;
my $P = 0;
my $k = 1;
my $dt = .01;
my @slopes;
my @values = ($S,$P);
my $time_steps = 500;

###################### FILES ########################

my $output_path = 'ADPAs.txt';
open OUTPUT_FILE, ">$output_path" or die "Can't open $output_path: $!\n";

################### MAIN PROGRAM ####################

Output_Header();
Output(@values);
foreach my $i (1 .. $time_steps) {
   @slopes = model(@values);
   my @increments = multiply($dt,@slopes);
   @values = add(\@values,\@increments);
   Output(@values);  
}

#################### SUBROUTINES ####################
sub model {
#  Reaction is S --> P
  my ($S,$P) = @_;
  my $S_prime = -$k * $S;
  my $P_prime =  $k * $S;
  return ($S_prime, $P_prime);
}

sub Output {
  print join("\t", @_), "\n";
  print OUTPUT_FILE join("\t",@_), "\n";
}

sub Output_Header {
  print OUTPUT_FILE "Substrate", "\t", "Product", "\n";
}

sub multiply {
  my ($scalar,@array) = @_;
  my @product;
  foreach my $i (0 .. @array - 1) {
     $product[$i] = $scalar * $array[$i];
  }
  return @product;
}

sub add {
  my ($array1ref,$array2ref) = @_;
  my @array1 = @$array1ref;
  my @array2 = @$array2ref;
  my @sum;
  foreach my $i (0 .. @array1 - 1) {
     $sum[$i] = $array1[$i] + $array2[$i];
     if ($sum[$i] < 0) {$sum[$i] = 0}
  }
  return @sum;
}
