	subroutine fractal_qton

*	See Schmidt, J Appl Cryst, 24, (1991), 414-435, Eqn (19)
*	See Hurd, Schaefer & Martin, 35, (1987), 2361-2364

	include 'frills_sources:function.inc'

	real*8		prefac,dm,ds,big_rg_small_rg,bkgd
	real*8		a,b,c,d,e,f,g,h,i

	prefac=p(1)
	dm=p(2)
	big_rg=p(3)
	ds=p(4)
	small_rg=p(5)
	bkgd=p(6)

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

        if (dm.gt.6) then
           write(6,*)'The mass fractal dimension must be <=6!'
           write(6,*)'Will set to 3.'
           dm=3.0
        end if

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

	a=dm/2.0
	b=(big_rg*big_rg)/(3.0*a)

10	format(' Will set to (6-Dm)=',f12.5)
*	IF C GOES NEGATIVE, SANDRA WILL CRASH WITH UNDEFINED EXPONENTIATION
*	SO (Ds+Dm)<=6
*	c=((ds-6.0)/-2.0)-a
	if ((ds.gt.(6.0-dm)).and.(small_rg.gt.0.0)) then
           ds=6.0-dm
           write(6,*)'The surface fractal dimension must be <=(6-Dm)!'
           write(6,10)ds
	endif
*	c=(ds-6.0+dm)/-2.0
	c=(6.0-ds-dm)/2.0

*	IF C=0 THEN SANDRA WILL CRASH WITH A FLOATING DIVIDE BY ZERO
	if (c.eq.0.) then
	   d=1.0E+37
	else
	   d=(small_rg*small_rg)/(3.0*c)
	endif

	do j=1,nv
	   e=x(j)*x(j)*b
	   f=x(j)*x(j)*d
	   g=(1.+e)**a
	   h=(1.+f)**c
	   i=g*h
	   ycal(j)=(prefac/i)+bkgd
	end do

	return
	end
