Ticket #84: CHEN.FOR

        subroutine chen

*       See Chen, Rouch & Tartaglia, Croat. Chem. Acta, 65(2), (1992), 353-366

        include 'frills_sources:function.inc'

        integer*4       lolim
        real*8          scale,r0,s,tau,df,bkgd
        real*8          pi
        real*8          rg,h,eta,twominustau,threeminustau,gamma2mt,gamma3mt
        real*8          i0,q,qeta,u,incogau,f1,f2,g1,g2,incogaxh,vol

        external incog,gammafunc

        scale=p(1)
        r0=p(2)
        s=p(3)
        tau=p(4)
        df=p(5)
        bkgd=p(6)

        pi=3.141592654

        lolim=1

*       For some reason the P O F or P O C extrapolation is introducing
*       some negative x-values at the start of the extrapolated data.  These
*       cause a %MTH-F-UNDEXP, undefined exponentiation error in the
*       calculation of u1.
        do j=1,nv
           if (x(j).lt.0.) then
              lolim=j
           end if
        end do

10      format(' Low Q points being ignored : ',i3)
        write(6,10)lolim

        if (df.le.0.) then
           write(6,*)' '
           write(6,*)'The fractal dimension must be >0!'
           write(6,*)'Will set to 0.1'
           df=0.1
        end if

        if (df.gt.6.) then
           write(6,*)' '
           write(6,*)'The fractal dimension must be <=6!'
           write(6,*)'Will set to 6.'
           df=6.0
        end if

        if (s.lt.1.) then
           write(6,*)' '
           write(6,*)'The aggregation number must be >=1!'
           write(6,*)'Will set to 1.'
           s=1.0
        end if

        if (tau.lt.1.) then
           write(6,*)' '
           write(6,*)'The polydispersity of the cluster must be >=1!'
           write(6,*)'A smaller number means GREATER polydispersity.'
           write(6,*)'Will set to 1.'
           tau=1.0
        end if

        if (tau.eq.2.) then
           write(6,*)'Adjusting polydispersity of the cluster to'
           write(6,*)'avoid a high performance arithmetic trap in the'
           write(6,*)'calculation of the incomplete Gamma function.'
           write(6,*)'Was 2.00, will set to 2.01'
           tau=2.01
        end if

        if (tau.eq.3.) then
           write(6,*)'Adjusting polydispersity of the cluster to'
           write(6,*)'avoid a high performance arithmetic trap in the'
           write(6,*)'calculation of the incomplete Gamma function.'
           write(6,*)'Was 3.00, will set to 3.01'
           tau=3.01
        end if

        if (r0.le.0.) then
           write(6,*)'The primary particle radius must be >=0!'
           write(6,*)'Will set to 5.00'
           r0=5.
        end if

*       Convert primary particle radius to equivalent Rg [Ang]
        rg=sqrt((3./5.)*r0*r0)

*       Volume of primary particle (in cm^3 !!!)
        vol=(4./3.)*pi*(r0*r0*r0)*1.e-24

        h=sqrt((df*(df+1.))/6.)

*       Cluster correlation length
        eta=h*rg*(s**(1./df))

20      format(' Equivalent Cluster Rg [Ang]: ',e12.5)
        write(6,20)rg
30      format(' Cluster corr. length [Ang] : ',e12.5)
        write(6,30)eta

        twominustau=2.-tau
        call gammafunc(twominustau,gamma2mt)

        threeminustau=3.-tau
        call gammafunc(threeminustau,gamma3mt)

        i0=s*gamma3mt/gamma2mt

40      format(' s * Gam(3-tau)/Gam(2-tau)  : ',e12.5)
        write(6,40)i0

*       do j=1,nv
        do j=lolim+1,nv
           q=x(j)
           qeta=q*eta
           u=(h*h*(1.+(qeta*qeta))/(qeta*qeta))**(df/2.)
           call incog(threeminustau,u,incogau)
           f1=gamma3mt-incogau
           f2=(1.+(qeta*qeta))**((-1.*df*threeminustau)/2.)
           g2=(qeta/h)**(-1.*df)
           call incog(twominustau,g2,incogaxh)
           g1=sin((df-1.)*pi/2.)*incogaxh/(df-1.)

           ycal(j)=(scale*vol*s/gamma2mt)*((f1*f2)+(g1*g2))+bkgd
        end do

        return
        end