00001
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef CLIPPER_FFTMAP
00046 #define CLIPPER_FFTMAP
00047
00048
00049 #include "coords.h"
00050
00051 #include <complex>
00052
00053
00054 namespace clipper
00055 {
00056
00057 typedef float ffttype;
00058
00059
00060 namespace datatypes {
00061 template<class T> class F_phi;
00062 }
00063
00064
00065 class FFTmap_base {
00066 public:
00067 enum FFTtype { Default, Measure, Estimate };
00068 protected:
00069 static Mutex mutex;
00070 };
00071
00073
00080 class FFTmap_p1 : public FFTmap_base
00081 {
00082 public:
00084 FFTmap_p1();
00086 FFTmap_p1( const FFTmap_p1& other ) { copy(other); }
00088 FFTmap_p1( const Grid_sampling& grid_sam, const FFTtype type = Default );
00090 const FFTmap_p1& operator =(const FFTmap_p1& other) { return copy(other); }
00092 void init( const Grid_sampling& grid_sam, const FFTtype type = Default );
00094 void reset();
00095
00097 const Grid_sampling& grid_real() const { return grid_sam_; }
00099 const Grid& grid_reci() const { return grid_reci_; }
00101 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()) ) ) ); }
00102
00104 void fft_h_to_x( const ftype& scale );
00106 void fft_x_to_h( const ftype& scale );
00107
00109 std::complex<ffttype> get_hkl( const HKL& hkl ) const;
00111 void set_hkl( const HKL& hkl, const std::complex<ffttype>& f );
00113 const std::complex<ffttype>& cplx_data( const Coord_grid& c ) const
00114 { return data_c[ grid_reci_.index( c ) ]; }
00116 std::complex<ffttype>& cplx_data( const Coord_grid& c )
00117 { return data_c[ grid_reci_.index( c ) ]; }
00119 const ffttype& real_data( const Coord_grid& c ) const
00120 { return data_r[ grid_real_.index( c ) ]; }
00122 ffttype& real_data( const Coord_grid& c )
00123 { return data_r[ grid_real_.index( c ) ]; }
00124
00126 static FFTtype& default_type() { return default_type_; }
00127
00128 void debug() const;
00129
00130 protected:
00131 const FFTmap_p1& copy( const FFTmap_p1& other );
00132
00133 enum FFTmode { NONE, RECI, REAL, OTHER };
00134
00135 FFTmode mode;
00136 FFTtype type_;
00137 Grid_sampling grid_sam_;
00138 Grid grid_reci_;
00139 Grid grid_real_;
00140 Grid grid_half_;
00141
00142 Matrix<char> req_kl, req_uv;
00143 std::vector<char> req_l, req_u;
00144
00145 std::vector<ffttype> datavec;
00146 ffttype* data_r;
00147 std::complex<ffttype>* data_c;
00148
00149 static FFTtype default_type_;
00150 };
00151
00152
00154
00165 class FFTmap : private FFTmap_p1
00166 {
00167 public:
00169 FFTmap();
00171 FFTmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00173 void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00175 void reset();
00176
00178 const Cell& cell() const { return cell_; }
00180 const Spacegroup& spacegroup() const { return spacegroup_; }
00182 const Grid_sampling& grid_sampling() const { return FFTmap_p1::grid_real(); }
00184 void fft_h_to_x();
00186 void fft_x_to_h();
00187
00189 template<class T> void get_recip_data( const HKL& rfl, datatypes::F_phi<T>& fphi ) const;
00191 template<class T> void set_recip_data( const HKL& rfl, const datatypes::F_phi<T>& fphi );
00192
00194 template<class T> void get_real_data( const Coord_grid& c, T& datum ) const;
00196 template<class T> void set_real_data( const Coord_grid& c, const T& datum );
00197
00199 datatypes::F_phi<ffttype> get_recip_data( const HKL& rfl ) const;
00201 const ffttype& get_real_data( const Coord_grid& c ) const
00202 { return real_data(c.unit(grid_real())); }
00203
00205 template<class H, class X> void fft_rfl_to_map( const H& h, X& x );
00207 template<class H, class X> void fft_map_to_rfl( const X& x, H& h );
00208
00209 void debug() const;
00210
00211 protected:
00212 Cell cell_;
00213 Spacegroup spacegroup_;
00214 std::vector<Isymop> isymop;
00215 };
00216
00217
00218
00219
00220
00232 template<class H, class X> void FFTmap::fft_rfl_to_map( const H& h, X& x )
00233 {
00234
00235 reset();
00236
00237
00238 typename H::HKL_reference_index ih;
00239 for ( ih = h.first_data(); !ih.last(); h.next_data( ih ) )
00240 set_recip_data( ih.hkl(), h[ih] );
00241
00242
00243 fft_h_to_x();
00244
00245
00246 typename X::Map_reference_index ix;
00247 for ( ix = x.first(); !ix.last(); ix.next() )
00248 get_real_data( ix.coord(), x[ix] );
00249 }
00250
00251
00264 template<class H, class X> void FFTmap::fft_map_to_rfl( const X& x, H& h )
00265 {
00266
00267 reset();
00268
00269
00270 typename X::Map_reference_index ix;
00271 for ( ix = x.first(); !ix.last(); ix.next() )
00272 set_real_data( ix.coord(), x[ix] );
00273
00274
00275 fft_x_to_h();
00276
00277
00278 typename H::HKL_reference_index ih;
00279 for ( ih = h.first(); !ih.last(); ih.next() )
00280 get_recip_data( ih.hkl(), h[ih] );
00281 }
00282
00283
00284 }
00285
00286 #endif