c     program main
c     call preamp(5.6e3,2.0e3,1.0e3,90.,5.0,2e-15)
c     print *, Zin(1e-4)
c     end

      subroutine preamp(R3val,R4val,R5val,betaval,Vccval,I0val)
      real R3val,R4val,R5val,betaval,Vccval,I0val
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      real qElectron,kBoltzmann,T
      R3=R3val
      R4=R4val
      R5=R5val
      beta=betaval
      Vcc=Vccval
      I0=I0val
      qElectron=1.602e-19
      kBoltzmann=1.380e-23
      T=293
      Vt=kBoltzmann*T/qElectron
      end

      subroutine modelsolver(Iin)
      real Iin
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision Vclast,Vblast
      real convergespeed
      data convergespeed/0.1/
      save convergespeed
   1  Vb=0.6
      Vc=Vcc
      Vblast=Vb
      Vclast=Vc
      do i=1,9999
c       Vc=-Iin*R3+Vb*(1+R3/R4)+I0*R3*exp(Vb/Vt)/beta
        Vc=R5*R3/(R5+R3)*(Vcc/R5+Vb/R3-I0*exp(Vb/Vt)*(1-1/beta))
        Vb=R3*R4/(beta*(R3+R4)-R3)*((beta-1)*Iin+(Vc-V0)/R5+beta*Vc/R3)
        if (abs(Vb-Vblast).lt.1e-10 .and.
     +      abs(Vc-Vclast).lt.1e-10) then
          goto 99
        else
          Vb=(1-convergespeed)*Vblast+convergespeed*Vb
          Vblast=Vb
          Vclast=Vc
        endif
      enddo
      if (convergespeed.gt.1e-10) then
        convergespeed=convergespeed/2
        goto 1
      endif
      print *, 'modelsolver failed to converge, giving up!'
      return
   99 continue
c     print *, 'modelsolver converged in',i,' steps'
      end

      real function Vout(Iin)
      real Iin
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      call modelsolver(Iin)
      Vout=Vc
      end

      real function Vin(Iin)
      real Iin
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      call modelsolver(Iin)
      Vin=Vb
      end

      real function Zin(Iin)
      real Iin
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision Vb0
      call modelsolver(0.)
      Vb0=Vb
      call modelsolver(Iin)
      Zin=(Vb-Vb0)/Iin
      end

      real function Zout(Iin)
      real Iin
      common /preamp_model/R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision R3,R4,R5,beta,Vcc,Vb,Vc,Vt,I0
      double precision Vc0
      call modelsolver(0.)
      Vc0=Vc
      call modelsolver(Iin)
      Zout=(Vc-Vc0)/Iin
      end