with ConsAKD; use ConsAKD;
     with FloaAKD; use FloaAKD;
     with InteAKD; use InteAKD;
     with MathAKD; use MathAKD;
     with TextAKD; use TextAKD;

     package body CNS221 is
     ----------------------------------------------------------------------
     ----------------------------------------------------------------------

     procedure Cable_Passive_Steady_State (
       Rm       : in     float := Rm_Def;
       Ri       : in     float := Ri_Def;
       Diameter : in     float := Diameter_Def ) is
     ----------------------------------------------------------------------
       Cm       : constant float := 1.0e-6; -- F/cm2
       V0       : constant float := 1.0;    -- Volt
       X        :          float;
       X_Max    :          float;           -- m
       dX       :          float;
       Vm       :          float;
       Lambda   :          float;           -- space constant
       Rm_New   :          float := Rm;
       Ri_New   :          float := Ri;
       Di_New   :          float := Diameter;
     begin
       for index in 1..4 loop
	 if index = 2 then
	   Rm_New := Rm / 2.0;
	 elsif index = 3 then
	   Rm_New := Rm;
	   Ri_New := Ri / 2.0;
	 elsif index = 4 then
	   Ri_New := Ri;
	   Di_New := Diameter / 2.0;
	 end if;
	 X := 0.0;
	 Vm := 0.0;
	 Lambda := ( Rm_New * Di_New ) / ( Ri_New * 4.0 );
	 X_Max := 3.0 * Lambda;
	 dX := X_Max / 1000.0;
	 while X <= X_Max loop
	   Vm := V0 * Exp ( -abs ( X ) / Lambda );
	   Put (  X, 2, 6, 4 ); Put ( ' ' );
	   Put ( Vm, 2, 6, 4 ); New_Line;
	   X := X + dX;
	 end loop;
       end loop;
     end Cable_Passive_Steady_State;

     function Current_Injection (
       Current : Current_Type;
       Time    : float       ;
       T_Peak  : float := 0.5 )
       return float is
     ----------------------------------------------------------------------
     begin
       case Current is
	 when Impulse =>
	   if Time = 0.0 then
	     return 1.0;
	   else
	     return 0.0;
	   end if;
	 when Step =>
	   if Time >= 0.0 then
	     return 1.0;
	   else
	     return 0.0;
	   end if;
	 when Alpha =>
	   return Time * Exp ( -Time / T_Peak );
       end case;
     end;

     procedure Demo is
     ----------------------------------------------------------------------
       Min : constant natural := 0;
       Max : constant natural := 3;
       Def :          natural := 1;
     begin
       loop
	 Put_Line ( " 0:  Quit" );
	 Put_Line ( " 1:  Inject parallel RC soma" );
	 Put_Line ( " 2:  Observe passive cable attenuation" );
	 Put_Line ( " 3:  Neurotransmitter pulse to parallel RC soma" );
	 case Ask_Nat ( "Option ", Def, Min, Max ) is
	   when 0 =>
	     exit;
	   when 1 =>
	     Soma_Parallel_RC;
	   when 2 =>
	     declare
	       Rm       : float := Rm_Def;
	       Ri       : float := Ri_Def;
	       Diameter : float := Diameter_Def;
	     begin
	       Rm := Ask ( "Rm ", Rm );
	       Ri := Ask ( "Ri ", Ri );
	       Diameter := Ask ( "Diameter ", Diameter );
	       Cable_Passive_Steady_State ( Rm, Ri, Diameter );
	     end;
	   when 3 =>
	     declare
	       Vm : float;
	       File : File_Type;
	     begin
	       Create ( File, Out_File, "Pulse.Dat" );
	       Put_Line ( "Output to Pulse.Dat." );
	       for Freq in 1..200 loop
		 Synapse_Soma_Parallel_RC ( Vm, float ( Freq ) );
		 Put ( Freq ); Put ( " " );
		 Put_Line ( Vm );
		 Put ( File, Freq ); Put ( File, " " );
		 Put_Line ( File, Vm );
	       end loop;
	       Close ( File );
	     end;
	   when others =>
	     null;
	 end case;
	 Def := Def + 1;
	 if Def > Max then
	   Def := Min;
	 end if;
       end loop;
     end Demo;

     procedure Soma_Parallel_RC (
       I_Amp : in     float := 1.0 ) is
     ----------------------------------------------------------------------
       Rm         : constant float := 20_000.0; -- Ohms cm2
       Cm         : constant float := 1.0e-6; -- F/cm2
       Time_Max   : constant float := 4.0;
       dt         : constant float := Time_Max / 1000.0;
       Time       :          float;
       I_Ext      :          float;
       dVm_dt     :          float;
       Vm         :          float;
     begin
       for Current in Current_Type'first..Current_Type'last loop
	 Time  := 0.0;
	 I_Ext := 0.0;
	 Vm    := 0.0;
	 while Time <= Time_Max loop
	   I_Ext := I_Amp * Current_Injection ( Current, Time );
	   dVm_dt := -Vm / ( Rm * Cm ) + I_Ext / Cm;
	   Vm := Vm + dVm_dt * dt;
	   Put ( Time , 2, 6, 4 );  Put ( ' ' );
	   Put ( I_Ext, 2, 6, 4 );  Put ( ' ' );
	   Put ( Vm   , 2, 6, 4 );  New_Line;
	   Time := Time + dt;
	 end loop;
       end loop;
     end Soma_Parallel_RC;

     procedure Synapse_Soma_Parallel_RC (
       V         :    out float;
       Frequency : in     float := 0.0 ) is
     ----------------------------------------------------------------------
     -- This model varies synaptic conductance g(t) as an alpha function.
     ----------------------------------------------------------------------
       E_Leak  : constant float := -70.0e-03; -- Volts
       T_Peak  : constant float :=   2.0e-03; -- seconds
       G_Peak  : constant float :=   0.5e-09; -- Siemens
       R_Leak  : constant float := 200.0e+06; -- Ohms
       C       : constant float := 100.0e-12; -- Farads
       E_Syn   : constant float :=   0.0    ; -- Volts
       dt      : constant float := T_Peak / 1.0e2; -- seconds, time delta
       Time    :          float :=   0.0    ; -- seconds
       Time_Max: constant float :=   0.1    ; -- seconds
       Nt      :          float := G_Peak * 2.718 / T_Peak;
       Nt_Last :          float             ; -- seconds
       Vm      :          float := E_Leak   ; -- Volts
       dVm_dt  :          float             ;
       G       :          float :=   0.0    ; -- Siemens
       dG_dt   :          float :=   0.0    ;
       dG_dt_2 :          float             ;
       I_Syn   :          float             ; -- Amperes
     begin
       if Frequency /= 0.0 then
	 Nt_Last := 1.0 / Frequency;
       end if;
       loop
	 if Frequency /= 0.0 then
	   if Nt_Last >= 1.0 / Frequency then
	     Nt_Last := 0.0;
	     Nt := G_Peak * 2.718 / T_Peak;
	   else
	     Nt := 0.0;
	   end if;
	 end if;
	 Put ( Time  ); Put ( "  " );
	 Put ( Nt    ); Put ( "  " );
	 Put ( G     ); Put ( "  " );
	 Put ( I_Syn ); Put ( "  " );
	 Put ( Vm    ); New_Line;
	 exit when Time >= Time_Max;
	 dG_dt_2 := -2.0 * dG_dt / T_Peak - G / ( T_Peak ** 2 ) + Nt / dt;
	 I_Syn := ( Vm - E_Syn ) * G;
	 dVm_dt := -I_Syn / C - ( Vm - E_Leak ) / ( R_Leak * C );
	 dG_dt := dG_dt + dt * dG_dt_2;
	 G := G + dt * dG_dt;
	 Vm := Vm + dt * dVm_dt;
	 Time := Time + dt;
	 if Frequency /= 0.0 then
	   Nt_Last := Nt_Last + dt;
	 else
	   Nt := 0.0;
	 end if;
       end loop;
       V := Vm;
     end Synapse_Soma_Parallel_RC;

     ----------------------------------------------------------------------
     ----------------------------------------------------------------------
     end CNS221;