      program gouychapman
      implicit doubleprecision (a-h,o-z)
      
      pi=3.1415927
      twopi=2*pi
      sq2=dsqrt(2.d0)
      sq3=dsqrt(3.d0)
      avog=6.02214d23
      tosi=1.0/avog/1000.
      dsi=1.d-10
      eps0=8.8542d-12
      boltz=1.3806d-23
      echg=1.6021917d-19

      tsi=298.15
      eps=78.46
      c1=1.0
      c2=1.0
      z1=1.0
      z2=-1.0
      d=3.0
      beta=1./(boltz*tsi)
      
      vred=-2.212
      fromv=boltz*tsi/echg
      v=vred*fromv
      
      write(*,*)"                         T = ",tsi
      write(*,*)"                   epsilon = ",eps
      write(*,*)"           c (cation) in M = ",c1
      write(*,*)"             c (anion) i M = ",c2
      write(*,*)"         valence of cation = ",z1
      write(*,*)"          valence of anion = ",z2
      write(*,*)"             d in Angstrom = ",d
      write(*,*)"                V in volts = ",v
      write(*,*)"                 V in kT/e = ",vred
      write(*,*)


      d=d*dsi
      
      dens1=c1/tosi
      dens2=c2/tosi
      write(*,*)"          Density in 1/m^3 = ",dens1
      write(*,*)"          Density in 1/d^3 = ",dens1*d**3
      write(*,*)
      
      q1=z1*echg
      q2=z2*echg

c      v=v/(beta*echg)
      dkappa=q1**2*dens1+q2**2*dens2
      dkappa=dkappa/(eps0*eps*boltz*tsi)
      dkappa=dsqrt(dkappa)
      write(*,*)"              kappa in 1/m = ",dkappa 
      write(*,*)"              kappa in 1/A = ",dkappa*dsi
      write(*,*)"              kappa in 1/d = ",dkappa*d
      write(*,*)
      
      cap=eps0*eps*dkappa/(1+d*dkappa/2.0)
      write(*,*)"               C in C/Vm^2 = ",cap
      write(*,*)
      
      sigma=v*cap
      write(*,*)"            sigma in C/m^2 = ",sigma
      write(*,*)"            sigma in e/d^2 = ",sigma*d*d/echg
      write(*,*)
      
      vd=v-d*sigma/(2*eps0*eps)
      write(*,*)"               Vd in volts = ",vd
      write(*,*)"                Vd in kT/e = ",vd/fromv

      write(*,*)
      open(2,file="cation.dat")
      open(12,file="anion.dat")
      open(22,file="charge.dat")
      do i=1,400
        x=d/2.0+(i-1)*0.05*dsi
        gcat=c1*(1-beta*q1*vd*dexp(-dkappa*(x-d/2.0)))
        gan=c2*(1-beta*q2*vd*dexp(-dkappa*(x-d/2.0)))
        write(2,*)x/d,gcat
        write(12,*)x/d,gan
        write(22,*)x/d,z1*gcat+z2*gan
      enddo 
      close(2)
      close(12)
      close(22)

      open(2,file="potential-reduced.dat")
      open(12,file="potential-volt.dat")
      write(2,*)0.0,v/fromv
      do i=1,400
        x=d/2.0+(i-1)*0.05*dsi
        vx=vd*dexp(-dkappa*(x-d/2.0))
        write(2,*)x/d,vx/fromv
        write(12,*)x/d,vx
      enddo
      close(2)
      close(12)


        

      stop
      end

        
