	subroutine fractal_ds

*	See Mildner & Hall, J Phys D Appl Phys, 19, (1986), 1535-1545
*	See Triolo et al, J Appl Cryst, 33, (2000), 863-866

	include 'frills_sources:function.inc'

	real*8		prefac,r0,ds,eta,bkgd
	real*8		a,b,c,d,e,f,g,h,o,q,pq

	external gammafunc

	prefac=p(1)
	r0=p(2)
	ds=p(3)
	eta=p(4)
	bkgd=p(5)

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

	if (ds.gt.6) then
	   write(6,*)' '
	   write(6,*)'The surface fractal dimension must be <=6!'
	   write(6,*)'Will set to 2.'
	   ds=2.0
	end if

	a=5.0-ds
	call gammafunc(a,g)

10	format(' Ds=',f12.5,'     Gamma(5-Ds)=',e12.5)
	write(6,10)ds,g

	b=(ds-5.0)/2.0

*	SANDRA will crash with %MTH-F-UNDEXP if eta goes negative...
	c=abs(eta)**a

	do j=1,nv
	   d=x(j)*eta
	   e=(1.0+(d*d))**b
	   f=dsin(-1*a*datan(d))
           h=x(j)*r0
           o=dsin(h)
           q=dcos(h)
*	Use a spherical form factor for P(q)
	   if (r0.eq.0.0) then
	      pq=1.0
	   else
              pq=(3.*(o-h*q)/(h*h*h))*(3.*(o-h*q)/(h*h*h))
	   endif
*	Use S(q)=f(ds,eta) rather than S(q)=1+f(ds,eta)
	   ycal(j)=((prefac*pq)*((g*c*e*f)/x(j)))+bkgd
	end do

	return
	end
