#!/usr/bin/perl -w

use strict;
use lib qw(../..);
use Experiment;
# GnuplotIF doesn't seem very useful except by using
# gnuplot commands directly.
use GnuplotIF;
use Getopt::Long;

# default is to run simulations and can also be set with --run
# won't run simulations if --norun specify on command line.
my $runSimulations = 1;
# default is to do all the gates
my $ICaLgates = 0;
my $IKrgates = 0;
my $IKsgate = 0;
my $INagates = 0;
my $Itogates = 0;
my $periodic = 1;
GetOptions("run!"=>\$runSimulations,
           "ICaLgates!"=>\$ICaLgates,
           "IKrgates!"=>\$IKrgates,
           "IKsgate!"=>\$IKsgate,
           "INagates!"=>\$INagates,
           "Itogates!"=>\$Itogates,
           "periodic!"=>\$periodic
           );

my $e = new Experiment;
$e->promptWarnings(1);

sub ICaLgates {
    my $modelC = "ICaL_gates-calcium_clamp.xml";
    my $dataC = "data/ICaL_gates-calcium_clamp.data";
    if ($runSimulations) {
        if ($e->makeDirectory($dataC) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($modelC,$dataC) < 0) {
            die "Error simulating the model";
        }
    }
    my $modelV = "ICaL_gates-voltage_clamp.xml";
    my $dataV = "data/ICaL_gates-voltage_clamp.data";
    if ($runSimulations) {
        if ($e->makeDirectory($dataV) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($modelV,$dataV) < 0) {
            die "Error simulating the model";
        }
    }

    my $plot = GnuplotIF->new();
    $plot->gnuplot_cmd("set multiplot");
    # Figure 2(a) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$dInf) = $e->getDataIndices($dataV,"voltage_clamp/Vm",
                                            "d_gate/d_infinity");
        $plot->gnuplot_cmd(
      "plot [-50:30] \"$dataV\" using $Vm:$dInf title \"d_inf\" with lines");
    }
    # Figure 2(b) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$fInf) = $e->getDataIndices($dataV,"voltage_clamp/Vm",
                                            "f_gate/f_infinity");
        $plot->gnuplot_cmd(
      "plot [-80:40] \"$dataV\" using $Vm:$fInf title \"f_inf\" with lines");
    }
    # Figure 2(c) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.0,0.0");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Cai,$fCa) = $e->getDataIndices($dataC,"calcium_clamp/Cai",
                                            "fCa_gate/fCa_infinity");
        $plot->gnuplot_cmd(
          "plot \"$dataC\" using $Cai:$fCa title \"fCa_inf\" with lines");
    }
    # Figure 2(d) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0.25");
        $plot->gnuplot_cmd("set size 0.5,0.25");
        my ($Vm,$tauD) = $e->getDataIndices($dataV,"voltage_clamp/Vm",
                                            "d_gate/tau_d");
        $plot->gnuplot_cmd(
      "plot [-80:40] \"$dataV\" using $Vm:$tauD title \"tau_d\" with lines");
    }
    # Figure 2(e) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0.0");
        $plot->gnuplot_cmd("set size 0.5,0.25");
        my ($Vm,$tauF) = $e->getDataIndices($dataV,"voltage_clamp/Vm",
                                            "f_gate/tau_f");
        $plot->gnuplot_cmd(
      "plot [-100:50] \"$dataV\" using $Vm:$tauF title \"tau_f\" with lines");
    }
    $plot->gnuplot_cmd("unset multiplot");
    $plot->gnuplot_pause(-1,"Figure 2 from ten Tusscher et al (2004)");
}

sub IKrgates {
    my $model = "IKr_gates-voltage_clamp.xml";
    my $data = "data/IKr_gates-voltage_clamp.data";
    if ($runSimulations) {
        if ($e->makeDirectory($data) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($model,$data) < 0) {
            die "Error simulating the model";
        }
    }
    my $plot = GnuplotIF->new();
    $plot->gnuplot_cmd("set multiplot");
    # Figure 6(a) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0.5");
        $plot->gnuplot_cmd("set size 1,0.5");
        my ($Vm,$X1,$X2) = 
            $e->getDataIndices($data,"voltage_clamp/Vm",
                               "Xr1_gate/Xr1_infinity",
                               "Xr2_gate/Xr2_infinity");
        $plot->gnuplot_cmd(
          "plot [-150:50] \"$data\" using $Vm:$X1 title \"Xr1\" with lines,"
                        ."\"$data\" using $Vm:$X2 title \"Xr1\" with lines"
                           );
    }
    # Figure 6(b) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "Xr1_gate/tau_Xr1");
        $plot->gnuplot_cmd(
       "plot [-150:50] \"$data\" using $Vm:$tau title \"tau_Xr1\" with lines");
    }
    # Figure 6(c) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "Xr2_gate/tau_Xr2");
        $plot->gnuplot_cmd(
       "plot [-150:50] \"$data\" using $Vm:$tau title \"tau_Xr2\" with lines");
    }
    $plot->gnuplot_cmd("unset multiplot");
    $plot->gnuplot_pause(-1,"Figure 6 from ten Tusscher et al (2004)");
}

