Develop and Download Open Source Software

Powder auto-indexing software Conograph

Browse Subversion Repository

/[conograph]/Conograph/trunk/src/utility_data_structure/TreeLattice.cc

Diff of /Conograph/trunk/src/utility_data_structure/TreeLattice.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 32 by rtomiyasu, Fri Feb 22 04:51:31 2013 UTC revision 33 by rtomiyasu, Wed Sep 7 04:38:51 2016 UTC
# Line 25  THE SOFTWARE. Line 25  THE SOFTWARE.
25   *   *
26   */   */
27  #include <algorithm>  #include <algorithm>
28    #include <cmath>
29  #include "FracMat.hh"  #include "FracMat.hh"
30  #include "TreeLattice.hh"  #include "TreeLattice.hh"
31    
# Line 206  void TreeLattice::putBud(set<Bud>& budtr Line 207  void TreeLattice::putBud(set<Bud>& budtr
207  }  }
208    
209    
210  bool TreeLattice::putQuadraticForm(SymMat<Double>& Q, multimap<Int4, VecDat3<Int4> >& qindex_hkl) const  //bool TreeLattice::putQuadraticForm(SymMat<Double>& Q, multimap<Int4, VecDat3<Int4> >& qindex_hkl) const
211  {  //{
212          static const VecDat3<Int4> hkl100(1,0,0);  //      static const VecDat3<Int4> hkl100(1,0,0);
213          static const VecDat3<Int4> hkl010(0,1,0);  //      static const VecDat3<Int4> hkl010(0,1,0);
214          static const VecDat3<Int4> hkl_1_10(-1,-1,0);  //      static const VecDat3<Int4> hkl_1_10(-1,-1,0);
215    //
216          qindex_hkl.clear();  //      qindex_hkl.clear();
217          if( m_root == NULL ) return false;  //      if( m_root == NULL ) return false;
218    //
219          if(m_root->Left() >= 0)  //      if(m_root->Left() >= 0)
220          {  //      {
221                  qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root->Left(), hkl100 ) );  //              qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root->Left(), hkl100 ) );
222          }  //      }
223          if(m_root->Right() >= 0)  //      if(m_root->Right() >= 0)
224          {  //      {
225                  qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root->Right(), hkl010 ) );  //              qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root->Right(), hkl010 ) );
226          }  //      }
227    //
228          if( m_root_on_left != NULL )  //      if( m_root_on_left != NULL )
229          {  //      {
230                  if(m_root_on_left->Right() >= 0)  //              if(m_root_on_left->Right() >= 0)
231                  {  //              {
232                          qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root_on_left->Right(), hkl_1_10 ) );  //                      qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root_on_left->Right(), hkl_1_10 ) );
233                  }  //              }
234    //
235                  m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);  //              m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);
236                  m_root_on_left->putQuadraticForm(hkl010, hkl_1_10, qindex_hkl);  //              m_root_on_left->putQuadraticForm(hkl010, hkl_1_10, qindex_hkl);
237                  if( m_root_on_right != NULL ) m_root_on_right->putQuadraticForm(hkl_1_10, hkl100, qindex_hkl);  //              if( m_root_on_right != NULL ) m_root_on_right->putQuadraticForm(hkl_1_10, hkl100, qindex_hkl);
238          }  //      }
239          else if( m_root_on_right != NULL )  //      else if( m_root_on_right != NULL )
240          {  //      {
241                  if(m_root_on_right->Left() >= 0)  //              if(m_root_on_right->Left() >= 0)
242                  {  //              {
243                          qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root_on_right->Left(), hkl_1_10 ) );  //                      qindex_hkl.insert( multimap<Int4, VecDat3<Int4> >::value_type( m_root_on_right->Left(), hkl_1_10 ) );
244                  }  //              }
245    //
246                  m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);  //              m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);
247                  m_root_on_right->putQuadraticForm(hkl_1_10, hkl100, qindex_hkl);  //              m_root_on_right->putQuadraticForm(hkl_1_10, hkl100, qindex_hkl);
248          }  //      }
249          else  //      else
250          {  //      {
251                  m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);  //              m_root->putQuadraticForm(hkl100, hkl010, qindex_hkl);
252          }  //      }
253    //
254    //      if( qindex_hkl.size() < 3 ) return false;
255    //
256    //      const vector<QData>& qdata = VCData::putPeakQData();
257    //      NRMat<Int4> icoef(3,3);
258    //      NRVec<Double> qvalue(3);
259    //
260    //      multimap<Int4, VecDat3<Int4> >::const_iterator it = qindex_hkl.begin();
261    //      icoef[0][0] = it->second[0]*it->second[0];
262    //      icoef[0][1] = it->second[1]*it->second[1];
263    //      icoef[0][2] = it->second[0]*it->second[1]*2;
264    //      qvalue[0] = qdata[(it++)->first].q;
265    //      icoef[1][0] = it->second[0]*it->second[0];
266    //      icoef[1][1] = it->second[1]*it->second[1];
267    //      icoef[1][2] = it->second[0]*it->second[1]*2;
268    //      qvalue[1] = qdata[(it++)->first].q;
269    //      icoef[2][0] = it->second[0]*it->second[0];
270    //      icoef[2][1] = it->second[1]*it->second[1];
271    //      icoef[2][2] = it->second[0]*it->second[1]*2;
272    //      qvalue[2] = qdata[(it++)->first].q;
273    //
274    //      const FracMat inv_mat = FInverse3( icoef );
275    //      assert( Q.size() == 2 );
276    //
277    //      Q(0,0) = ( inv_mat.mat[0][0]*qvalue[0] + inv_mat.mat[0][1]*qvalue[1] + inv_mat.mat[0][2]*qvalue[2] ) / inv_mat.denom;
278    //      Q(1,1) = ( inv_mat.mat[1][0]*qvalue[0] + inv_mat.mat[1][1]*qvalue[1] + inv_mat.mat[1][2]*qvalue[2] ) / inv_mat.denom;
279    //      Q(0,1) = ( inv_mat.mat[2][0]*qvalue[0] + inv_mat.mat[2][1]*qvalue[1] + inv_mat.mat[2][2]*qvalue[2] ) / inv_mat.denom;
280    //
281    //      return true;
282    //}
283    
         if( qindex_hkl.size() < 3 ) return false;  
