subroutine pmt5r
      print *, 'pmt5r initialized'
      end

      real function phasfit(f)
      vector par(7)
      real logf
      logf=log(f)/log(10.)
      phasfit=par(1)
     +       +par(2)*exp(-0.5*((logf-par(3))/par(4))**2)
     +       +par(5)*exp(-0.5*((logf-par(6))/par(7))**2)
      end

      complex function Zcomplex(i)
      integer i
      real R1,R2,C
      parameter (R1=5.1e6,R2=0.836e6,C=70e-12)
      complex i_
      parameter (i_=(0.,1.))
      vector y1,y2
      vector f1,f2
      vector ph1,ph2
      complex Z1,Z2
      real omega
      omega=2*3.1415926*f1(i)
      Z1=complex(y1(i)*cos(ph1(i)),y1(i)*sin(ph1(i)))
      Z2=complex(y2(i)*cos(ph2(i)),y2(i)*sin(ph2(i)))
      Zcomplex=(Z1*R1-Z2*R2)/(Z2-Z1-i_*omega*C*(Z1*R1-Z2*R2))
      end

      complex function gain(i)
      integer i,n
      real R1,R2,C
      parameter (R1=5.1e6,R2=0.836e6,C=70e-12)
      complex i_
      parameter (i_=(0.,1.))
      vector y1,y2
      vector f1,f2
      vector ph1,ph2
      complex g1,g2
      complex Zcomplex
      external Zcomplex
      complex Z,Z1,Z2
      real omega
      omega=2*3.1415926*f1(i)
      Z1=complex(y1(i)*cos(ph1(i)),y1(i)*sin(ph1(i)))
      Z2=complex(y2(i)*cos(ph2(i)),y2(i)*sin(ph2(i)))
      Z=Zcomplex(i)
      g1=Z1*(Z+R1*(1+i_*omega*C*Z))
      g2=Z2*(Z+R2*(1+i_*omega*C*Z))
      gain=g1
      end

      real function Zreal(i)
      integer i
      complex Zcomplex
      external Zcomplex
      Zreal=real(Zcomplex(i))
      end

      real function Zimag(i)
      integer i
      complex Zcomplex
      external Zcomplex
      Zimag=aimag(Zcomplex(i))
      end

      real function greal(i)
      integer i
      complex gain
      external gain
      greal=real(gain(i))
      end

      real function gimag(i)
      integer i
      complex gain
      external gain
      gimag=aimag(gain(i))
      end