c c ---------------------------------------------------------------------- c subroutine epscen(ih,ik,il,centre,sysabs,eps,phicen,nsym,rsym) include 'fffear.fh' c EPSCEN - Centric and Epsilon tests c .. slow but reliable version integer ih,ik,il,nsym real eps,phicen,rsym(4,4,192) logical centre,sysabs integer i,jh,jk,jl c real rfrac external rfrac c centre=.false. sysabs=.false. eps=1.0 phicen=0.0 do 100 i=2,nsym c calculate symetry equivalent jh=nint(rsym(1,1,i)*ih+rsym(2,1,i)*ik+rsym(3,1,i)*il) jk=nint(rsym(1,2,i)*ih+rsym(2,2,i)*ik+rsym(3,2,i)*il) jl=nint(rsym(1,3,i)*ih+rsym(2,3,i)*ik+rsym(3,3,i)*il) c eps zone if rfl maps to self c sysabs if rfl maps to self with different phase c centric if maps to friedel opposite if (ih.eq.jh.and.ik.eq.jk.and.il.eq.jl) then eps=eps+1.0 if (cos(2.0*pi*(ih*rsym(1,4,i)+ik*rsym(2,4,i)+il*rsym(3,4,i))) + .le.0.999) sysabs=.true. endif if (ih.eq.-jh.and.ik.eq.-jk.and.il.eq.-jl) then centre=.true. phicen=pi*rfrac(ih*rsym(1,4,i)+ik*rsym(2,4,i)+il*rsym(3,4,i)) endif 100 continue return end c c ---------------------------------------------------------------------- c logical function hemi(ih,ik,il) include 'fffear.fh' c HEMI - reflection in triclinic hemisphere (CCP4) integer ih,ik,il hemi=(il.gt.0.or.(il.eq.0.and.(ih.gt.0.or.(ih.eq.0.and.ik.gt.0)))) return end c c ---------------------------------------------------------------------- c complex function phase(phi) include 'fffear.fh' c PHASE - calculate a complex phase factor real phi c phase=cmplx(cos(phi),sin(phi)) return end c c ---------------------------------------------------------------------- c real function phaseof(f) include 'fffear.fh' c PHASEOF - calculate phase of a complex number complex f c if (abs(f).gt.0.01) then phaseof=atan2(aimag(f),real(f)) else phaseof=0.0 endif return end c c ---------------------------------------------------------------------- c real function atanh(x) include 'fffear.fh' c arc hyperbolic tangent real x c if (x.gt.0.999999) then atanh=54.0 else if (x.lt.-0.999999) then atanh=-54.0 else atanh=log((1+x)/(1-x))/2.0 endif return end c c ---------------------------------------------------------------------- c real function sim(x) include 'fffear.fh' c Calculate Sim & Srinivasan non-centric figure of merit as c I1(X)/I0(X), where I1 and I0 are the modified 1st and zero c order bessel functions. c References: Sim, G. A. (1960) Acta Cryst. 13, 511-512; c Srinivasan, R. (1966) Acta Cryst. 20, 143-144; c Abramowitz & Stegun, Handbook of Mathematical Functions, 378. c This routine was obtained from W. Kabsch. real x,t c if (abs(x).gt.0.0001) then t=abs(x)/3.75 if (t.gt.1.0) then t=1.0/t sim=((((((((0.01787654-t*0.00420059)*t+ (-0.02895312))*t+ + 0.02282967)*t+ (-0.01031555))*t+0.00163801)*t+ + (-0.00362018))*t+ (-0.03988024))*t+0.39894228)* + sign(1.0,x)/ ((((((((-0.01647633+t*0.00392377)*t+ + 0.02635537)*t+ (-0.02057706))*t+0.00916281)*t+ + (-0.00157565))*t+0.00225319)*t+0.01328592)*t+0.39894228) else t=t*t sim=((((((t*0.00032411+0.00301532)*t+0.02658733)*t+ + 0.15084934)*t+0.51498869)*t+0.87890594)*t+0.5)*x/ + ((((((t*0.0045813+0.0360768)*t+0.2659732)*t+ + 1.2067492)*t+3.0899424)*t+3.5156229)*t+1.0) endif else sim=0.0 endif return end c c ---------------------------------------------------------------------- c real function siminv(fom) include 'fffear.fh' c SIMINV - inverse SIM function c slow but accurate version real fom,x,x1,x2,sim integer i external sim x1=0.0 x2=1.0 do 100 i=1,8 if (sim(x2).gt.fom) goto 110 x1=x2 x2=2.0*x2 100 continue siminv=x2 return 110 do 200 i=1,12 x=(x1+x2)/2 if (sim(x).gt.fom) then x2=x else x1=x endif 200 continue siminv=x return end c c ---------------------------------------------------------------------- c real function siminvq(fom) include 'fffear.fh' c SIMINV - inverse SIM function c tabulated 0..1 step 0.04 c value returned for fom=1.0 is truncated to x=54 real fom,tabmax integer ntab parameter (ntab=25,tabmax=1.0) real table(0:ntab) real rtabulated c external rtabulated c data table / 0.0000, 0.0801, 0.1605, 0.2417, 0.3242, 0.4083, + 0.4945, 0.5835, 0.6759, 0.7724, 0.8741, 0.9821, 1.0979, + 1.2235, 1.3616, 1.5157, 1.6913, 1.8964, 2.1436, 2.4549, + 2.8713, 3.4790, 4.4888, 6.5391,12.6432,54.3844 / c siminvq=rtabulated(fom,tabmax,table,ntab) return end c c ---------------------------------------------------------------------- c real function cumprc(e) include 'fffear.fh' c CUMPRC - cumulative probability distribution for Centric E magnitude real e,tabmax integer ntab parameter (ntab=25,tabmax=2.5) real table(0:ntab) real rtabulated c external rtabulated c data table / 0.0399, 0.1193, 0.1975, 0.2738, 0.3474, 0.4178, + 0.4845, 0.5469, 0.6049, 0.6581, 0.7065, 0.7501, 0.7889, + 0.8232, 0.8531, 0.8790, 0.9012, 0.9200, 0.9358, 0.9489, + 0.9597, 0.9685, 0.9756, 0.9813, 0.9858, 0.9893 / c cumprc=rtabulated(e,tabmax,table,ntab) return end c c ---------------------------------------------------------------------- c real function cumprnc(e) include 'fffear.fh' c CUMPRNC - cumulative probability distribution for NonCentric E magnitude real e,tabmax integer ntab parameter (ntab=25,tabmax=2.5) real table(0:ntab) real rtabulated c external rtabulated c data table / 0.0000, 0.0198, 0.0582, 0.1131, 0.1812, 0.2591, + 0.3428, 0.4286, 0.5130, 0.5931, 0.6666, 0.7322, 0.7891, + 0.8371, 0.8765, 0.9081, 0.9329, 0.9518, 0.9659, 0.9761, + 0.9835, 0.9886, 0.9920, 0.9944, 0.9959, 0.9968 / c cumprnc=rtabulated(e,tabmax,table,ntab) return end c c ---------------------------------------------------------------------- c real function rtabulated(x,tabmax,table,ntab) include 'fffear.fh' c RTABULATED - look up a tabulated function integer ntab,iy real x,y,tabmax,table(0:ntab) c y=x/tabmax*real(ntab) iy=min(max(int(y),0),ntab-1) rtabulated=table(iy)+(y-iy)*(table(iy+1)-table(iy)) return end c c ---------------------------------------------------------------------- c real function recipwt(s,r,modewt) include 'fffear.fh' c RECIPWT - weight for real space avaraging a map (from SFALL.FOR) c Fourier transform of function 1-R/r c S is 2.0*sin(theta)/lambda c R is radius of the sphere c MODEWT=1 for spherical top hat, MODEWT=2 for spherical cone integer modewt real s,r,x,phiwt,psi,psi2 c x = 6.2831853*s*r if (x.lt.0.01) then recipwt=1.0 else if (modewt.eq.0) then phiwt = (sin(x)-cos(x)*x)*3.0/x**3 recipwt = (phiwt) else if (modewt.eq.1) then phiwt = (sin(x)-cos(x)*x)*3.0/x**3 psi = (2.0*x*sin(x)-(x**2-2.0)*cos(x)-2.0)*3.0/x**4 recipwt = 4.0*(phiwt - psi) else phiwt = (sin(x)-cos(x)*x)*3.0/x**3 psi2 = ((x**2*3.0-6.0)*sin(x)-(x**3-6.0*x)*cos(x))*3.0/x**5 recipwt = 2.5*(phiwt - psi2) endif c return end c c ---------------------------------------------------------------------- c