real function pshape(x)
      real x
      vector pspar(3)
      real Erf
      external Erf
      real root2pi
      root2pi=sqrt(2*3.1415926)
      pshape = 1-Erf((x+pspar(1)*pspar(2)**2)/(sqrt(2.)*pspar(2)))
      if (pshape.gt.0) then
        pshape=pshape*pspar(1)/2
     +        *exp(x*pspar(1)+0.5*(pspar(1)*pspar(2))**2)
        pshape=pspar(3)/(root2pi*pspar(2))*exp(-0.5*(x/pspar(2))**2)
     +        +(1-pspar(3))*pshape
      endif
      end

      real function fitpeaks(x)
      real x
      vector pspar(3)
      real par
      common /pawpar/par(12)
      real peak
      real peaks
      real lambda
      real pshape
      external pshape
      real q
      integer p,maxp
      parameter (maxp=15)
      integer s,maxs
      parameter (maxs=15)
      do p=1,3
        pspar(p)=par(p+6)
      enddo
      expfact=exp(-par(4))
      q=(x-par(2))/par(3)
      fitpeaks=0
      do p=0,maxp
        peaks=0
        lambda=par(5)*p**1.9
        expfacts=exp(-lambda)
        do s=0,maxs
          pspar(2)=sqrt(par(8)**2+(p+s)*par(6)**2)
          peak=expfact*(par(4)**p)/gamma(p+1.)
     +        *expfacts*(lambda**s)/gamma(s+1.)
     +        *pshape(q-p-s)
          if (peak.ge.peaks*1e-5) then
            peaks=peaks+peak
          else
            goto 5 ! exit
          endif
        enddo
    5   continue
        if (peaks.ge.fitpeaks*1e-5) then
          fitpeaks=fitpeaks+peaks
        else
          goto 6 ! exit
        endif
      enddo
    6 continue
      fitpeaks=fitpeaks*par(1)
      end

      real function fitpeaksf(x)
      real x
      vector par(9)
      common /pawpar/ppar(12)
      real ppar
      real fitpeaks
      external fitpeaks
      integer p
      do p=1,9
        ppar(p)=par(p)
      enddo
      fitpeaksf=fitpeaks(x)
      end