!*****************
! File Description:  
!   
! Creator: Hoang-Trong Minh Tuan @ George Mason University, 2010
! 
! Rev: 0.1 
! Update: 
!    0.2: (Nov - 8, 2010 - Tuan) rename the file from data_spark.f95 to data_ecc.f95
!    0.3: rename current_state_sfu --> sfu_state 
!    0.4: move gating variables to here (from tscell_gpu_utility, ecc_gpu_utility, ecc_utility)
!    1.0: rename data to reflect using the new algorithm (Mar 3rd, 2012 - Tuan)
!         sfu_state --> CRU_states (1.0)
!*******************/

MODULE data_ecc
  USE precisioninfo
  USE data_global
  USE math_utility
#include "macro.h"

#if _TEST_ORPHAN_CRU_ == _YES
  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE:: aKm2_orphan, aKp2_orphan
  INTEGER, DIMENSION(:), ALLOCATABLE :: U_Ro_orphan
  INTEGER :: irow_R_orphan
#endif

  !+ compact form data
  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: lcc_compKv_act, lcc_compKv_inact, &
       lcc_compKcal_inact, lcc_compKfree,ryr_compKfree_mode1, ryr_compKcal_act_mode1, &
       ryr_compKfree_mode2, ryr_compKcal_act_mode2, ryr_compKfree, ryr_compKRos
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: lcc_indxKv_act, lcc_indxKv_inact, &
       lcc_indxKcal_inact, lcc_indxKfree, ryr_indxKfree_mode1, ryr_indxKcal_act_mode1,&
       ryr_indxKfree_mode2, ryr_indxKcal_act_mode2, ryr_indxKfree, ryr_indxKRos
  INTEGER:: lcc_Ncol_Kv_act, lcc_Ncol_Kv_inact, lcc_Ncol_Kcal_inact, lcc_Ncol_Kfree, ryr_Ncol_Kfree_mode1, &
       ryr_Ncol_Kcal_act_mode1, ryr_Ncol_Kfree_mode2, ryr_Ncol_Kcal_act_mode2, ryr_Ncol_Kfree, ryr_Ncol_KRos

  !+ vector tells number of open LCC/RyR in each homogeneous cluster space
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: U_Lo, U_Ro
  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE:: U_RyR_State
  !+ number of rows in each vector above
  INTEGER :: irow_L, irow_R

  !+ RYRgate(NSFU), DHPRgate(NSFU): # open channels at NSFU release sites
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE ::  RYRgate, DHPRgate

  !+ transition matrices for an assembled channel (deterministic approach)
  ! we ignore the number of channels, i.e. assuming the total is 1, and we're interested in the fraction of channel in each state
  REAL(KIND=dp), DIMENSION(:,:), ALLOCATABLE :: aKv1_nj,aKv2_nj,aKm1_nj,aKp1_nj,aKm2_nj,aKp2_nj
  !+ vector that tell which state is the open state
  INTEGER, DIMENSION(:), ALLOCATABLE :: U_Ro_nj, U_Lo_nj
  !+ the length of each vector
  INTEGER :: irow_L_nj, irow_R_nj
  
  REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: pL_nj, pR_nj
  !+ CRU_states(NSFU): keep index of the current site's state at NSFU release site
  ! the index still based on heterogenous cluster, yet the computation will use the index of homogeneous cluster
  ! based on the mapping strategy from heterogeneous cluster index to homogeneous cluster index
  INTEGER, DIMENSION(:), ALLOCATABLE:: CRU_states, LCCcluster_states, RYRcluster_states
  !+ gating variables
  REAL(KIND=dp) :: m_Na, h_Na, j_Na, a_Ktof, &
       a_Ktos, a_Kss, in_Ktof, in_Kss, in_Ktos
  !+ parameters for gating variables
  REAL(KIND=dp) :: a_h, a_j, b_h, b_j, a_m, b_m

  REAL(KIND=dp) :: kva, kvi  ! Vm-dependent transition rates LCC

  INTERFACE state_combine
     MODULE PROCEDURE state_combine_single_test, state_combine_all
  END INTERFACE
  SAVE
