| 1 | subroutine fractal_qton |
|---|
| 2 | |
|---|
| 3 | * See Schmidt, J Appl Cryst, 24, (1991), 414-435, Eqn (19) |
|---|
| 4 | * See Hurd, Schaefer & Martin, 35, (1987), 2361-2364 |
|---|
| 5 | |
|---|
| 6 | include 'frills_sources:function.inc' |
|---|
| 7 | |
|---|
| 8 | real*8 prefac,dm,ds,big_rg_small_rg,bkgd |
|---|
| 9 | real*8 a,b,c,d,e,f,g,h,i |
|---|
| 10 | |
|---|
| 11 | prefac=p(1) |
|---|
| 12 | dm=p(2) |
|---|
| 13 | big_rg=p(3) |
|---|
| 14 | ds=p(4) |
|---|
| 15 | small_rg=p(5) |
|---|
| 16 | bkgd=p(6) |
|---|
| 17 | |
|---|
| 18 | if (dm.le.0) then |
|---|
| 19 | write(6,*)'The mass fractal dimension must be >0!' |
|---|
| 20 | write(6,*)'Will set to 3.' |
|---|
| 21 | dm=3.0 |
|---|
| 22 | end if |
|---|
| 23 | |
|---|
| 24 | if (dm.gt.6) then |
|---|
| 25 | write(6,*)'The mass fractal dimension must be <=6!' |
|---|
| 26 | write(6,*)'Will set to 3.' |
|---|
| 27 | dm=3.0 |
|---|
| 28 | end if |
|---|
| 29 | |
|---|
| 30 | if (ds.le.0) then |
|---|
| 31 | write(6,*)'The surface fractal dimension must be >0!' |
|---|
| 32 | write(6,*)'Will set to 2.' |
|---|
| 33 | ds=2.0 |
|---|
| 34 | end if |
|---|
| 35 | |
|---|
| 36 | if (ds.gt.6) then |
|---|
| 37 | write(6,*)'The surface fractal dimension must be <=6!' |
|---|
| 38 | write(6,*)'Will set to 2.' |
|---|
| 39 | ds=2.0 |
|---|
| 40 | end if |
|---|
| 41 | |
|---|
| 42 | a=dm/2.0 |
|---|
| 43 | b=(big_rg*big_rg)/(3.0*a) |
|---|
| 44 | |
|---|
| 45 | 10 format(' Will set to (6-Dm)=',f12.5) |
|---|
| 46 | * IF C GOES NEGATIVE, SANDRA WILL CRASH WITH UNDEFINED EXPONENTIATION |
|---|
| 47 | * SO (Ds+Dm)<=6 |
|---|
| 48 | * c=((ds-6.0)/-2.0)-a |
|---|
| 49 | if ((ds.gt.(6.0-dm)).and.(small_rg.gt.0.0)) then |
|---|
| 50 | ds=6.0-dm |
|---|
| 51 | write(6,*)'The surface fractal dimension must be <=(6-Dm)!' |
|---|
| 52 | write(6,10)ds |
|---|
| 53 | endif |
|---|
| 54 | * c=(ds-6.0+dm)/-2.0 |
|---|
| 55 | c=(6.0-ds-dm)/2.0 |
|---|
| 56 | |
|---|
| 57 | * IF C=0 THEN SANDRA WILL CRASH WITH A FLOATING DIVIDE BY ZERO |
|---|
| 58 | if (c.eq.0.) then |
|---|
| 59 | d=1.0E+37 |
|---|
| 60 | else |
|---|
| 61 | d=(small_rg*small_rg)/(3.0*c) |
|---|
| 62 | endif |
|---|
| 63 | |
|---|
| 64 | do j=1,nv |
|---|
| 65 | e=x(j)*x(j)*b |
|---|
| 66 | f=x(j)*x(j)*d |
|---|
| 67 | g=(1.+e)**a |
|---|
| 68 | h=(1.+f)**c |
|---|
| 69 | i=g*h |
|---|
| 70 | ycal(j)=(prefac/i)+bkgd |
|---|
| 71 | end do |
|---|
| 72 | |
|---|
| 73 | return |
|---|
| 74 | end |
|---|