sub IKsgate {
    my $model = "IKs_gate-voltage_clamp.xml";
    my $data = "data/IKs_gate-voltage_clamp.data";
    if ($runSimulations) {
        if ($e->makeDirectory($data) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($model,$data) < 0) {
            die "Error simulating the model";
        }
    }
    my $plot = GnuplotIF->new();
    $plot->gnuplot_cmd("set multiplot");
    # Figure 5(a) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0.5");
        $plot->gnuplot_cmd("set size 1,0.5");
        my ($Vm,$Xs) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                          "Xs_gate/Xs_infinity");
        $plot->gnuplot_cmd(
          "plot [-100:100] \"$data\" using $Vm:$Xs title \"Xs_inf\" with lines"
                           );
    }
    # Figure 5(b) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0");
        $plot->gnuplot_cmd("set size 1,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "Xs_gate/tau_Xs");
        $plot->gnuplot_cmd(
       "plot [-50:50] \"$data\" using $Vm:$tau title \"tau_Xs\" with lines");
    }
    $plot->gnuplot_cmd("unset multiplot");
    $plot->gnuplot_pause(-1,"Figure 5 from ten Tusscher et al (2004)");
}

sub INagates {
    my $model = "INa_gates-voltage_clamp.xml";
    my $data = "data/INa_gates-voltage_clamp.data";
    if ($runSimulations) {
        if ($e->makeDirectory($data) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($model,$data) < 0) {
            die "Error simulating the model";
        }
    }
    my $plot = GnuplotIF->new();
    $plot->gnuplot_cmd("set multiplot");
    # Figure 1(a) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$m) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                         "m_gate/m_infinity");
        $plot->gnuplot_cmd(
          "plot [-80:20] \"$data\" using $Vm:(\$$m*\$$m*\$$m) "
            ."title \"m_inf**3\" with lines"
                           );
    }
    # Figure 1(b) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($Vm,$h,$j) = 
            $e->getDataIndices($data,"voltage_clamp/Vm","h_gate/h_infinity",
                               "j_gate/j_infinity");
        $plot->gnuplot_cmd(
                           "plot [-150:0] \"$data\" using $Vm:(\$$h*\$$j) "
                           ."title \"h*j\" with lines");
    }
    # Figure 1(c) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.0,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "m_gate/tau_m");
        $plot->gnuplot_cmd(
                           "plot [-100:40] \"$data\" using $Vm:$tau "
                           ."title \"tau_m\" with lines");
    }
    # Figure 1(d) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.33,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "h_gate/tau_h");
        $plot->gnuplot_cmd(
                           "plot [-120:0] \"$data\" using $Vm:$tau "
                           ."title \"tau_h\" with lines");
    }
    # Figure 1(e) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.66,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($Vm,$tau) = $e->getDataIndices($data,"voltage_clamp/Vm",
                                           "j_gate/tau_j");
        $plot->gnuplot_cmd(
                           "plot [-130:0] \"$data\" using $Vm:$tau "
                           ."title \"tau_j\" with lines");
    }
    $plot->gnuplot_cmd("unset multiplot");
    $plot->gnuplot_pause(-1,"Figure 1 from ten Tusscher et al (2004)");
}