284    
285          const vector<QData>& qdata = VCData::putPeakQData();  void TreeLattice::putQuadraticForm(SymMat<VCData>& Q) const
286          NRMat<Int4> icoef(3,3);  {
287          NRVec<Double> qvalue(3);          assert(Q.size() == 2);
288            if( m_root == NULL ) return;
289          multimap<Int4, VecDat3<Int4> >::const_iterator it = qindex_hkl.begin();          if( m_root->IsBud() ) return;
         icoef[0][0] = it->second[0]*it->second[0];  
         icoef[0][1] = it->second[1]*it->second[1];  
         icoef[0][2] = it->second[0]*it->second[1]*2;  
         qvalue[0] = qdata[(it++)->first].q;  
         icoef[1][0] = it->second[0]*it->second[0];  
         icoef[1][1] = it->second[1]*it->second[1];  
         icoef[1][2] = it->second[0]*it->second[1]*2;  
         qvalue[1] = qdata[(it++)->first].q;  
         icoef[2][0] = it->second[0]*it->second[0];  
         icoef[2][1] = it->second[1]*it->second[1];  
         icoef[2][2] = it->second[0]*it->second[1]*2;  
         qvalue[2] = qdata[(it++)->first].q;  
   
         const FracMat inv_mat = FInverse3( icoef );  
         assert( Q.size() == 2 );  
   
         Q(0,0) = ( inv_mat.mat[0][0]*qvalue[0] + inv_mat.mat[0][1]*qvalue[1] + inv_mat.mat[0][2]*qvalue[2] ) / inv_mat.denom;  
         Q(1,1) = ( inv_mat.mat[1][0]*qvalue[0] + inv_mat.mat[1][1]*qvalue[1] + inv_mat.mat[1][2]*qvalue[2] ) / inv_mat.denom;  
         Q(0,1) = ( inv_mat.mat[2][0]*qvalue[0] + inv_mat.mat[2][1]*qvalue[1] + inv_mat.mat[2][2]*qvalue[2] ) / inv_mat.denom;  
290    
291          return true;          set<Bud> budtray;
292            this->putRootBuds(budtray);
293            Q(0,0) = budtray.begin()->Q1();
294            Q(1,1) = budtray.begin()->Q2();
295            Q(0,1) = ( Q(0,0) + Q(1,1) - budtray.begin()->Q3() )/2;
296  }  }
297    
298    
# Line 286  void TreeLattice::print(ostream& os, con Line 301  void TreeLattice::print(ostream& os, con
301          if( m_root == NULL ) return;          if( m_root == NULL ) return;
302          if( m_root->IsBud() ) return;          if( m_root->IsBud() ) return;
303                    
304      os << "* count : " << this->putCountOfQ() << endl;      os << "* Number of used peak positions: " << this->putCountOfQ() << endl;
305            SymMat<VCData> Q(2);
306          if( HeadIsTail )      this->putQuadraticForm(Q);
307          {  
308          os << "* This topograph contains a super-basis : ";      const Double det_Inv_Q = 1.0/( Q(0,0).Value()*Q(1,1).Value() - Q(0,1).Value()*Q(0,1).Value() );
309          set<Bud> budtray;      SymMat<Double> Inv_Q(2);
310          this->putRootBuds(budtray);      Inv_Q(0,0) =  Q(1,1).Value()*det_Inv_Q;
311                  os << "(" << budtray.begin()->Q3() << ",";      Inv_Q(1,1) =  Q(0,0).Value()*det_Inv_Q;
312          os << budtray.begin()->Q1() << ",";      Inv_Q(0,1) = -Q(0,1).Value()*det_Inv_Q;
313          os << budtray.begin()->Q2() << ")";  
314          os << endl;      Double a = sqrt(Q(0,0).Value());
315        Double b = sqrt(Q(1,1).Value());
316        os << "* (a*, b*, \\gamma*): (" << a << ", " << b << ", " << acos(Q(0,1).Value()/(a*b))*180.0/M_PI << ")\n";
317    
318        a = sqrt(Inv_Q(0,0));
319        b = sqrt(Inv_Q(1,1));
320        os << "* (a, b, \\gamma): (" << a << ", " << b << ", " << acos(Inv_Q(0,1)/(a*b))*180.0/M_PI << ")\n";
321    
322        if( HeadIsTail )
323            {
324          os.width(15);          os.width(15);
325          if( m_root->Upper() < 0 ) os << -1;          if( m_root->Upper() < 0 ) os << -1;
326          else os << m_root->Upper() + 1; // HeadIsTail is true.          else os << m_root->Upper() + 1; // HeadIsTail is true.
# Line 310  void TreeLattice::print(ostream& os, con Line 333  void TreeLattice::print(ostream& os, con
333          }          }
334          else if( m_root_on_left != NULL || m_root_on_right != NULL )          else if( m_root_on_left != NULL || m_root_on_right != NULL )
335      {      {
         os << "* This topograph contains a super-basis : ";  
         set<Bud> budtray;  
         this->putRootBuds(budtray);  
                 os << "(" << budtray.begin()->Q3() << ",";  
         os << budtray.begin()->Q1() << ",";  
         os << budtray.begin()->Q2() << ")";  
         os << endl;  
   
336          os.width(15);          os.width(15);
337          if( m_root_on_left != NULL )          if( m_root_on_left != NULL )
338          {          {
Legend:
Removed from v.32  
changed lines
  Added in v.33
Back to OSDN">Back to OSDN
 ViewVC Help
Powered by ViewVC 1.1.26  
About OSDN
Find Software
Develop Software
Help