Ticket #84: DS.FOR

File DS.FOR, 1.3 KB (added by jhjcho, 12 years ago)
Line 
1        subroutine fractal_ds
2
3*       See Mildner & Hall, J Phys D Appl Phys, 19, (1986), 1535-1545
4*       See Triolo et al, J Appl Cryst, 33, (2000), 863-866
5
6        include 'frills_sources:function.inc'
7
8        real*8          prefac,r0,ds,eta,bkgd
9        real*8          a,b,c,d,e,f,g,h,o,q,pq
10
11        external gammafunc
12
13        prefac=p(1)
14        r0=p(2)
15        ds=p(3)
16        eta=p(4)
17        bkgd=p(5)
18
19        if (ds.le.0) then
20           write(6,*)' '
21           write(6,*)'The surface fractal dimension must be >0!'
22           write(6,*)'Will set to 2.'
23           ds=2.0
24        end if
25
26        if (ds.gt.6) then
27           write(6,*)' '
28           write(6,*)'The surface fractal dimension must be <=6!'
29           write(6,*)'Will set to 2.'
30           ds=2.0
31        end if
32
33        a=5.0-ds
34        call gammafunc(a,g)
35
3610      format(' Ds=',f12.5,'     Gamma(5-Ds)=',e12.5)
37        write(6,10)ds,g
38
39        b=(ds-5.0)/2.0
40
41*       SANDRA will crash with %MTH-F-UNDEXP if eta goes negative...
42        c=abs(eta)**a
43
44        do j=1,nv
45           d=x(j)*eta
46           e=(1.0+(d*d))**b
47           f=dsin(-1*a*datan(d))
48           h=x(j)*r0
49           o=dsin(h)
50           q=dcos(h)
51*       Use a spherical form factor for P(q)
52           if (r0.eq.0.0) then
53              pq=1.0
54           else
55              pq=(3.*(o-h*q)/(h*h*h))*(3.*(o-h*q)/(h*h*h))
56           endif
57*       Use S(q)=f(ds,eta) rather than S(q)=1+f(ds,eta)
58           ycal(j)=((prefac*pq)*((g*c*e*f)/x(j)))+bkgd
59        end do
60
61        return
62        end