sub Itogates {
    my $modelN = "Ito_gates-voltage_clamp-endo.xml";
    my $dataN = "data/Ito_gates-voltage_clamp-endo.data";
    if ($runSimulations) {
        if ($e->makeDirectory($dataN) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($modelN,$dataN) < 0) {
            die "Error simulating the model";
        }
    }
    my $modelP = "Ito_gates-voltage_clamp-epi.xml";
    my $dataP = "data/Ito_gates-voltage_clamp-epi.data";
    if ($runSimulations) {
        if ($e->makeDirectory($dataP) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($modelP,$dataP) < 0) {
            die "Error simulating the model";
        }
    }

    my $plot = GnuplotIF->new();
    $plot->gnuplot_cmd("set multiplot");
    # Figure 3(a) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($VmN,$rN) = $e->getDataIndices($dataN,"voltage_clamp/Vm",
                                           "r_gate/r_infinity");
        my ($VmP,$rP) = $e->getDataIndices($dataP,"voltage_clamp/Vm",
                                           "r_gate/r_infinity");
        $plot->gnuplot_cmd(
          "plot [-80:80] \"$dataN\" using $VmN:$rN "
            ."title \"r_inf-endo\" with lines,"
            ."\"$dataP\" using $VmP:$rP "
            ."title \"r_inf-epi\" with points"
                           );
    }
    # Figure 3(b) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.5,0.5");
        $plot->gnuplot_cmd("set size 0.5,0.5");
        my ($VmN,$sN) = $e->getDataIndices($dataN,"voltage_clamp/Vm",
                                           "endo_s_gate/s_infinity");
        my ($VmP,$sP) = $e->getDataIndices($dataP,"voltage_clamp/Vm",
                                           "epi_M_s_gate/s_infinity");
        $plot->gnuplot_cmd(
          "plot [-80:40] \"$dataN\" using $VmN:$sN "
            ."title \"s_inf-endo\" with lines,"
            ."\"$dataP\" using $VmP:$sP "
            ."title \"s_inf-epi\" with points"
                           );
    }
    # Figure 3(c) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.0,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($VmN,$tN) = $e->getDataIndices($dataN,"voltage_clamp/Vm",
                                           "r_gate/tau_r");
        my ($VmP,$tP) = $e->getDataIndices($dataP,"voltage_clamp/Vm",
                                           "r_gate/tau_r");
        $plot->gnuplot_cmd(
          "plot [-100:65] \"$dataN\" using $VmN:$tN "
            ."title \"tau_r-endo\" with lines,"
            ."\"$dataP\" using $VmP:$tP "
            ."title \"tau_r-epi\" with points"
                           );
    }
    # Figure 3(d) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.33,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($Vm,$tau) = $e->getDataIndices($dataP,"voltage_clamp/Vm",
                                           "epi_M_s_gate/tau_s");
        $plot->gnuplot_cmd(
                           "plot [-100:40] \"$dataP\" using $Vm:$tau "
                           ."title \"tau_s-epi\" with lines");
    }
    # Figure 1(e) of ten Tusscher (2004) article
    {
        $plot->gnuplot_cmd("set origin 0.66,0.0");
        $plot->gnuplot_cmd("set size 0.33,0.5");
        my ($Vm,$tau) = $e->getDataIndices($dataN,"voltage_clamp/Vm",
                                           "endo_s_gate/tau_s");
        $plot->gnuplot_cmd(
                           "plot [-100:50] \"$dataN\" using $Vm:$tau "
                           ."title \"tau_s-endo\" with lines");
    }
    $plot->gnuplot_cmd("unset multiplot");
    $plot->gnuplot_pause(-1,"Figure 3 from ten Tusscher et al (2004)");
}

sub periodic {
    my $endoM = "periodic-stimulus-endocardial.xml";
    my $endoD = "data/periodic-stimulus-endocardial.data";
    if ($runSimulations) {
        if ($e->makeDirectory($endoD) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($endoM,$endoD) < 0) {
            die "Error simulating the model";
        }
    }
    my $epiM = "periodic-stimulus-epicardial.xml";
    my $epiD = "data/periodic-stimulus-epicardial.data";
    if ($runSimulations) {
        if ($e->makeDirectory($epiD) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($epiM,$epiD) < 0) {
            die "Error simulating the model";
        }
    }
    my $midM = "periodic-stimulus-M.xml";
    my $midD = "data/periodic-stimulus-M.data";
    if ($runSimulations) {
        if ($e->makeDirectory($midD) < 0) {
            die "Error making data directory";
        }
        if ($e->simulate($midM,$midD) < 0) {
            die "Error simulating the model";
        }
    }

    # the full trace
    {
        my $plot = GnuplotIF->new();
        my ($epiT,$epiV) = $e->getDataIndices($epiD,"time",
                                              "membrane_potential/V");
        my ($endoT,$endoV) = $e->getDataIndices($endoD,"time",
                                                "membrane_potential/V");
        my ($midT,$midV) = $e->getDataIndices($midD,"time",
                                              "membrane_potential/V");
        $plot->gnuplot_cmd(
                           "plot "
                           ."\"$endoD\" using $endoT:$endoV "
                             ."title \"endo\" with lines,"
                           ."\"$epiD\" using $epiT:$epiV "
                             ."title \"epi\" with lines,"
                           ."\"$midD\" using $midT:$midV "
                             ."title \"M\" with lines"
                           );
        $plot->gnuplot_pause(-1,"Figure 8(a) from my Noble review");
    }

    # something like figure 8(a) from my Noble review
    {
        my $plot = GnuplotIF->new();
        my ($epiT,$epiV) = $e->getDataIndices($epiD,"time",
                                              "membrane_potential/V");
        my ($endoT,$endoV) = $e->getDataIndices($endoD,"time",
                                                "membrane_potential/V");
        my ($midT,$midV) = $e->getDataIndices($midD,"time",
                                              "membrane_potential/V");
        $plot->gnuplot_cmd(
                           "plot [19.9e3:20.4e3] "
                           ."\"$endoD\" using $endoT:$endoV "
                             ."title \"endo\" with lines,"
                           ."\"$epiD\" using $epiT:$epiV "
                             ."title \"epi\" with lines,"
                           ."\"$midD\" using $midT:$midV "
                             ."title \"M\" with lines"
                           );
        $plot->gnuplot_pause(-1,"Figure 8(a) from my Noble review");
    }
}

if ($ICaLgates) {
    &ICaLgates();
}
if ($IKrgates) {
    &IKrgates();
}
if ($IKsgate) {
    &IKsgate();
}
if ($INagates) {
    &INagates();
}
if ($Itogates) {
    &Itogates();
}
if ($periodic) {
    &periodic();
}


__END__