1 /* clipper_util.cpp: implementation file for clipper helper functions */
2 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
3 //L
4 //L This library is free software and is distributed under the terms
5 //L and conditions of version 2.1 of the GNU Lesser General Public
6 //L Licence (LGPL) with the following additional clause:
7 //L
8 //L `You may also combine or link a "work that uses the Library" to
9 //L produce a work containing portions of the Library, and distribute
10 //L that work under terms of your choice, provided that you give
11 //L prominent notice with each copy of the work that the specified
12 //L version of the Library is used in it, and that you include or
13 //L provide public access to the complete corresponding
14 //L machine-readable source code for the Library including whatever
15 //L changes were used in the work. (i.e. If you make changes to the
16 //L Library you must distribute those, but you do not need to
17 //L distribute source or object code to those portions of the work
18 //L not covered by this licence.)'
19 //L
20 //L Note that this clause grants an additional right and does not impose
21 //L any additional restriction, and so does not affect compatibility
22 //L with the GNU General Public Licence (GPL). If you wish to negotiate
23 //L other terms, please contact the maintainer.
24 //L
25 //L You can redistribute it and/or modify the library under the terms of
26 //L the GNU Lesser General Public License as published by the Free Software
27 //L Foundation; either version 2.1 of the License, or (at your option) any
28 //L later version.
29 //L
30 //L This library is distributed in the hope that it will be useful, but
31 //L WITHOUT ANY WARRANTY; without even the implied warranty of
32 //L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
33 //L Lesser General Public License for more details.
34 //L
35 //L You should have received a copy of the CCP4 licence and/or GNU
36 //L Lesser General Public License along with this library; if not, write
37 //L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
38 //L The GNU Lesser General Public can also be obtained by writing to the
39 //L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
40 //L MA 02111-1307 USA
41
42
43 #include "clipper_util.h"
44
45
46 namespace clipper {
47
48 float Util::nanf_; //!< declare static
49 double Util::nand_; //!< declare static
50 ftype Util::nan_ ; //!< declare static
51 ftype Util::onepi_ = M_PI; //!< one*pi
52 ftype Util::twopi_ = 2.0*M_PI; //!< two*pi
53 ftype Util::twopi2_ = 2.0*M_PI*M_PI; //!< two*pi*pi
54 ftype Util::eightpi2_ = 8.0*M_PI*M_PI; //!< two*pi*pi
55 ftype Util::d2rad_ = M_PI/180.0; //!< degree-radian conversion
56 ftype Util::sim_a = 1.639294; //!< sim fn param
57 ftype Util::sim_b = 3.553967; //!< sim fn param
58 ftype Util::sim_c = 2.228716; //!< sim fn param
59 ftype Util::sim_d = 3.524142; //!< sim fn param
60 ftype Util::sim_e = 7.107935; //!< sim fn param
61 ftype Util::sim_A = -1.28173889903; //!< invsim fn param
62 ftype Util::sim_B = 0.69231689903; //!< invsim fn param
63 ftype Util::sim_C = -1.33099462667; //!< invsim fn param
64 ftype Util::sim_g = 2.13643992379; //!< invsim fn param
65 ftype Util::sim_p = 0.04613803811; //!< invsim fn param
66 ftype Util::sim_q = 1.82167089029; //!< invsim fn param
67 ftype Util::sim_r = -0.74817947490; //!< invsim fn param
68
69
Util()70 Util::Util()
71 {
72 ((U32*)&nanf_)->i = CLIPPER_NULL_MASK_32;
73 ((U64*)&nand_)->i = CLIPPER_NULL_MASK_64;
74 if ( sizeof(ftype) == 4 ) {
75 ((U32*)&nan_)->i = CLIPPER_NULL_MASK_32;
76 } else if ( sizeof(ftype) == 8 ) {
77 ((U64*)&nan_)->i = CLIPPER_NULL_MASK_64;
78 } else {
79 /* fail on build */;
80 }
81 }
82
83 /*! \param x The argument \return I1(x)/I0(x) */
sim(const ftype & x)84 ftype Util::sim( const ftype& x )
85 {
86 if (x >= 0.0) return (((x + sim_a)*x + sim_b)*x)
87 / (((x + sim_c)*x + sim_d)*x + sim_e);
88 else return -(-(-(-x + sim_a)*x + sim_b)*x)
89 / (-(-(-x + sim_c)*x + sim_d)*x + sim_e);
90 }
91
92 /*! \param x I1(y)/I0(y) \return y */
invsim(const ftype & x)93 ftype Util::invsim( const ftype& x )
94 {
95 ftype x0 = fabs(x);
96 ftype a0 = -7.107935*x0;
97 ftype a1 = 3.553967-3.524142*x0;
98 ftype a2 = 1.639294-2.228716*x0;
99 ftype a3 = 1.0-x0;
100 ftype w = a2/(3.0*a3);
101 ftype p = a1/(3.0*a3)-w*w;
102 ftype q = -w*w*w+0.5*(a1*w-a0)/a3;
103 ftype d = sqrt(q*q+p*p*p);
104 ftype q1 = q + d;
105 ftype q2 = q - d;
106 ftype r1 = pow(fabs(q1), 1.0/3.0);
107 ftype r2 = pow(fabs(q2), 1.0/3.0);
108 if (x >= 0.0) return (((q1>0.0)? r1 : -r1) + ((q2>0.0)? r2 : -r2) - w);
109 else return -(((q1>0.0)? r1 : -r1) + ((q2>0.0)? r2 : -r2) - w);
110 }
111
sim_integ(const ftype & x0)112 ftype Util::sim_integ( const ftype& x0 )
113 {
114 const ftype x = fabs(x0);
115 const ftype z = (x+sim_p)/sim_q;
116 return sim_A*log(x+sim_g) + 0.5*sim_B*log(z*z+1.0) + sim_r*atan(z) + x + 1.0;
117 }
118
sim_deriv(const ftype & x)119 ftype Util::sim_deriv( const ftype& x )
120 {
121 if (x >= 0.0) return (((((sim_c-sim_a)*x+(2.0*sim_d-2.0*sim_b))*x+(3.0*sim_e+sim_a*sim_d-sim_b*sim_c))*x+(2.0*sim_a*sim_e))*x+(sim_b*sim_e)) / pow( (((x + sim_c)*x + sim_d)*x + sim_e), 2.0 );
122 else return (-(-(-(-(sim_c-sim_a)*x+(2.0*sim_d-2.0*sim_b))*x+(3.0*sim_e+sim_a*sim_d-sim_b*sim_c))*x+(2.0*sim_a*sim_e))*x+(sim_b*sim_e)) / pow( (-(-(-x + sim_c)*x + sim_d)*x + sim_e), 2.0 );
123 }
124
sim_deriv_recur(const ftype & x0)125 ftype Util::sim_deriv_recur( const ftype& x0 )
126 {
127 const ftype x = fabs(x0);
128 const ftype m = sim(x);
129 if ( x > 1.0e-4 )
130 return ( -m/x + ( 1.0 - m*m ) );
131 else
132 return ( 0.5 - m*m );
133 }
134
bessel_i0(const ftype & x0)135 ftype Util::bessel_i0( const ftype& x0 )
136 {
137 ftype i0 = 0.0, t, x;
138 x = fabs(x0);
139 t=x/3.75;
140 if (t < 1.0) {
141 t=t*t;
142 i0= ((((((t*0.0045813+0.0360768)*t+0.2659732)*t+
143 1.2067492)*t+3.0899424)*t+3.5156229)*t+1.0);
144 } else {
145 i0= (1.0/sqrt(x))*((((((((t*0.00392377-0.01647633)*t+
146 0.02635537)*t-0.02057706)*t+0.00916281)*t-
147 0.00157565)*t+0.00225319)*t+0.01328592)*t+0.39894228)*exp(x);
148 }
149 return i0;
150 }
151
152 /*! \param x Angle in degrees \return Angle in radians */
d2rad(const ftype & x)153 ftype Util::d2rad( const ftype& x )
154 { return x*d2rad_; }
155
156 /*! \param x Angle in radians \return Angle in degrees */
rad2d(const ftype & x)157 ftype Util::rad2d( const ftype& x )
158 { return x/d2rad_; }
159
160
161 // template instantiations
162 template ftype32 Util::max<ftype32>( const ftype32& a, const ftype32& b );
163 template ftype32 Util::min<ftype32>( const ftype32& a, const ftype32& b );
164 template ftype32 Util::bound<ftype32>( const ftype32& a, const ftype32& b, const ftype32& c );
165 template void Util::swap<ftype32>( ftype32& a, ftype32& b );
166 template void Util::swap<ftype32>( ftype32& a, ftype32& b, ftype32& c );
167 template ftype32 Util::sqr<ftype32>( const ftype32& a );
168 template ftype32 Util::isqrt<ftype32>( const ftype32& a );
169
170 template ftype64 Util::max<ftype64>( const ftype64& a, const ftype64& b );
171 template ftype64 Util::min<ftype64>( const ftype64& a, const ftype64& b );
172 template ftype64 Util::bound<ftype64>( const ftype64& a, const ftype64& b, const ftype64& c );
173 template void Util::swap<ftype64>( ftype64& a, ftype64& b );
174 template void Util::swap<ftype64>( ftype64& a, ftype64& b, ftype64& c );
175 template ftype64 Util::sqr<ftype64>( const ftype64& a );
176 template ftype64 Util::isqrt<ftype64>( const ftype64& a );
177
178 template int Util::max<int>( const int& a, const int& b );
179 template int Util::min<int>( const int& a, const int& b );
180 template int Util::bound<int>( const int& a, const int& b, const int& c );
181 template void Util::swap<int>( int& a, int& b );
182 template void Util::swap<int>( int& a, int& b, int& c );
183 template int Util::sqr<int>( const int& a );
184 template int Util::isqrt<int>( const int& a );
185
186
187 } // namespace clipper
188