CONTAINS

  !!>****f* data_ecc/cleanup_data_ecc
  !! SYNOPSIS  
  !!   Free data
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN
  !!   
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  SUBROUTINE cleanup_data_ecc()
    if (allocated(lcc_compKv_act)) then
       deallocate(lcc_compKv_act)
    endif
    if (allocated(lcc_compKv_inact)) then
       deallocate(lcc_compKv_inact)
    endif
    if (allocated(lcc_compKfree)) then
       deallocate(lcc_compKfree)
    endif
    if (allocated(lcc_compKcal_inact)) then
       deallocate(lcc_compKcal_inact)
    endif
    if (allocated(lcc_indxKv_act)) then
       deallocate(lcc_indxKv_act)
    endif
    if (allocated(lcc_indxKv_inact)) then
       deallocate(lcc_indxKv_inact)
    endif
    if (allocated(lcc_indxKfree)) then
       deallocate(lcc_indxKfree)
    endif
    if (allocated(lcc_indxKcal_inact)) then
       deallocate(lcc_indxKcal_inact)
    endif
    if (allocated(ryr_compKcal_act_mode1)) then
       deallocate(ryr_compKcal_act_mode1)
    endif
    if (allocated(ryr_compKcal_act_mode2)) then
       deallocate(ryr_compKcal_act_mode2)
    endif
    if (allocated(ryr_compKfree_mode1)) then
       deallocate(ryr_compKfree_mode1)
    endif
    if (allocated(ryr_compKfree_mode2)) then
       deallocate(ryr_compKfree_mode2)
    endif
    if (allocated(ryr_compKfree)) then
       deallocate(ryr_compKfree)
    endif
    if (allocated(ryr_compKRos)) then
       deallocate(ryr_compKRos)
    endif
    if (allocated(ryr_indxKcal_act_mode1)) then
       deallocate(ryr_indxKcal_act_mode1)
    endif
    if (allocated(ryr_indxKcal_act_mode2)) then
       deallocate(ryr_indxKcal_act_mode2)
    endif
    if (allocated(ryr_indxKfree_mode1)) then
       deallocate(ryr_indxKfree_mode1)
    endif
    if (allocated(ryr_indxKfree_mode2)) then
       deallocate(ryr_indxKfree_mode2)
    endif
    if (allocated(ryr_indxKfree)) then
       deallocate(ryr_indxKfree)
    endif
    if (allocated(ryr_indxKRos)) then
       deallocate(ryr_indxKRos)
    endif
    if (allocated(U_Ro)) then
       deallocate(U_Ro)
    endif
    if (allocated(U_RyR_State)) then
       deallocate(U_RyR_State)
    endif
    if (allocated(U_Lo)) then
       deallocate(U_Lo)
    endif
    if (allocated(U_Ro_nj)) deallocate(U_Ro_nj)
    if (allocated(U_Lo_nj)) deallocate(U_Lo_nj)
    if (allocated(aKv1_nj)) deallocate(aKv1_nj)
    if (allocated(aKv2_nj)) deallocate(aKv2_nj)
    if (allocated(aKm1_nj)) deallocate(aKm1_nj)
    if (allocated(aKp1_nj)) deallocate(aKp1_nj)
    if (allocated(aKm2_nj)) deallocate(aKm2_nj)
    if (allocated(aKp2_nj)) deallocate(aKp2_nj)
  END SUBROUTINE cleanup_data_ecc
  
  !!>****f* data_ecc/init_CRU_states
  !! SYNOPSIS  
  !!    set the initial state for a Calcium Release Unit (CRU = SFU = stochastic functional unit)
  !!    By default, we expect all channels start with close-state
  !! INPUTS
  !!   lcc_initialstate :: initial state for LCCs (expected to be state index in Markov that map to Close)
  !!   ryr_initialstate :: initial state for RyRs (expected to be state index in Markov that map to Close)
  !!   N_L              :: number of LCCs
  !!   mL               :: number of states in LCC Markov
  !!   N_R              :: number of RyRs
  !!   mR               :: number of states in RyR Markov
  !! OUTPUTS 
  !!   CRU_states       :: state index for all CRUs
  !!   LCCcluster_states:: state index for all LCC clusters
  !!   RYRcluster_states:: state index for all RYR clusters
  !! RETURN
  !!   
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !!_****
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  SUBROUTINE init_CRU_states(lcc_initialstate, ryr_initialstate, N_L, mL, N_R, mR,  &
       CRU_states, LCCcluster_states, RYRcluster_states)
    !! work for EC-coupling where LCC + RyR form a heterogeneous cluster
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: lcc_initialstate, ryr_initialstate, N_L, mL, N_R, mR
    INTEGER, DIMENSION(:), ALLOCATABLE:: CRU_states, LCCcluster_states, RYRcluster_states
    ! local
    INTEGER :: balls, bins, irows, icols, curr_state, state1, state2
    INTEGER, DIMENSION(:,:), ALLOCATABLE :: LC_in_state, RyR_in_state
    INTEGER :: ii, irow_R

    balls = N_L
    bins = mL
    CALL count_ball2bin(balls, bins, LC_in_state, irows, icols)
    !The triggered voltage is applied only when the system reaches the steady-state, i.e.
    !the initial state is 2 (close state), instead of 1 (inactivate state)
    DO ii=1,irows
       IF (LC_in_state(ii,lcc_initialstate) .EQ. N_L) THEN
          state1 = ii
          EXIT
       ENDIF
    END DO
    LCCcluster_states = state1
    
