Program jfactor
      av=6.e23
      pi=3.14159
      bconst=1.38e-16
      temp=295.
      open(11,file='results.txt')

      write(*,*) ' input tors. rigidity (in 10e-19 erg*cm, like 2.4)'
      read (*,*) c
      write(*,*) ' input persistence length (in nm, like 48.0)'
      read (*,*) a
      write(*,*) ' input DNA helical repeat (like 10.51)'
      read (*,*) gamma
      write(*,*) ' input the starting DNA length (in bp, like 191)'
      read (*,*) nbp0

      c=c*1.e-19
      write(*,12) c,a, gamma
      write(*,*)
   12 format(' C(erg*cm) = ',e9.2,'   a(nm) = ',f4.1,'  gamma=', f5.2)
      b=a*2.e-8
      do 1 nbp=nbp0,nbp0+15
      rl=nbp*3.4e-9
      ex=exp(-b/rl*pi**2+0.514*rl/b)
      term1=4.*(pi*b)**3/(av*rl**6)
      rj=term1*ex

c     twist component

           tw0=nbp/gamma
           ntw0=tw0
           dnkuhn=rl/b

           fn=0.752*dnkuhn*exp(-6.242/dnkuhn)
           varwr=(0.00385*dnkuhn+0.113*fn)*dnkuhn/(1.+fn)

           vartw=bconst*temp*(10.*rl)/(4.*pi**2*c)
           varlk=varwr+vartw
           dnorm=sqrt(1./(2.*pi*varlk))

           sum=0.

           do 3 lk=ntw0-3,ntw0+3
             sum=sum+dnorm*exp(-(lk-tw0)**2/(2.*varlk))
    3      continue

      rj=rj*sum


      write(11,14) nbp, rj*1.e9
   14 format(1x,i3,3x,f8.5)
      write(*,5) nbp, rj*1.e9

    1 continue
    5 format('  nbp =',i3, 5x,'   j(in nM) = ',f8.5)
      pause
      end