2
Header file for P1 fft map
4
//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
6
//L This library is free software and is distributed under the terms
7
//L and conditions of version 2.1 of the GNU Lesser General Public
8
//L Licence (LGPL) with the following additional clause:
10
//L `You may also combine or link a "work that uses the Library" to
11
//L produce a work containing portions of the Library, and distribute
12
//L that work under terms of your choice, provided that you give
13
//L prominent notice with each copy of the work that the specified
14
//L version of the Library is used in it, and that you include or
15
//L provide public access to the complete corresponding
16
//L machine-readable source code for the Library including whatever
17
//L changes were used in the work. (i.e. If you make changes to the
18
//L Library you must distribute those, but you do not need to
19
//L distribute source or object code to those portions of the work
20
//L not covered by this licence.)'
22
//L Note that this clause grants an additional right and does not impose
23
//L any additional restriction, and so does not affect compatibility
24
//L with the GNU General Public Licence (GPL). If you wish to negotiate
25
//L other terms, please contact the maintainer.
27
//L You can redistribute it and/or modify the library under the terms of
28
//L the GNU Lesser General Public License as published by the Free Software
29
//L Foundation; either version 2.1 of the License, or (at your option) any
32
//L This library is distributed in the hope that it will be useful, but
33
//L WITHOUT ANY WARRANTY; without even the implied warranty of
34
//L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35
//L Lesser General Public License for more details.
37
//L You should have received a copy of the CCP4 licence and/or GNU
38
//L Lesser General Public License along with this library; if not, write
39
//L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
40
//L The GNU Lesser General Public can also be obtained by writing to the
41
//L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
45
#ifndef CLIPPER_FFTMAP
46
#define CLIPPER_FFTMAP
57
typedef float ffttype;
61
template<class T> class F_phi;
65
//! FFTmap_p1: low level P1 map used for calculating FFTs
66
/*! This is a pure real P1 map, with an extra section in reciprocal
67
space to allow generation of the full set of resiprocal space
68
magnitudes. Access is by Coord_grid in both spaces, and indices
69
must be non-negative and in range. The first and last sections
70
along the half-length direction only have half the elements
71
stored, the contents of the other half is ignored.
76
enum FFTtype { Default, Measure, Estimate }; //!< optimisation options
80
FFTmap_p1( const FFTmap_p1& other ) { copy(other); }
81
//! Constructor: takes grid
82
FFTmap_p1( const Grid_sampling& grid_sam, const FFTtype type = Default );
83
//! Assignment operator
84
const FFTmap_p1& operator =(const FFTmap_p1& other) { return copy(other); }
85
//! initialiser: takes grid
86
void init( const Grid_sampling& grid_sam, const FFTtype type = Default );
90
//! Return real space grid.
91
const Grid_sampling& grid_real() const { return grid_sam_; }
92
//! Return reciprocal space grid (i.e. half real grid + 1 section).
93
const Grid& grid_reci() const { return grid_reci_; }
94
//! Test whether a coordinate is in the valid part of the recip. grid.
95
bool uniq_reci( const Coord_grid& c ) const { return ( (c.w()>0 && c.w()<grid_half_.nw()) || (c.w()<=grid_half_.nw() && ( (c.v()>0 && c.v()<grid_half_.nv()) || (c.v()<=grid_half_.nv() && c.u()<=grid_half_.nu()) ) ) ); }
97
//! Transform to real space
98
void fft_h_to_x( const ftype& scale );
99
//! Transform to reciprocal space
100
void fft_x_to_h( const ftype& scale );
102
//! get reciprocal space data: slow form with hemisphere check
103
std::complex<ffttype> get_hkl( const HKL& hkl ) const;
104
//! set reciprocal space data: slow form with hemisphere check
105
void set_hkl( const HKL& hkl, const std::complex<ffttype>& f );
106
//! get reciprocal space data
107
const std::complex<ffttype>& cplx_data( const Coord_grid& c ) const
108
{ return data_c[ grid_reci_.index( c ) ]; }
109
//! set reciprocal space data
110
std::complex<ffttype>& cplx_data( const Coord_grid& c )
111
{ return data_c[ grid_reci_.index( c ) ]; }
112
//! get real space data
113
const ffttype& real_data( const Coord_grid& c ) const
114
{ return data_r[ grid_real_.index( c ) ]; }
115
//! set real space data
116
ffttype& real_data( const Coord_grid& c )
117
{ return data_r[ grid_real_.index( c ) ]; }
119
//! set/get default optimisation type
120
static FFTtype& default_type() { return default_type_; }
125
const FFTmap_p1& copy( const FFTmap_p1& other ); //!< copy function
127
enum FFTmode { NONE, RECI, REAL, OTHER }; //!< space enumeration
129
FFTmode mode; //!< real or reciprocal space?
130
FFTtype type_; //!< optimisation options
131
Grid_sampling grid_sam_; //!< unit cell grid
132
Grid grid_reci_; //!< reciprocal space grid
133
Grid grid_real_; //!< real space grid
134
Grid grid_half_; //!< half grid (for marking unique)
136
Matrix<char> req_kl, req_uv; //!< reci section lookup
137
std::vector<char> req_l, req_u; //!< real section lookup
139
std::vector<ffttype> datavec; //!< vector for the data
140
ffttype* data_r; //!< pointer to real data
141
std::complex<ffttype>* data_c; //!< pointer to complex data
143
static FFTtype default_type_; //!< default optimisation options
147
//! FFTmap: P1 map with symmetry used for calculating FFTs
148
/*! The FFTmap is represented in P1 in memory. However, it also has
149
a spacegroup, and the contained data remains consistent with this
150
spacegroup at all times. It has three states - unassigned,
151
real-space, and reciprocal space. In real space it contains real
152
map data. In reciprocal space it holds a hemisphere of complex
153
structure factors, with the Friedels duplicated on the zero
156
The user should be able to ignore all the issues of spacegroup
157
symmetry, Friedel opposites, and storage order.
159
class FFTmap : private FFTmap_p1
164
//! Constructor: takes spacegroup, cell, grid
165
FFTmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
167
void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
172
const Cell& cell() const { return cell_; }
173
//! get the spacegroup
174
const Spacegroup& spacegroup() const { return spacegroup_; }
175
//! get the cell grid
176
const Grid_sampling& grid_sampling() const { return FFTmap_p1::grid_real(); }
177
//! Transform to real space
179
//! Transform to reciprocal space
182
//! get reciprocal space data
183
template<class T> void get_recip_data( const HKL& rfl, datatypes::F_phi<T>& fphi ) const;
184
//! set reciprocal space data
185
template<class T> void set_recip_data( const HKL& rfl, const datatypes::F_phi<T>& fphi );
187
//! get real space data
188
template<class T> void get_real_data( const Coord_grid& c, T& datum ) const;
189
//! set real space data
190
template<class T> void set_real_data( const Coord_grid& c, const T& datum );
192
//! get reciprocal space data (No error checking)
193
datatypes::F_phi<ffttype> get_recip_data( const HKL& rfl ) const;
194
//! get real space data (No error checking)
195
const ffttype& get_real_data( const Coord_grid& c ) const
196
{ return real_data(c.unit(grid_real())); }
198
//! calculate map-like object from reflection-like object
199
template<class H, class X> void fft_rfl_to_map( const H& h, X& x );
200
//! calculate reflection-like object from map-like object
201
template<class H, class X> void fft_map_to_rfl( const X& x, H& h );
206
Cell cell_; //!< unit cell
207
Spacegroup spacegroup_; //!< spacegroup
208
std::vector<Isymop> isymop; //!< Integerised symops
212
// template implementations
215
/*! Fill this FFTmap object from a reflection object, transform it,
216
and fill the given map object from the FFTmap. This will work for
217
any reflection data object which implements a HKL_reference_index,
218
and every map data object which implements a Map_reference_index.
220
For the results to be sensible, the spacegroup, cell and grids
221
should match. (The map will be zeroed if necessary).
223
\param h The source reflection data object.
224
\param x The target map object.
226
template<class H, class X> void FFTmap::fft_rfl_to_map( const H& h, X& x )
231
// copy from reflection data
232
typename H::HKL_reference_index ih;
233
for ( ih = h.first_data(); !ih.last(); h.next_data( ih ) )
234
set_recip_data( ih.hkl(), h[ih] );
239
// and copy into the map
240
typename X::Map_reference_index ix;
241
for ( ix = x.first(); !ix.last(); ix.next() )
242
get_real_data( ix.coord(), x[ix] );
246
/*! Fill this FFTmap object from a map object, transform it, and
247
fill the given reflection object from the FFTmap. This will work
248
for any reflection data object which implements a
249
HKL_reference_index, and every map data object which implements a
252
For the results to be sensible, the spacegroup, cell and grids
253
should match. (The map will be zeroed if necessary).
255
\param x The source map object.
256
\param h The target reflection data object.
258
template<class H, class X> void FFTmap::fft_map_to_rfl( const X& x, H& h )
264
typename X::Map_reference_index ix;
265
for ( ix = x.first(); !ix.last(); ix.next() )
266
set_real_data( ix.coord(), x[ix] );
272
typename H::HKL_reference_index ih;
273
for ( ih = h.first(); !ih.last(); ih.next() )
274
get_recip_data( ih.hkl(), h[ih] );
278
} // namespace clipper