#if _DEBUG_ == 3
    PRINT *, state1
#endif    

    balls = N_R
    bins = mR
    CALL count_ball2bin(balls, bins, RyR_in_state, irows, icols)
    irow_R = irows
    
    DO ii=1,irow_R
       !IF (RyR_in_state(ii,ryr_initialstate) .EQ. N_R) THEN
       IF (RyR_in_state(ii,1) .EQ. 49) THEN
          state2 = ii
          EXIT
       ENDIF
    END DO
    RYRcluster_states = state2

    !PRINT *, 'irow_R', irow_R
    
    !CRU_states = state_combine(state1, state2, irow_R)
    CRU_states = state_combine_all(LCCcluster_states, RYRcluster_states, irow_R)

!!$    PRINT *, CRU_states(1:5)
!!$    PRINT *, LCCcluster_states(1:5)
!!$    PRINT *, RYRcluster_states(1:5)
!!$    STOP
    
#if _DEBUG_ == 3
    PRINT *, CRU_states(1:10)
    stop
#endif

    PRINT *, "++end init_SFU_state()"
  END SUBROUTINE init_CRU_states

  !!>****f* data_ecc/init_data_ECC
  !! SYNOPSIS  
  !!   Initialize variables for gating variables
  !!   and related data to EC-coupling
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  SUBROUTINE init_data_ECC()
    IMPLICIT NONE
    PRINT *, "start init_data_ECC()"
    
    m_na    = 0.0015d0
    h_na    = 0.9849d0
    j_na    = 0.9905d0
    a_ktof  = 0.0021d0
    in_ktof  = 1.d0
    a_ktos  = 2.9871d-04
    in_ktos  = 0.9994d0
    a_kss   = 0.002d0
    in_kss   = 1.d0

    PRINT *, "++end init_data_ECC()"
  END SUBROUTINE init_data_ECC

  !!>****f* data_ecc/init_data_ECC
  !! SYNOPSIS  
  !!   This subroutine break down from the heterogeneous cluster state index
  !!   to homogeneous cluster state index
  !! INPUTS 
  !! OUTPUTS 
  !! RETURN 
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  !attributes(device, host)
  SUBROUTINE state_breakdown(cru_state, irow_R, lcc_state, ryr_state)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: cru_state, irow_R
    INTEGER, INTENT(INOUT) :: lcc_state, ryr_state
    
    lcc_state = (cru_state-1)/irow_R + 1
    ryr_state = cru_state - (lcc_state-1)*irow_R
    !PRINT *, lcc_state, ryr_state, cru_state, irow_R
  END SUBROUTINE state_breakdown
  SUBROUTINE state_breakdown_all(cru_state, irow_R, lcc_state, ryr_state)
    IMPLICIT NONE
    INTEGER, DIMENSION(:), INTENT(IN) :: cru_state
    INTEGER, INTENT(IN) :: irow_R
    INTEGER, DIMENSION(:), INTENT(INOUT) :: lcc_state, ryr_state
    
    lcc_state = (cru_state-1)/irow_R + 1
    ryr_state = cru_state - (lcc_state-1)*irow_R
    !PRINT *, lcc_state, ryr_state, cru_state, irow_R
  END SUBROUTINE state_breakdown_all

  !!>****f* data_ecc/state_combine_single
  !! SYNOPSIS  
  !!   This subroutine find the  heterogeneous cluster state index
  !!   for a single CRU based on component homogeneous cluster state indices
  !! INPUTS
  !!   LCCcluster_state :: state index for LCC cluster
  !!   RYRcluster_state :: state index for RYR cluster
  !!   irow_R           :: number of states in RYR cluster
  !! OUTPUTS 
  !! RETURN
  !!   cru_state        :: state index of CRU (LCCs + RYRs)
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
  !attributes(device, host)
  FUNCTION state_combine_single(LCCcluster_state, RYRcluster_state, irow_R) RESULT(cru_state)
    INTEGER :: cru_state
    INTEGER, INTENT(IN) :: LCCcluster_state, RYRcluster_state
    INTEGER, INTENT(IN) :: irow_R
    
    cru_state = (LCCcluster_state-1) * irow_R + RYRcluster_state

    RETURN
  END FUNCTION state_combine_single

  FUNCTION state_combine_single_test(LCCcluster_state, RYRcluster_state, irow_R) RESULT(cru_state)
    INTEGER :: cru_state
    INTEGER, INTENT(IN) :: LCCcluster_state, RYRcluster_state
    INTEGER, INTENT(IN) :: irow_R
    
    cru_state = (LCCcluster_state-1) * irow_R + RYRcluster_state

    RETURN
  END FUNCTION state_combine_single_test

  !!>****f* data_ecc/state_combine
  !! SYNOPSIS  
  !!   This subroutine find the  heterogeneous cluster state index
  !!   for all CRUs based on component homogeneous cluster state indices
  !!   at each CRU
  !! INPUTS
  !!   LCCcluster_states :: state index for LCC cluster at all CRUs
  !!   RYRcluster_states :: state index for RYR cluster at all CRUs
  !!   irow_R           :: number of states in RYR cluster
  !! OUTPUTS 
  !! RETURN
  !!   cru_state        :: state index at all CRUs (each one is composed of LCCs + RYRs)
  !! SIDE EFFECTS 
  !! HISTORY 
  !! BUGS 
  !! TODO 
  !! SEE ALSO 
  !!_****
