| 1 | subroutine fractal_dm |
|---|
| 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,dm,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 | dm=p(3) |
|---|
| 16 | eta=p(4) |
|---|
| 17 | bkgd=p(5) |
|---|
| 18 | |
|---|
| 19 | if (dm.le.0) then |
|---|
| 20 | write(6,*)' ' |
|---|
| 21 | write(6,*)'The mass fractal dimension must be >0!' |
|---|
| 22 | write(6,*)'Will set to 3.' |
|---|
| 23 | dm=3.0 |
|---|
| 24 | end if |
|---|
| 25 | |
|---|
| 26 | if (dm.gt.6) then |
|---|
| 27 | write(6,*)' ' |
|---|
| 28 | write(6,*)'The mass fractal dimension must be <=6!' |
|---|
| 29 | write(6,*)'Will set to 3.' |
|---|
| 30 | dm=3.0 |
|---|
| 31 | end if |
|---|
| 32 | |
|---|
| 33 | a=dm-1.0 |
|---|
| 34 | call gammafunc(a,g) |
|---|
| 35 | |
|---|
| 36 | 10 format(' Dm=',f12.5,' Gamma(Dm-1)=',e12.5) |
|---|
| 37 | write(6,10)dm,g |
|---|
| 38 | |
|---|
| 39 | b=(1.0-dm)/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(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(dm,eta) rather than S(q)=1+f(dm,eta) |
|---|
| 58 | ycal(j)=((prefac*pq)*((g*c*e*f)/x(j)))+bkgd |
|---|
| 59 | end do |
|---|
| 60 | |
|---|
| 61 | return |
|---|
| 62 | end |
|---|
