layer.cpp Source File

00001 #include "lattice.h"
00002 #include "gui.h"
00003 
00011 Layer::Layer() : TCanvas() {}
00012 
00013 Layer::Layer(MainFrame* main, Configuration* conf, const char* name, Int_t ww, Int_t wh, Int_t winid)
00014     : TCanvas(name, ww, wh, winid)
00015 {
00016     pLattice = 0;
00017     pEFormula = 0;
00018     pConfig = conf;
00019     pMainFrame = main;
00020     pInfoText = new TText();
00021     mcThread = 0;
00022     runSimulation = kFALSE;
00023     
00024     
00025     SetBorderMode(0);
00026 }
00027 
00028 Layer::~Layer()
00029 {
00030     delete pLattice;
00031 }
00032 
00033 void Layer::ExecuteEvent(Int_t event, Int_t px, Int_t py)
00034 {
00035     if (!gPad->IsEditable()) return;
00036     
00037     switch (event) {
00038         case kButton3Down:
00039             ;
00040         default:
00041             TCanvas::ExecuteEvent(event, px, py);
00042     }
00043 }
00044 
00045 void Layer::SetLattice(Int_t type, Int_t nx, Int_t ny, Double_t c)
00046 {
00047     if (pLattice) delete pLattice;
00048     
00049     cd();
00050     switch (type) {
00051         case kLSquare:
00052             pLattice = new SquareLattice(this, nx, ny, c);
00053         break;
00054         case kLHexagonal:
00055             pLattice = new HexLattice(this, nx, ny, c);
00056         break;
00057     }
00058     intLimit = GetX2() > GetY2() ? GetX2() : GetY2();
00059     pConfig->SetInteractionLimit(intLimit);
00060     Draw();
00061 }
00062 
00063 Bool_t Layer::MoveDefect(Double_t& deltaE)
00064 {
00065     Atom *orig, *neighbor;
00066     TObjArray *defects = pLattice->GetDefects();
00067     Int_t defCount = defects->GetLast()+1;
00068     Double_t oldPartEnergy;
00069     Int_t rnd;
00070     
00071     if ((defCount == pLattice->GetNAtoms()) || (defCount < 2)) return kFALSE;
00072     
00073     neighbor = 0;
00074     while (!neighbor) {
00075         rnd = (Int_t)(((Double_t)(defCount))*rand()/(RAND_MAX+1.0));
00076         orig = (*pLattice)[rnd];
00077         neighbor = pLattice->GetAtomNeighbor(orig);
00078     }
00079     
00080     oldPartEnergy = CompPartEnergy(orig);
00081     oldAtomIndex = (orig)->index;
00082     (orig)->RemoveDefect();
00083     newAtomIndex = neighbor->index;
00084     neighbor->InsertDefect();
00085     
00086     deltaE = CompPartEnergy(neighbor) - oldPartEnergy;
00087     return kTRUE;
00088 }
00089 
00090 void Layer::UndoMoveDefect()
00091 {
00092     Int_t tmpIndex;
00093     Atom** atoms;
00094     
00095     atoms = pLattice->GetAtoms();
00096     atoms[newAtomIndex]->RemoveDefect();
00097     atoms[oldAtomIndex]->InsertDefect();
00098     tmpIndex = oldAtomIndex;
00099     oldAtomIndex = newAtomIndex;
00100     newAtomIndex = tmpIndex;
00101 }
00102 
00103 Double_t Layer::CompEnergy()
00104 {
00105     Int_t dim;
00106     Int_t i, j, max, l, lmax;
00107     Double_t curXPos, curYPos;
00108     Double_t x; // First Variable in TEFormula, should be distance between defects
00109     Double_t y; // Second Variable, should be angle between defects
00110     Double_t dx, dy, lx, ly;
00111     Double_t energy;
00112     TObjArray *defects;
00113     const Atom *atom;
00114     
00115     defects = pLattice->GetDefects();
00116     
00117     energy = 0;
00118     max = defects->GetLast()+1;
00119     lx = GetX2();
00120     ly = GetY2();
00121     dim = pEFormula->GetNdim();
00122     
00123     for (i = 0; i < max; i++) {
00124         curXPos = (*pLattice)[i]->GetX();
00125         curYPos = (*pLattice)[i]->GetY();
00126         for (j = 0; j < max; j++) {
00127             atom = (*pLattice)[j];
00128             dx = atom->GetX() - curXPos;
00129             dy = atom->GetY() - curYPos;
00130             x = sqrt(dx*dx + dy*dy);
00131             if (x < intLimit) {
00132                 if (i != j) {
00133                     if (dim == 2) {
00134                         y = atan(dy/dx);
00135                         energy += pEFormula->Eval(x, y);
00136                     }
00137                     else if (dim == 1)
00138                         energy += pEFormula->Eval(x);
00139                 }
00140                 // ---- periodic boundary conditions! ----------------
00141                 lmax = 2;
00142                 while (x < intLimit) {
00143                     // one up
00144                     dy += ly;
00145                     x = sqrt(dx*dx + dy*dy);
00146                     if (dim == 2) {
00147                         y = atan(dy/dx);
00148                         energy += pEFormula->Eval(x, y);
00149                     }
00150                     else if (dim == 1)
00151                         energy += pEFormula->Eval(x);
00152                     // lmax/2 steps left
00153                     for (l = 0; l < lmax/2; l++) {
00154                         dx -= lx;
00155                         x = sqrt(dx*dx + dy*dy);
00156                         if (dim == 2) {
00157                             y = atan(dy/dx);
00158                             energy += pEFormula->Eval(x, y);
00159                         }
00160                         else if (dim == 1)
00161                             energy += pEFormula->Eval(x);
00162                     }
00163                     // lmax steps down
00164                     for (l = 0; l < lmax; l++) {
00165                         dy -= ly;
00166                         x = sqrt(dx*dx + dy*dy);
00167                         if (dim == 2) {
00168                             y = atan(dy/dx);
00169                             energy += pEFormula->Eval(x, y);
00170                         }
00171                         else if (dim == 1)
00172                             energy += pEFormula->Eval(x);
00173                     }
00174                     // lmax steps right
00175                     for (l = 0; l < lmax; l++) {
00176                         dx += lx;
00177                         x = sqrt(dx*dx + dy*dy);
00178                         if (dim == 2) {
00179                             y = atan(dy/dx);
00180                             energy += pEFormula->Eval(x, y);
00181                         }
00182                         else if (dim == 1)
00183                             energy += pEFormula->Eval(x);
00184                     }
00185                     // lmax steps up
00186                     for (l = 0; l < lmax; l++) {
00187                         dy += ly;
00188                         x = sqrt(dx*dx + dy*dy);
00189                         if (dim == 2) {
00190                             y = atan(dy/dx);
00191                             energy += pEFormula->Eval(x, y);
00192                         }
00193                         else if (dim == 1)
00194                             energy += pEFormula->Eval(x);
00195                     }
00196                     // lmax/2 steps left
00197                     for (l = 0; l < lmax/2 - 1; l++) {
00198                         dx -= lx;
00199                         x = sqrt(dx*dx + dy*dy);
00200                         if (dim == 2) {
00201                             y = atan(dy/dx);
00202                             energy += pEFormula->Eval(x, y);
00203                         }
00204                         else if (dim == 1)
00205                             energy += pEFormula->Eval(x);
00206                     }
00207                     dx -= lx;
00208                     x = sqrt(dx*dx + dy*dy);
00209                     if (dim == 2)
00210                         y = atan(dy/dx);
00211                     // circle completed
00212                     lmax += 2;
00213                 }
00214             }
00215         }   
00216     }
00217     
00218     curEnergy = energy / 2;
00219     return curEnergy;
00220 }
00221 
00222 Double_t Layer::CompPartEnergy(Atom* mainAtom)
00223 {
00224     Int_t dim;
00225     Int_t i, max, l, lmax;
00226     Double_t curXPos, curYPos;
00227     Double_t x; // First Variable in TEFormula, should be distance between defects
00228     Double_t y; // Second Variable, should be angle between defects
00229     Double_t dx, dy, lx, ly;
00230     Double_t partEnergy;
00231     TObjArray *defects;
00232     const Atom* atom;
00233     
00234     defects = pLattice->GetDefects();
00235     
00236     partEnergy = 0;
00237     max = defects->GetLast()+1;
00238     lx = GetX2();
00239     ly = GetY2();
00240     dim = pEFormula->GetNdim();
00241     curXPos = mainAtom->GetX();
00242     curYPos = mainAtom->GetY();
00243     for (i = 0; i < max; i++) {
00244         atom = (*pLattice)[i];
00245         dx = atom->GetX() - curXPos;
00246         dy = atom->GetY() - curYPos;
00247         x = sqrt(dx*dx + dy*dy);
00248         if (x < intLimit) {
00249             if (atom->index != mainAtom->index) {
00250                 if (dim == 2) {
00251                     y = atan(dy/dx);
00252                     partEnergy += pEFormula->Eval(x, y);
00253                 }
00254                 else if (dim == 1)
00255                     partEnergy += pEFormula->Eval(x);
00256             }
00257             // ---- periodic boundary conditions! ----------------
00258             lmax = 2;
00259             while (x < intLimit) {
00260                 // one up
00261                 dy += ly;
00262                 x = sqrt(dx*dx + dy*dy);
00263                 if (dim == 2) {
00264                     y = atan(dy/dx);
00265                     partEnergy += pEFormula->Eval(x, y);
00266                 }
00267                 else if (dim == 1)
00268                     partEnergy += pEFormula->Eval(x);
00269                 // lmax/2 steps left
00270                 for (l = 0; l < lmax/2; l++) {
00271                     dx -= lx;
00272                     x = sqrt(dx*dx + dy*dy);
00273                     if (dim == 2) {
00274                         y = atan(dy/dx);
00275                         partEnergy += pEFormula->Eval(x, y);
00276                     }
00277                     else if (dim == 1)
00278                         partEnergy += pEFormula->Eval(x);
00279                 }
00280                 // lmax steps down
00281                 for (l = 0; l < lmax; l++) {
00282                     dy -= ly;
00283                     x = sqrt(dx*dx + dy*dy);
00284                     if (dim == 2) {
00285                         y = atan(dy/dx);
00286                         partEnergy += pEFormula->Eval(x, y);
00287                     }
00288                     else if (dim == 1)
00289                         partEnergy += pEFormula->Eval(x);
00290                 }
00291                 // lmax steps right
00292                 for (l = 0; l < lmax; l++) {
00293                     dx += lx;
00294                     x = sqrt(dx*dx + dy*dy);
00295                     if (dim == 2) {
00296                         y = atan(dy/dx);
00297                         partEnergy += pEFormula->Eval(x, y);
00298                     }
00299                     else if (dim == 1)
00300                         partEnergy += pEFormula->Eval(x);
00301                 }
00302                 // lmax steps up
00303                 for (l = 0; l < lmax; l++) {
00304                     dy += ly;
00305                     x = sqrt(dx*dx + dy*dy);
00306                     if (dim == 2) {
00307                         y = atan(dy/dx);
00308                         partEnergy += pEFormula->Eval(x, y);
00309                     }
00310                     else if (dim == 1)
00311                         partEnergy += pEFormula->Eval(x);
00312                 }
00313                 // lmax/2 steps left
00314                 for (l = 0; l < lmax/2 - 1; l++) {
00315                     dx -= lx;
00316                     x = sqrt(dx*dx + dy*dy);
00317                     if (dim == 2) {
00318                         y = atan(dy/dx);
00319                         partEnergy += pEFormula->Eval(x, y);
00320                     }
00321                     else if (dim == 1)
00322                         partEnergy += pEFormula->Eval(x);
00323                 }
00324                 dx -= lx;
00325                 x = sqrt(dx*dx + dy*dy);
00326                 if (dim == 2)
00327                     y = atan(dy/dx);
00328                 // circle completed
00329                 lmax += 2;
00330             }
00331         }
00332     }
00333     
00334     return partEnergy;
00335 }
00336 
00337 void Layer::MonteCarlo(void* arg)
00338 {
00339     ULong_t moves;
00340     Int_t imgCount, bMoved;
00341     Double_t prob, rnd, b, deltaE;
00342     char* imgFilename;
00343     const char *str;
00344     MainFrame* pMainFrame = (MainFrame*)arg;
00345     Layer* layer = pMainFrame->GetLayer();
00346     Configuration* pConfig = layer->GetConfig();
00347     
00348     imgCount = 0;
00349     moves = 0;
00350     bMoved = 0;
00351     
00352     TThread::Lock();
00353     pMainFrame->SetInfoStatus();
00354     layer->CompEnergy();
00355     
00356     pMainFrame->SetInfoEnergy(layer->GetCurEnergy());
00357     
00358     if (pConfig->AutoSave()) {
00359         str = pConfig->GetImagePrefix();
00360         imgFilename = new char[strlen(str) + 10];
00361         switch (pConfig->GetImageFormat()) {
00362             case 1:
00363                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.ps", str, imgCount);
00364                 break;
00365             case 2:
00366                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.eps", str, imgCount);
00367                 break;
00368             case 3:
00369                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.gif", str, imgCount);
00370                 break;
00371             case 4:
00372                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.C", str, imgCount);
00373                 break;
00374             case 5:
00375                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.cxx", str, imgCount);
00376                 break;
00377             case 6:
00378                 snprintf(imgFilename, strlen(str) + 10, "%s%04d.root", str, imgCount);
00379                 break;
00380         }
00381         
00382         gPad->SaveAs(gSystem->ConcatFileName(pConfig->GetImagePath(), imgFilename));
00383         delete imgFilename;
00384         imgCount++;
00385         if (pMainFrame->IsInfoEnabled());
00386             pMainFrame->SetInfoImages(imgCount);
00387     }
00388     TThread::UnLock();
00389     
00390     while (layer->RunSimulation()) {
00391         TThread::Lock();
00392         layer->GetLattice()->SetInfo(kFALSE);
00393         if (layer->MoveDefect(deltaE)) { // moving defect
00394             if (deltaE > 0) {
00395                 b = 1/(layer->GetBoltzConst()*layer->GetTemp());
00396                 prob = exp(-b*deltaE);
00397                 rnd = ((double)rand())/((double)RAND_MAX);
00398                 if (prob < rnd) {   // boltzmann
00399                     layer->UndoMoveDefect();
00400                     bMoved = 0;
00401                 }
00402                 else {
00403                     layer->SetCurEnergy(layer->GetCurEnergy() + deltaE);
00404                     moves++;
00405                     bMoved = 1;
00406                     if (pMainFrame->IsInfoEnabled()) {
00407                         pMainFrame->SetInfoMoves(moves);
00408                         pMainFrame->SetInfoStatus();
00409                     }
00410                 }
00411             }
00412             else {
00413                 layer->SetCurEnergy(layer->GetCurEnergy() + deltaE);
00414                 moves++;
00415                 bMoved = 1;
00416                 if (pMainFrame->IsInfoEnabled()) {
00417                     pMainFrame->SetInfoMoves(moves);
00418                     pMainFrame->SetInfoStatus();
00419                 }
00420             }
00421         }
00422         else bMoved = 0;
00423         layer->GetLattice()->SetInfo(kTRUE);
00424         TThread::UnLock();
00425         
00426         if (pConfig->UpdateCanvas() && bMoved)
00427             if ((moves % pConfig->GetUpdateInterval()) == 0) {
00428                 if (pMainFrame->IsInfoEnabled())
00429                     pMainFrame->SetInfoEnergy(layer->GetCurEnergy());
00430                 gPad->Modified();
00431                 gPad->Update();         }
00432             
00433         if (pConfig->AutoSave() && bMoved)
00434             if ((moves % pConfig->GetSaveInterval()) == 0) {
00435                 TThread::Lock();
00436                 str = pConfig->GetImagePrefix();
00437                 imgFilename = new char[strlen(str) + 10];
00438                 switch (pConfig->GetImageFormat()) {
00439                     case 1:
00440                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.ps", str, imgCount);
00441                         break;
00442                     case 2:
00443                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.eps", str, imgCount);
00444                         break;
00445                     case 3:
00446                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.gif", str, imgCount);
00447                         break;
00448                     case 4:
00449                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.C", str, imgCount);
00450                         break;
00451                     case 5:
00452                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.cxx", str, imgCount);
00453                         break;
00454                     case 6:
00455                         snprintf(imgFilename, strlen(str) + 10, "%s%04d.root", str, imgCount);
00456                         break;
00457                 }
00458                 
00459                 gPad->SaveAs(gSystem->ConcatFileName(pConfig->GetImagePath(), imgFilename));
00460                 delete imgFilename;
00461                 imgCount++;
00462                 if (pMainFrame->IsInfoEnabled());
00463                     pMainFrame->SetInfoImages(imgCount);
00464                 TThread::UnLock();
00465             }
00466     }
00467     
00468     if (pMainFrame->IsInfoEnabled())
00469         pMainFrame->SetInfoStatus();
00470         
00471     TThread::Exit();
00472 }
00473 
00474 Int_t Layer::MonteCarloStart(MainFrame *main)
00475 {
00476     runSimulation = kTRUE;
00477     
00478     if (!mcThread) {
00479         mcThread = new TThread("MonteCarlo", (void(*) (void*))&MonteCarlo, (void*)main);
00480         mcThread->Run();
00481         return 0;
00482     }
00483     
00484     return 1;
00485 }
00486 
00487 Int_t Layer::MonteCarloStop()
00488 {
00489     char* title;
00490     Long_t id;
00491     Int_t size;
00492     
00493     size = strlen(mcThread->GetName())+1;
00494     id = mcThread->GetId();
00495     title = new char[size];
00496     strncpy(title, mcThread->GetName(), size);
00497     
00498     runSimulation = kFALSE;
00499     TThread::Join(mcThread->GetId());
00500     
00501     if (mcThread) {
00502         delete mcThread;
00503         mcThread = 0;
00504     }
00505     
00506     cout << " Thread " << title << "." << id << " stopped" << endl << endl;
00507     delete title;
00508     return 0;
00509 }