!!$  FUNCTION state_combine_all(LCCcluster_states, RYRcluster_states, irow_R) 
!!$  !subroutine state_combine_all(LCCcluster_states, RYRcluster_states, irow_R, cru_states) 
!!$    INTEGER, DIMENSION(:) :: state_combine_all
!!$    INTEGER, DIMENSION(:) :: LCCcluster_states, RYRcluster_states
!!$    INTEGER, INTENT(IN) :: irow_R
!!$    
!!$    state_combine_all = (LCCcluster_states-1) * irow_R + RYRcluster_states
!!$    
!!$    RETURN
!!$  END FUNCTION state_combine_all
  FUNCTION state_combine_all(LCCcluster_states, RYRcluster_states, irow_R) RESULT(cru_states)
  !subroutine state_combine_all(LCCcluster_states, RYRcluster_states, irow_R, cru_states) 
    INTEGER, DIMENSION(:) :: cru_states
    INTEGER, DIMENSION(:) :: LCCcluster_states, RYRcluster_states
    INTEGER, INTENT(IN) :: irow_R
    
    cru_states = (LCCcluster_states-1) * irow_R + RYRcluster_states
    
    RETURN
  END FUNCTION state_combine_all

  FUNCTION find_v_ryr(i_1ryr, Vol_ref, NSFU, N_R) RESULT(v_ryr)
    IMPLICIT NONE
    REAL(KIND=dp) :: v_ryr
    !i_1ryr = single channel current
    !Vol_ref = reference volume
    REAL(KIND=dp), INTENT(IN) :: i_1ryr, Vol_ref
    INTEGER, INTENT(IN) :: NSFU !number release sites
    INTEGER, INTENT(IN) :: N_R  !number RyRs/release site

    !local
    REAL(KIND=dp) :: Ca_jsr, Ca_ds

    Ca_jsr = 1.0d3 !uM
    Ca_ds = 0.1d0 !uM
    
    v_ryr = i_1ryr*1d6 / (zCa*zF*Vol_ref*(Ca_jsr-Ca_ds))  ![1/sec]
    v_ryr = v_ryr * N_R * NSFU
    RETURN 
  END FUNCTION find_v_ryr
END MODULE data_ecc
