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
00066
00073 class FFTmap_p1
00074 {
00075 public:
00076 enum FFTtype { Default, Measure, Estimate };
00077
00078 FFTmap_p1();
00080 FFTmap_p1( const FFTmap_p1& other ) { copy(other); }
00082 FFTmap_p1( const Grid_sampling& grid_sam, const FFTtype type = Default );
00084 const FFTmap_p1& operator =(const FFTmap_p1& other) { return copy(other); }
00086 void init( const Grid_sampling& grid_sam, const FFTtype type = Default );
00088 void reset();
00089
00091 const Grid_sampling& grid_real() const { return grid_sam_; }
00093 const Grid& grid_reci() const { return grid_reci_; }
00095 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()) ) ) ); }
00096
00098 void fft_h_to_x( const ftype& scale );
00100 void fft_x_to_h( const ftype& scale );
00101
00103 std::complex<ffttype> get_hkl( const HKL& hkl ) const;
00105 void set_hkl( const HKL& hkl, const std::complex<ffttype>& f );
00107 const std::complex<ffttype>& cplx_data( const Coord_grid& c ) const
00108 { return data_c[ grid_reci_.index( c ) ]; }
00110 std::complex<ffttype>& cplx_data( const Coord_grid& c )
00111 { return data_c[ grid_reci_.index( c ) ]; }
00113 const ffttype& real_data( const Coord_grid& c ) const
00114 { return data_r[ grid_real_.index( c ) ]; }
00116 ffttype& real_data( const Coord_grid& c )
00117 { return data_r[ grid_real_.index( c ) ]; }
00118
00120 static FFTtype& default_type() { return default_type_; }
00121
00122 void debug() const;
00123
00124 protected:
00125 const FFTmap_p1& copy( const FFTmap_p1& other );
00126
00127 enum FFTmode { NONE, RECI, REAL, OTHER };
00128
00129 FFTmode mode;
00130 FFTtype type_;
00131 Grid_sampling grid_sam_;
00132 Grid grid_reci_;
00133 Grid grid_real_;
00134 Grid grid_half_;
00135
00136 Matrix<char> req_kl, req_uv;
00137 std::vector<char> req_l, req_u;
00138
00139 std::vector<ffttype> datavec;
00140 ffttype* data_r;
00141 std::complex<ffttype>* data_c;
00142
00143 static FFTtype default_type_;
00144 };
00145
00146
00148
00159 class FFTmap : private FFTmap_p1
00160 {
00161 public:
00163 FFTmap();
00165 FFTmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00167 void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling grid_sam, const FFTtype type = Default );
00169 void reset();
00170
00172 const Cell& cell() const { return cell_; }
00174 const Spacegroup& spacegroup() const { return spacegroup_; }
00176 const Grid_sampling& grid_sampling() const { return FFTmap_p1::grid_real(); }
00178 void fft_h_to_x();
00180 void fft_x_to_h();
00181
00183 template<class T> void get_recip_data( const HKL& rfl, datatypes::F_phi<T>& fphi ) const;
00185 template<class T> void set_recip_data( const HKL& rfl, const datatypes::F_phi<T>& fphi );
00186
00188 template<class T> void get_real_data( const Coord_grid& c, T& datum ) const;
00190 template<class T> void set_real_data( const Coord_grid& c, const T& datum );
00191
00193 datatypes::F_phi<ffttype> get_recip_data( const HKL& rfl ) const;
00195 const ffttype& get_real_data( const Coord_grid& c ) const
00196 { return real_data(c.unit(grid_real())); }
00197
00199 template<class H, class X> void fft_rfl_to_map( const H& h, X& x );
00201 template<class H, class X> void fft_map_to_rfl( const X& x, H& h );
00202
00203 void debug() const;
00204
00205 protected:
00206 Cell cell_;
00207 Spacegroup spacegroup_;
00208 std::vector<Isymop> isymop;
00209 };
00210
00211
00212
00213
00214
00226 template<class H, class X> void FFTmap::fft_rfl_to_map( const H& h, X& x )
00227 {
00228
00229 reset();
00230
00231
00232 typename H::HKL_reference_index ih;
00233 for ( ih = h.first_data(); !ih.last(); h.next_data( ih ) )
00234 set_recip_data( ih.hkl(), h[ih] );
00235
00236
00237 fft_h_to_x();
00238
00239
00240 typename X::Map_reference_index ix;
00241 for ( ix = x.first(); !ix.last(); ix.next() )
00242 get_real_data( ix.coord(), x[ix] );
00243 }
00244
00245
00258 template<class H, class X> void FFTmap::fft_map_to_rfl( const X& x, H& h )
00259 {
00260
00261 reset();
00262
00263
00264 typename X::Map_reference_index ix;
00265 for ( ix = x.first(); !ix.last(); ix.next() )
00266 set_real_data( ix.coord(), x[ix] );
00267
00268
00269 fft_x_to_h();
00270
00271
00272 typename H::HKL_reference_index ih;
00273 for ( ih = h.first(); !ih.last(); ih.next() )
00274 get_recip_data( ih.hkl(), h[ih] );
00275 }
00276
00277
00278 }
00279
00280 #endif