CLASS  1.1
src/IsotopicVector.cxx
Aller à la documentation de ce fichier.
00001 #include "IsotopicVector.hxx"
00002 #include "ZAI.hxx"
00003 #include "LogFile.hxx"
00004 #include "Defines.hxx"
00005 
00006 
00007 #include <cmath>
00008 #include <iostream>
00009 #include <fstream>
00010 #include <string>
00011 #include <algorithm>
00012 
00013 //________________________________________________________________________
00014 //________________________________________________________________________
00015 //
00016 //
00017 //
00018 //                              IsotopicVector
00019 //
00020 //
00021 //________________________________________________________________________
00022 //________________________________________________________________________
00023 
00024 
00025 
00026 
00027 //________________________________________________________________________
00028 //__________________________Operator Overlaoding__________________________
00029 //________________________________________________________________________
00030 
00031 //____________________________General Operator____________________________
00032 //________________________________________________________________________
00033 
00034 ClassImp(IsotopicVector)
00035 
00036 double  Norme(IsotopicVector IV1)
00037 {
00038 DBGL;
00039         IsotopicVector IV; 
00040         return Distance(IV1, IV);
00041 DBGL;
00042 }
00043 
00044 double Distance(IsotopicVector IV1, IsotopicVector IV2 )
00045 {
00046 DBGL;
00047         double d2 = 0;
00048         
00049         IsotopicVector IVtmp = IV1 + IV2;
00050         map<ZAI ,double> IVtmpIsotopicQuantity = IVtmp.GetIsotopicQuantity();
00051         
00052         for(map<ZAI ,double >::iterator it = IVtmpIsotopicQuantity.begin(); it != IVtmpIsotopicQuantity.end(); it++)
00053         {
00054                 double Z1 = IV1.GetZAIIsotopicQuantity((*it).first);
00055                 double Z2 = IV2.GetZAIIsotopicQuantity((*it).first);
00056                 d2 += pow(Z1-Z2 , 2 );
00057         }
00058         
00059 DBGL;
00060         return sqrt(d2);
00061 }
00062 
00063 
00064 IsotopicVector operator+(IsotopicVector const& IVa, IsotopicVector const& IVb)
00065 {
00066         IsotopicVector IVtmp;
00067         IVtmp = IVa;
00068         return IVtmp += IVb;
00069 }
00070 
00071 //________________________________________________________________________
00072 IsotopicVector operator-(IsotopicVector const& IVa, IsotopicVector const& IVb)
00073 {
00074         IsotopicVector IVtmp;
00075         IVtmp = IVa;
00076         return IVtmp -= IVb;
00077 }
00078 
00079 
00080 //________________________________________________________________________
00081 IsotopicVector operator*(ZAI const& zai, double F)
00082 {
00083         IsotopicVector IVtmp;
00084                 
00085         IVtmp.Add( zai, F);
00086                 
00087         return IVtmp;
00088 }
00089 //________________________________________________________________________
00090 IsotopicVector operator/(ZAI const& zai, double F)
00091 {
00092         IsotopicVector IVtmp;
00093                 
00094         IVtmp.Add( zai, 1./F);
00095                 
00096         return IVtmp;
00097 }
00098 
00099 
00100 //________________________________________________________________________
00101 IsotopicVector operator*(IsotopicVector const& IVA, double F)
00102 {
00103 DBGL;
00104         IsotopicVector IV = IVA;
00105         IV.Multiply(F);
00106         return IV;
00107 }
00108 
00109 //________________________________________________________________________
00110 IsotopicVector operator/(IsotopicVector const& IVA, double F)
00111 {
00112 DBGL;
00113         IsotopicVector IV = IVA;
00114         IV.Multiply(1./F);
00115         return IV;
00116 }
00117 
00118 
00119 //____________________________InClass Operator____________________________
00120 
00121 //________________________________________________________________________
00122 IsotopicVector& IsotopicVector::operator+=(const IsotopicVector& IVa)
00123 {
00124         Add(IVa);
00125         return *this;
00126 }
00127 
00128 //________________________________________________________________________
00129 IsotopicVector& IsotopicVector::operator-=(const IsotopicVector& IVa)
00130 {
00131         Remove(IVa);
00132         return *this;
00133 }
00134 
00135 
00136 
00137 //________________________________________________________________________
00138 //________________________Constructor & Destructor________________________
00139 //________________________________________________________________________
00140 IsotopicVector::IsotopicVector()
00141 {
00142 DBGL;
00143 DBGL;
00144 }
00145 
00146 
00147 //________________________________________________________________________
00148 IsotopicVector::~IsotopicVector()
00149 {
00150 DBGL;
00151 }
00152 
00153 
00154 
00155 //________________________________________________________________________
00156 //_____________________________General Method_____________________________
00157 //________________________________________________________________________
00158 void IsotopicVector::Clear()
00159 {
00160 DBGL;
00161         fIsotopicQuantityNeeded.clear();
00162         fIsotopicQuantity.clear();
00163 DBGL;
00164 }
00165 //________________________________________________________________________
00166 void IsotopicVector::ClearNeed()
00167 {
00168 DBGL;
00169         fIsotopicQuantityNeeded.clear();
00170 DBGL;
00171 }
00172 
00173 //________________________________________________________________________
00174 void IsotopicVector::Multiply(double factor)
00175 {
00176 DBGL;
00177         for(map<ZAI ,double >::iterator it = fIsotopicQuantity.begin(); it != fIsotopicQuantity.end(); it++)
00178                 (*it).second = (*it).second * factor;
00179 DBGL;
00180         
00181 }
00182 
00183 //________________________________________________________________________
00184 void IsotopicVector::Add(const ZAI& zai, double quantity)
00185 {
00186 DBGL;
00187         if( ceil(quantity*1e6) - quantity*1e6 >  quantity*1e6 - floor(quantity*1e6) )
00188                 quantity = floor(quantity*1e6)*1/1e6;
00189         else    quantity = ceil(quantity*1e6)*1/1e6;
00190         
00191         pair<map<ZAI, double>::iterator, bool> IResult;
00192         if(quantity > 0)
00193         {
00194                 IResult = fIsotopicQuantity.insert( pair<ZAI ,double>(zai, quantity));
00195                 if(IResult.second == false)
00196                         IResult.first->second += quantity;
00197         }
00198 
00199 DBGL;
00200 }
00201 //________________________________________________________________________
00202 
00203 void IsotopicVector::Add(const IsotopicVector& isotopicvector)
00204 {
00205 DBGL;
00206         map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
00207         
00208         for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
00209                 Add( (*it).first, (*it).second);
00210 DBGL;
00211 
00212 }
00213 //________________________________________________________________________
00214 
00215 void IsotopicVector::Add(const map<ZAI ,double>& quantity)
00216 {
00217 DBGL;
00218         map<ZAI ,double> isotopicquantity = quantity;
00219         
00220         for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
00221                 Add( (*it).first, (*it).second);
00222 DBGL;
00223 
00224 }
00225 
00226 
00227 //________________________________________________________________________
00228 void IsotopicVector::Remove(const ZAI& zai, double quantity)
00229 {
00230 DBGL;
00231 
00232         map<ZAI ,double>::iterator it;
00233         it = fIsotopicQuantity.find(zai);
00234 
00235         if( ceil(quantity*1e6) - quantity*1e6 >  quantity*1e6 - floor(quantity*1e6) )
00236                 quantity = floor(quantity*1e6)*1/1e6;
00237         else    quantity = ceil(quantity*1e6)*1/1e6;
00238         
00239         if(quantity > 0)
00240         {       
00241                 if ( it != fIsotopicQuantity.end() ) 
00242                 {
00243                         if (it->second > quantity) 
00244                                 it->second = it->second - quantity;
00245                         else 
00246                         {
00247                                 Need(zai, quantity - it->second );
00248                                 it->second = 0;
00249                         }
00250                 }
00251                 else
00252                 {
00253                         Need(zai, quantity);
00254                 }
00255         }
00256 DBGL;
00257 
00258 }
00259 
00260 //________________________________________________________________________
00261 void IsotopicVector::Remove(const IsotopicVector& isotopicvector)
00262 {
00263 DBGL;
00264         map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
00265         for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
00266                 Remove( (*it).first, (*it).second);
00267 DBGL;
00268 }
00269 
00270 //________________________________________________________________________
00271 void IsotopicVector::Need(const ZAI& zai, double quantity)
00272 {
00273 DBGL;
00274         if( ceil(quantity*1e6) - quantity*1e6 >  quantity*1e6 - floor(quantity*1e6) )
00275                 quantity = floor(quantity*1e6)*1/1e6;
00276         else    quantity = ceil(quantity*1e6)*1/1e6;
00277         
00278 
00279         pair<map<ZAI, double>::iterator, bool> IResult;
00280         if(quantity > 0)
00281         {
00282                 IResult = fIsotopicQuantityNeeded.insert( pair<ZAI ,double>(zai, quantity));
00283                 if(IResult.second == false)
00284                         IResult.first->second += quantity;
00285         }
00286 DBGL;
00287         
00288 
00289 }
00290 
00291 //________________________________________________________________________
00292 void IsotopicVector::Need(const IsotopicVector& isotopicvector)
00293 {
00294 DBGL;
00295         map<ZAI ,double> isotopicquantity = isotopicvector.GetIsotopicQuantity();
00296         
00297         for(map<ZAI ,double >::iterator it = isotopicquantity.begin(); it != isotopicquantity.end(); it++)
00298                 Need( (*it).first, (*it).second);
00299 DBGL;
00300 }
00301 
00302 
00303 //________________________________________________________________________
00304 double  IsotopicVector::GetZAIIsotopicQuantity(const ZAI& zai) const
00305 {
00306 DBGL;
00307         map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
00308         
00309         map<ZAI ,double>::iterator it;
00310         it = IsotopicQuantity.find(zai);
00311         DBGL;
00312         
00313         if ( it != IsotopicQuantity.end() ) 
00314         {
00315                 return it->second;      
00316         }       
00317         else    
00318         {       
00319                 return 0;
00320         }
00321 }
00322 
00323 //________________________________________________________________________
00324 double  IsotopicVector::GetZAIIsotopicQuantity(const int z, const int a, const int i) const
00325 {
00326 DBGL;
00327         ZAI zai(z, a, i);
00328         return GetZAIIsotopicQuantity(zai);
00329 }
00330 
00331 IsotopicVector  IsotopicVector::GetAtomicComposition(int z) const
00332 {
00333 DBGL;
00334         IsotopicVector IV;
00335         map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity();
00336         map<ZAI ,double >::iterator it;
00337         for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
00338 
00339                 if( (*it).first.Z() == z  )  
00340                         IV += (*it).first * (*it).second;
00341 
00342         return IV;
00343 DBGL;
00344 }
00345 
00346 vector<int> IsotopicVector::GetAtomicSpecies() const
00347 {
00348 DBGL;
00349         vector<int> AtomicSpecies;
00350         
00351         map<ZAI ,double > IsotopicQuantity = GetIsotopicQuantity();
00352         map<ZAI ,double >::iterator it;
00353         for( it = IsotopicQuantity.begin(); it != IsotopicQuantity.end(); it++)
00354                 if( (int)AtomicSpecies.size() ==0 || (*it).first.Z() != AtomicSpecies.back() )  
00355                         AtomicSpecies.push_back((*it).first.Z());
00356 
00357 DBGL;
00358         return AtomicSpecies;
00359 }
00360 
00361 
00362 //________________________________________________________________________
00363 void IsotopicVector::Write(string filename, long int time) const 
00364 {
00365         ofstream IVfile(filename.c_str(), ios_base::app);               // Open the File
00366         if(!IVfile)
00367         cout << "!!Warning!! !!!IsotopicVector!!! \n Can't open \"" << filename << "\"\n" << endl;
00368 
00369         IVfile << time << " ";
00370         
00371         map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
00372         for(map<ZAI ,double >::iterator it = IsotopicQuantity.begin();
00373                                         it != IsotopicQuantity.end(); it++)
00374         {
00375                 IVfile << (*it).first.Z() << " ";
00376                 IVfile << (*it).first.A() << " ";
00377                 IVfile << (*it).first.I() << " ";
00378                 IVfile << (*it).second << " ";
00379                 
00380                 if((*it).first.Z()>70)
00381                 {
00382                         char Z[33];
00383                         sprintf(Z,"%d",(*it).first.Z());
00384                         char A[33];
00385                         sprintf(A,"%d",(*it).first.A());
00386                         char I[33];
00387                         sprintf(I,"%d",(*it).first.I());
00388                 
00389                 
00390                         string filenameZAI = filename+ "_DIR" +"/" + Z + "_" + A + "_"+ I;
00391                         string cmd = " mkdir " + filename+ "_DIR" +" 2>/dev/null";
00392                         int test = system(cmd.c_str());
00393                         if (test!=0) cout << "!!Warning!! !!!IsotopicVector!!! \n Can't create \"" << filename<<  "_DIR directory " << "\"\n" << endl;
00394                         ofstream IVfileZAI(filenameZAI.c_str(), ios_base::app);         // Open the File
00395 
00396                         if(!IVfileZAI)
00397                                 cout << "!!Warning!! !!!IsotopicVector!!! \n Can't open \"" << filenameZAI << "\"\n" << endl;
00398                 
00399                         IVfileZAI << time << " " << (*it).second << endl;;
00400                         IVfileZAI.close();              
00401                 }
00402         
00403         }
00404         IVfile << endl;
00405 }
00406 //________________________________________________________________________
00407 void IsotopicVector::Print(string option) const 
00408 {
00409 DBGL;
00410         cout << "**************************" << endl;
00411         cout << "*Isotopic Vector Property*" << endl;
00412         cout << "**************************" << endl << endl;
00413 
00414         bool QuantityPrint = false;
00415         bool DBPrint = false;
00416 
00417         QuantityPrint = true;
00418         
00419         if(QuantityPrint)
00420         {
00421                 cout << "*Isotopic Vector Quantity*" << endl;
00422                 map<ZAI ,double> IsotopicQuantity = GetIsotopicQuantity();
00423                 for(map<ZAI ,double >::iterator it = IsotopicQuantity.begin();
00424                                                 it != IsotopicQuantity.end(); it++)
00425                 {
00426                         cout << (*it).first.Z() << " ";
00427                         cout << (*it).first.A() << " ";
00428                         cout << (*it).first.I() << " ";
00429                         cout << ": " << (*it).second;
00430                         cout << endl;
00431                 }
00432                 cout << endl;
00433                 cout << "*Isotopic Vector Quantity Needed*" << endl;
00434                 map<ZAI ,double> IsotopicQuantityNeeded = GetIsotopicQuantityNeeded();
00435                 for(map<ZAI ,double >::iterator it = IsotopicQuantityNeeded.begin();
00436                                                 it != IsotopicQuantityNeeded.end(); it++)
00437                 {
00438                         cout << (*it).first.Z() << " ";
00439                         cout << (*it).first.A() << " ";
00440                         cout << (*it).first.I() << " ";
00441                         cout << ": " << (*it).second;
00442                         cout << endl;
00443                 }
00444                 cout << endl;
00445         }
00446         if(DBPrint)
00447         {
00448                 cout << "****Isotopic Vector DB****" << endl;
00449         }
00450 DBGL;
00451 }
00452 
00453 
00454 
00455 
00456 
00457 
 Tout Classes Fichiers Fonctions Variables Macros