c   CV simulation - reversible Redox in solution - macro electrode
c   negative scan, D_A = D_B
c........1.........2.........3.........4.........5.........6.........7..
c23456789|123456789|123456789|123456789|123456789|123456789|123456789|12

      implicit real*8 (a-h,l,o-z)	! except "l" for lambda
      parameter(npa=10000)
      dimension gamma_mod(npa), CX(npa)

c ---------------------------------------------------------------
      pi = acos(-1.0d0)
      write (6,*) 'pi =', pi
c ---------!! input section !! change these parameters ---------------

      Ei = 0.5d0		! initial potential / V
      Ev = -0.5d0		! vertex potential / V
      vCV = 100.0d0		! scan rate / mV s-1

      nx = 1000			! number of mesh in x direction (integer < npa-1)
      Delta_E = 1.0d0		! potential step / mV

      E0 = 0.0d0		! formal potential / V
      c_ox_inf = 1.0d-3	! initial Ox concentration / mol dm-3
      D_ox = 1.0d-9		! Diffusion coefficient of Ox / m2 s-1
      eps = 5.0d-3		! electrode radius / m
      Temp = 25.0d0+273.15d0	! temperature / K

      open (2,file="output.csv")

      if(nx+1 > npa) then
       stop 1042
      end if

c --- unit conversion
      vCV = vCV*1.0d-3			! mV s-1 -> V s-1
      Delta_E = Delta_E*1.0d-3	! mV -> V
      c_ox_inf = c_ox_inf*1.0d3	! mol dm-3 -> mol m-3

c--- physical constant & parameter
      Fara = 96485.0d0		! Faraday constant / C mol-1
      Rgas = 8.3145d0		! gas constant / J mol-1 K-1

      F_RT = Fara/(Rgas*Temp)

c - Dimensionless --------------------------------------------
      d_ls_ox = D_ox/D_ox		! dimensionless diffusion coefficient of O

      theta_i = F_RT*(Ei-E0)		! dimensionless initial potential
      
      theta_v = F_RT*(Ev-E0)		! dimensionless vertex potential
      Delta_theta = F_RT*Delta_E	! dimensionless potential step
      sigma = eps**2.0d0/D_ox*F_RT*vCV	! dimensionless scan rate

      write (6,*) 'd_ls_ox =', d_ls_ox
      write (6,*) 'theta_i, theta_v =', theta_i, theta_v
      write (6,*) 'Delta_theta, sigma =', Delta_theta, sigma

c - Calculate Delta_T, m_T, Delta_X --------------------------
      Delta_T = Delta_theta/sigma
      T_max = 2.0d0*(abs(theta_i - theta_v))/sigma
      mT = int(T_max/Delta_T)				! whole step number

      T_switch = abs(theta_i - theta_v)/sigma
      mT_v = int(T_switch/Delta_T)			! step number at vertex potential

      write (6,*) 'Delta_T, T_max =', Delta_T, T_max
      write (6,*) 'mT, mT_v =', mT, mT_v

      X_max = 6.0d0*sqrt(T_max)
      Delta_X = X_max/real(nX,kind(0d0))

      write (6,*) 'X_max, Delta_X =', X_max, Delta_X

c - Calculate Thomas coefficients ----------------------------
      lambda = d_ls_ox*Delta_T/(Delta_X)**2.0d0
      alpha = -lambda
      beta = 1.0d0 + 2.0d0*lambda
      gamma = -lambda

c --- modified gamma
      gamma_mod = 0.0d0 ! *1

      gamma_mod(1) = 0.0d0 ! *2
      do i = 2, nX
       gamma_mod(i) = gamma/(beta-gamma_mod(i-1)*alpha)
      end do

c ------------------------------------------------------------
c --- begin simulation ---------------------------------------
c ------------------------------------------------------------

      write (2,*) 'potential /V', ',', 'current density / mA cm-2' ! *3

      theta = theta_i ! *4
      CX = 1.0d0

c --- CV start ---------------------------------------------
      do k=1,mT+1 ! *5

c      --- potential *6
       if (k == 1) then
        theta = theta_i
       else if (2 <= k .and. k <= mT_v+1) then
        theta = theta - Delta_theta
       else
        theta = theta + Delta_theta
       end if

c      --- modified delta *7
       CX(1) = 1.0d0/(1.0d0 + exp(-theta))
       do i=2, nX
        CX(i) = (CX(i) - CX(i-1)*alpha)
     &               /(beta - gamma_mod(i-1)*alpha)
       end do

c      --- back substitution *8
       CX(nX+1) = 1.0d0
       do i=nX,1, -1
        CX(i) = CX(i) - gamma_mod(i)*CX(i+1)
       end do

c ------ output current -----------------------------------------
       cur_ls = -(-CX(3) + 4.0d0*CX(2) -3.0d0*CX(1)) ! *9
     &                   /(2.0d0*Delta_X)

       potential = theta/F_RT+E0 ! *10
       cur = cur_ls*(Fara*c_ox_inf*D_ox/eps)*1.0d3/1.0d4	! mA cm-2

       write(2,*) potential, ',', cur

c --- CV finish -------------------------------------------------
      end do ! *5

      close(2)

      stop
      end