c     hh.f
c     hodgkin-huxley program for a space clamped axon at
c     6.3 centigrade

      implicit none
      real*8 n,m,h,an,am,ah,bn,bm,bh,dndt,dmdt,dhdt,dvdt
      real*8 gk,gna,jk,jna,jl,jmemb,v,t,jstim,tp
      real*8 deltat,tpstep,vrest,cmemb,tmax,vna,vk
      integer i,j
      character*80 filea,fileb

      write(*,*) 'Enter the name of the conductance output file'
      read(*,'(a)') filea
      write(*,*) 'Enter the name of the current output file'
      read(*,'(a)') fileb
      open(10,file=filea,status='new')
      open(11,file=fileb,status='new')

c     define some constants

      deltat= 1.0e-6 ! integration step
      tpstep= 1.0e-4 ! write out a line this often
      vrest= -65.0e-3 ! resting potential
      cmemb = 1.0e-6 ! membrane capacitance
      tmax = 5.0e-3 ! time to quit
      vna = 50.0e-3 ! nernst potential for sodium
      vk = -77.0e-3 ! nernst potential for potassium

c******************************************************************

c     main program for calculating the response of a space
c     clamped axon to a stimulus

      t=0.d0
      v=vrest
      tp=0.d0

      call calcinitvals(n,m,h,an,am,ah,bn,bm,bh,v,vrest)

      do while(t.lt.tmax)

       call calccurr(gk,gna,jk,jna,jl,jmemb,v,vrest,n,m,h,vk,vna)
       if(t.ge.tp) then   !write out the results only once in awhile
           write(10,100) t,v,gna,gk
           write(11,111) t,v,jmemb,jna,jk,jl
 100       format(4(1x,e13.5))
 111       format(6(e13.5))
           tp =tp+tpstep
       else
           continue
       endif

       if((t.ge.(5.0e-4)).and.(t.lt.(6.0e-4)))then !apply stimulus
           jstim=1.0e-4
       else
           jstim=0.0
       endif
       

       dvdt=(-jmemb+jstim)/cmemb
       v= v+(dvdt*deltat)

       call calcab(an,am,ah,bn,bm,bh,v,vrest)

       dndt=an*(1-n)-bn*n
       dmdt=am*(1-m)-bm*m
       dhdt=ah*(1-h)-bh*h
       n=n+dndt*deltat
       m=m+dmdt*deltat
       h=h+dhdt*deltat
       t=t+deltat
      enddo
      
        stop
        end

c*******************************************************************
      subroutine calcab(an,am,ah,bn,bm,bh,v,vrest)
c     calculate the alpha and betas for the hh equations
c     a's and b's are in volts and seconds

      implicit none
      real*8 an,am,ah,bn,bm,bh,vrest,v
      
      an = 10.*(-1000.*(v-vrest)+10.)/
     1     (exp((-1000.*(v-vrest)+10.)/10.)-1.)
      am = 100.*(-1000.*(v-vrest)+25.)/
     1     (exp((-1000.*(v-vrest)+25.)/10.)-1.)
      ah = 70.*exp(-1000.*(v-vrest)/20.)
      bn = 125.*exp(-1000.*(v-vrest)/80.)
      bm = 4000.*exp(-1000.*(v-vrest)/18.)
      bh = 1000./(exp((-1000.*(v-vrest)+30.)/10.)+1.)

      return
      end

c*******************************************************************
c
      subroutine calcinitvals(n,m,h,an,am,ah,bn,bm,bh,v,vrest)

c     calculate initial values of n,m,h

      implicit none
      real*8 n,m,h,an,am,ah,bn,bm,bh,v,vrest

      call calcab(an,am,ah,bn,bm,bh,v,vrest)
      n=an/(an+bn)
      m=am/(am+bm)
      h=ah/(ah+bh)

      return
      end

c*******************************************************************
c
      subroutine calccurr(gk,gna,jk,jna,jl,jmemb,v,vrest,
     1                     n,m,h,vk,vna)

c     calculate conductances in S/cm**2 and current densities
c     in A/cm**2
      implicit none
      real*8 gk,gna,jk,jna,jl,jmemb,vna,vrest,vk,n,m,v,h

      gk=36.0e-3*(n**4)
      gna = 120.0e-3*h*(m**3)
      jk = gk*(v-vk)
      jna = gna*(v-vna)
      jl = 3.0e-4*(v-(vrest+10.6e-3))
      jmemb = jk+jna+jl

      return
      end