	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
