Minor changes in libTFitPofB - the bigger problems have to be debugged later

src/external/TFitPofB-lib/classes/Makefile.libTFitPofB

@@ -37,14 +37,14 @@ endif
 
 # -- Linux
 ifeq ($(OS),LINUX)
-CXX           = g++-4.4.0
+CXX           = g++-4.4.2
 CXXFLAGS      = -O3 -fopenmp -Wall -Wno-trigraphs -fPIC
 PMUSRPATH     = ../../../include
 MNPATH        = $(ROOTSYS)/include
 LOCALPATH     = ../include
 FFTW3PATH     = /usr/include
 INCLUDES      = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FFTW3PATH)
-LD            = g++-4.4.0
+LD            = g++-4.4.2
 LDFLAGS       = -O
 SOFLAGS       = -shared -fopenmp
 SHLIB         = libTFitPofB.so

src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp

@@ -10,7 +10,7 @@
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2009 by Bastian M. Wojek                                *
+ *   Copyright (C) 2009 by Bastian M. Wojek,  Alexander Maisuradze         *
  *                                                                         *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -41,7 +41,7 @@ using namespace std;
 #define TWOPI 6.28318530717958647692
 
 const double fluxQuantum(2.067833667e7); // 10e14 times CGS units %% in CGS units should be 10^-7
-                                         // in this case this is Gauss per square nm
+                                         // in this case this is Gauss times square nm
 const double sqrt3(sqrt(3.0));
 
 double getXi(const double hc2) { // get xi given Hc2 in Gauss
@@ -82,11 +82,13 @@ TBulkVortexFieldCalc::~TBulkVortexFieldCalc() {
 
 TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const string& wisdom, const unsigned int steps) {
   fWisdom = wisdom;
-  fSteps = steps;
+  if (steps % 2) {
+    fSteps = steps + 1;
+  } else {
+    fSteps = steps;
+  }
   fParam.resize(3);
   fGridExists = false;
-  if (fSteps%2)
-    fSteps++;
 
   int init_threads(fftw_init_threads());
   if (init_threads)
@@ -195,22 +197,21 @@ TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const string& wisdo
 
 void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
   // SetParameters - method has to be called from the user before the calculation!!
-  double field(abs(fParam[0])), lambda(abs(fParam[1])), xi(abs(fParam[2]));
+  double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
   double Hc2(getHc2(xi));
 
-  double latConstTr(sqrt(fluxQuantum/field*sqrt(4.0/3.0)));
+  double latConstTr(2.0*sqrt(fluxQuantum/(field*sqrt3)));
   double xisq_2_scaled(2.0/3.0*pow(xi*PI/latConstTr,2.0)), lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0));
 
-  int NFFT(fSteps);
+  const int NFFT(fSteps);
-  int NFFT_2(fSteps/2);
+  const int NFFT_2(fSteps/2);
-  int NFFTsq(fSteps*fSteps);
+  const int NFFTsq(fSteps*fSteps);
 
    // fill the field Fourier components in the matrix
 
-  int m;
-
   // ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC
   if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) {
+  int m;
 #pragma omp parallel for default(shared) private(m) schedule(dynamic)
     for (m = 0; m < NFFTsq; m++) {
       fFFTout[m] = field;

src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp

@@ -97,7 +97,11 @@ void TFitPofBStartupHandler::OnEndDocument()
  */
 void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attributes)
 {
-  if (!strcmp(str, "data_path")) {
+  if (!strcmp(str, "LEM")) {
+    fKey = eLEM;
+  } else if (!strcmp(str, "VortexLattice")) {
+    fKey = eVortex;
+  } else if (!strcmp(str, "data_path")) {
     fKey = eDataPath;
   } else if (!strcmp(str, "energy")) {
     fKey = eEnergy;
@@ -111,8 +115,8 @@ void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attrib
     fKey = eNSteps;
   } else if (!strcmp(str, "N_VortexGrid")) {
     fKey = eGridSteps;
-  } else if (!strcmp(str, "N_VortexFourier")) {
+  } else if (!strcmp(str, "VortexSymmetry")) {
-    fKey = eVortexFourierComp;
+    fKey = eVortexSymmetry;
   }
 }
 
@@ -140,6 +144,12 @@ void TFitPofBStartupHandler::OnEndElement(const char *str)
 void TFitPofBStartupHandler::OnCharacters(const char *str)
 {
   switch (fKey) {
+    case eLEM:
+      fLEM = true;
+      break;
+    case eVortex:
+      fVortex = true;
+      break;
     case eDataPath:
       // set the data path to the given path
       fDataPath = str;
@@ -168,9 +178,12 @@ void TFitPofBStartupHandler::OnCharacters(const char *str)
       // convert str to int and assign it to the GridSteps-member
       fGridSteps = atoi(str);
       break;
-     case eVortexFourierComp:
+    case eVortexSymmetry:
-      // convert str to int and assign it to the VortexFourierComp-member
+      // convert str to int and assign it to the VortexSymmetry-member (square = 0, triangular = 1)
-      fVortexFourierComp = atoi(str);
+      if (!strcmp(str, "square"))
+        fVortexSymmetry = 2;
+      else
+        fVortexSymmetry = 1;
       break;
     default:
       break;
@@ -256,68 +269,92 @@ void TFitPofBStartupHandler::CheckLists()
 {
   // check if anything was set, and if not set some default stuff
 
-  // check if any data path is given
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check data path ...";
-  if (!fDataPath.size()) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: This is not going to work, you have to set a valid data path where to find the rge-files in the xml-file!" << endl;
-    exit(-1);
-  }
-
-  // check if any energies are given
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check energy list ..." << endl;
-  if (!fEnergyList.size()) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: Energy list empty! Setting the default list." << endl;
-    char eChar[5];
-    for(unsigned int i(0); i<33; i++) {
-      for(unsigned int j(0); j<10; j++) {
-        sprintf(eChar, "%02u_%u", i, j);
-        fEnergyList.push_back(string(eChar));
-      }
-    }
-  }
-
   // check if delta_t is given, if not set default
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check specified time resolution ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check specified time resolution ... " << endl;
   if(!fDeltat) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the time resolution. Setting the default." << endl;
+    cout << "TFitPofBStartupHandler::CheckLists: You did not specify the time resolution. Setting the default (10 ns)." << endl;
     fDeltat = 0.01;
+  } else {
+    cout << fDeltat << " us" << endl;
   }
 
   // check if delta_B is given, if not set default
   cout << endl << "TFitPofBStartupHandler::CheckLists: check specified field resolution ..." << endl;
   if(!fDeltaB) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the field resolution. Setting the default." << endl;
+    cout << "TFitPofBStartupHandler::CheckLists: You did not specify the field resolution. Setting the default (0.1 G)." << endl;
-    fDeltaB = 0.05;
+    fDeltaB = 0.1;
+  } else {
+    cout << fDeltaB << " G" << endl;
   }
 
   // check if any wisdom-file is specified
   cout << endl << "TFitPofBStartupHandler::CheckLists: check wisdom-file ..." << endl;
   if (!fWisdomFile.size()) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify a wisdom file. Setting the default." << endl;
+    cout << "TFitPofBStartupHandler::CheckLists: You did not specify a wisdom file. No FFTW plans will be loaded or saved." << endl;
-    fWisdomFile = "WordsOfWisdom.dat";
+    fWisdomFile = "";
+  } else {
+    cout << fWisdomFile << endl;
   }
 
-  // check if any number of steps for the theory function is specified
+
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for theory ..." << endl;
+  if (fLEM) {
-  if (!fNSteps) {
+
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the theory. Setting the default." << endl;
+    // check if any data path is given
-    fNSteps = 3000;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: check data path ..." << endl;
+    if (!fDataPath.size()) {
+      cout << "TFitPofBStartupHandler::CheckLists: This is not going to work, you have to set a valid data path where to find the rge-files in the xml-file!" << endl;
+      exit(-1);
+    } else {
+      cout << fDataPath << endl;
+    }
+
+    // check if any energies are given
+    cout << endl << "TFitPofBStartupHandler::CheckLists: check energy list ..." << endl;
+    if (!fEnergyList.size()) {
+      cout << "TFitPofBStartupHandler::CheckLists: Energy list empty! Setting the default list ( 0.0:0.1:32.9 keV)." << endl;
+      char eChar[5];
+      for(unsigned int i(0); i<33; i++) {
+        for(unsigned int j(0); j<10; j++) {
+          sprintf(eChar, "%02u_%u", i, j);
+          fEnergyList.push_back(string(eChar));
+        }
+      }
+    } else {
+      for (unsigned int i (0); i < fEnergyList.size(); i++)
+        cout << fEnergyList[i] << " ";
+      cout << endl;
+    }
+
+    // check if any number of steps for the theory function is specified
+    cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for theory ..." << endl;
+    if (!fNSteps) {
+      cout << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the theory. Setting the default (3000)." << endl;
+      fNSteps = 3000;
+    } else {
+      cout << fNSteps << endl;
+    }
   }
 
-  // check if any number of steps for the theory function is specified
+  if (fVortex) {
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl;
-  if (!fGridSteps) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the grid. Setting the default." << endl;
-    fGridSteps = 200;
-  }
 
-  // check if any number of steps for the theory function is specified
+    // check if any number of steps for the theory function is specified
-  cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex lattice Fourier components ..." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl;
-  if (!fVortexFourierComp) {
+    if (!fGridSteps) {
-    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of Vortex lattice Fourier components. Setting the default." << endl;
+      cout << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the grid. Setting the default (256)." << endl;
-    fVortexFourierComp = 17;
+      fGridSteps = 256;
-  }
+    } else {
+      cout << fGridSteps << endl;
+    }
 
+    // check if any number of steps for the theory function is specified
+    cout << endl << "TFitPofBStartupHandler::CheckLists: check symmetry of the Vortex lattice ..." << endl;
+    if (!fVortexSymmetry) {
+      cout << "TFitPofBStartupHandler::CheckLists: You did not specify the symmetry of the vortex lattice. Setting the default (triangular)." << endl;
+      fVortexSymmetry = 1;
+    } else {
+      cout << (fVortexSymmetry == 2 ? "square" : "triangular") << endl;
+    }
+  }
 }
 
 // end ---------------------------------------------------------------------

src/external/TFitPofB-lib/classes/TPofBCalc.cpp

@@ -391,7 +391,7 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   if (lastZerosStart >= fPBSize)
     lastZerosStart = fPBSize - 1;
 
-//  cout << endl << fBmin << " " << fBmax << " " << firstZerosEnd << " " << lastZerosStart << " " << numberOfPoints << endl;
+//  cout << endl << fBmin << " " << fBmax << " " << firstZerosEnd << " " << lastZerosStart << " " << numberOfSteps << endl;
 
   if (!vortexLattice->GridExists()) {
     vortexLattice->CalculateGrid();
@@ -401,7 +401,7 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   unsigned int fill_index, counter(0);
 
   for (unsigned int j(0); j < numberOfStepsSq_2; j++) {
-    fill_index = static_cast<unsigned int>(ceil(abs((vortexFields[j]/fDB))));
+    fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[j]/fDB))));
     if (fill_index >= fPBSize)
       fPB[fPBSize - 1] += 1.0;
     else

src/external/TFitPofB-lib/classes/TPofTCalc.cpp

@@ -80,7 +80,7 @@ TPofTCalc::TPofTCalc (const TPofBCalc *PofB, const string &wisdom, const vector<
 
   fTBin = 1.0/(gBar*double(fNFFT-1)*par[2]);
 
-  int NFFT_2p1(fNFFT/2 + 1);
+  const int NFFT_2p1(fNFFT/2 + 1);
 
   // allocating memory for the time- and polarisation vectors
 

src/external/TFitPofB-lib/classes/TVortex.cpp

@@ -133,7 +133,7 @@ double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) con
   if(fFirstCall){
     fPar = par;
 
-    for (unsigned int i(0); i < 3; i++){
+    for (unsigned int i(0); i < 3; i++) {
       fParForVortex[i] = fPar[i+1];
     }
     fFirstCall = false;

src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h

@@ -14,7 +14,7 @@
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007 by Andreas Suter, Bastian M. Wojek                 *
+ *   Copyright (C) 2007-2008 by Andreas Suter, Bastian M. Wojek            *
  *                                                                         *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -66,13 +66,15 @@ class TFitPofBStartupHandler : public TQObject {
     virtual const string GetWisdomFile() const { return fWisdomFile; }
     virtual const unsigned int GetNSteps() const { return fNSteps; }
     virtual const unsigned int GetGridSteps() const { return fGridSteps; }
-    virtual const unsigned int GetVortexFourierComp() const { return fVortexFourierComp; }
+    virtual const unsigned int GetVortexSymmetry() const { return fVortexSymmetry; }
 
   private:
-    enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps, eVortexFourierComp};
+    enum EKeyWords {eEmpty, eComment, eLEM, eVortex, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps, eVortexSymmetry};
 
     EKeyWords       fKey;
 
+    bool            fLEM;
+    bool            fVortex;
     string          fDataPath;
     vector<string>  fEnergyList;
     double          fDeltat;
@@ -80,7 +82,7 @@ class TFitPofBStartupHandler : public TQObject {
     string          fWisdomFile;
     unsigned int    fNSteps;
     unsigned int    fGridSteps;
-    unsigned int    fVortexFourierComp;
+    unsigned int    fVortexSymmetry;
 
   ClassDef(TFitPofBStartupHandler, 1)
 };

src/external/TFitPofB-lib/test/Makefile.testVortex-v2

@@ -6,12 +6,12 @@ ROOTLIBS      = $(shell $(ROOTSYS)/bin/root-config --libs)
 
 #---------------------------------------------------
 
-CXX           = g++-4.4.0
+CXX           = g++-4.4.2
 CXXFLAGS      = -O3 -Wall
 LOCALINCLUDE  = ../include
 ROOTINCLUDE   = $(ROOTSYS)/include
 INCLUDES      = -I$(LOCALINCLUDE) -I$(ROOTINCLUDE)
-LD            = g++-4.4.0
+LD            = g++-4.4.2
 LDFLAGS       = -O3 -L../classes -lTFitPofB -lfftw3_threads -lfftw3 -lm -lpthread -fopenmp -lPMusr
 
 # the output from the root-config script:

src/external/TFitPofB-lib/test/TFitPofB_startup.xml

@@ -2,39 +2,56 @@
 <TFitPofB>
     <comment>
         TFitPofB_startup.xml 
-        Defines path/prefix and energies (keV, format: %02u_%1u) of TrimSP-rge-files,
+        Here some external parameters used in conjunction with the TFitPofB-library are set.
-        path/name to the FFTW-wisdom-file and time/field binning (us/G)
+        There are different sections:
-        N_theory determines the number of points in "real space" where the theory function will be calculated
+        * Common parameters which are used for any kind of function
-        N_VortexGrid determines the number of points in which each axis of the "real space" will be divided when calculating the vortex B(x,y)
+        * Parameters for low energy uSR data analysis with given B(z)
-        N_VortexFourier defines the number of Fourier components which should be summed up in order to get B(x,y)
+        * Parameters for bulk uSR data analysis (currently only one model to calculate field distributions for a triangular vortex lattice)
-      </comment>
+
-    <data_path>/home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000-</data_path>
+        Common parameters:
+        - wisdom : sets the path to an FFTW-wisdom file - if there is no valid wisdom at the given path, wisdom handling will be disabled
+        - delta_t : time resolution of P(t) in microseconds
+        - delta_B : field resolution of P(B) in Gauss
+
+        LEM parameters (if the LEM tag is present it is expected that the data_path tag contains a valid path):
+        - data_path, energy_list : Defines path/prefix and energies (keV, format: %02u_%1u) of TrimSP-rge-files,
+        - N_theory : determines the number of points in "real space" where the theory function B(z) will be calculated (if necessary)
+
+        bulk parameters:
+        - N_VortexGrid : determines the number of points used for the calculation of the vortex lattice field distribution (the grid will be N*N)
+        - VortexSymmetry : specify the vortex lattice symmetry - either "square" or "triangular"
+    </comment>
     <wisdom>/home/l_wojek/analysis/WordsOfWisdom.dat</wisdom>
     <delta_t>0.01</delta_t>
     <delta_B>0.1</delta_B>
-    <N_theory>5000</N_theory>
+    <VortexLattice>
-    <N_VortexGrid>200</N_VortexGrid>
+        <N_VortexGrid>256</N_VortexGrid>
-    <N_VortexFourier>17</N_VortexFourier>
+        <VortexSymmetry>triangular</VortexSymmetry>
-    <energy_list>
+    </VortexLattice>
-        <energy>03_0</energy>
+    <LEM>
-        <energy>03_6</energy>
+        <data_path>/home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000-</data_path>
-        <energy>05_0</energy>
+        <N_theory>5000</N_theory>
-        <energy>05_3</energy>
+        <energy_list>
-        <energy>07_0</energy>
+          <energy>03_0</energy>
-        <energy>07_7</energy>
+          <energy>03_6</energy>
-        <energy>08_0</energy>
+          <energy>05_0</energy>
-        <energy>09_0</energy>
+          <energy>05_3</energy>
-        <energy>10_0</energy>
+          <energy>07_0</energy>
-        <energy>10_2</energy>
+          <energy>07_7</energy>
-        <energy>12_0</energy>
+          <energy>08_0</energy>
-        <energy>14_1</energy>
+          <energy>09_0</energy>
-        <energy>16_0</energy>
+          <energy>10_0</energy>
-        <energy>16_4</energy>    
+          <energy>10_2</energy>
-        <energy>18_0</energy>
+          <energy>12_0</energy>
-        <energy>19_7</energy>
+          <energy>14_1</energy>
-        <energy>20_0</energy>
+          <energy>16_0</energy>
-        <energy>22_0</energy>
+          <energy>16_4</energy>
-        <energy>24_0</energy>
+          <energy>18_0</energy>
-        <energy>24_6</energy>
+          <energy>19_7</energy>
+          <energy>20_0</energy>
+          <energy>22_0</energy>
+          <energy>24_0</energy>
+          <energy>24_6</energy>
      </energy_list>
+   </LEM>
 </TFitPofB>

src/external/TFitPofB-lib/test/testVortex-v2.cpp

@@ -9,7 +9,7 @@ using namespace std;
 
 int main(){
 
-  unsigned int NFFT(256);
+  unsigned int NFFT(512);
 
   vector<double> parForVortex;
   parForVortex.resize(3);
@@ -20,47 +20,47 @@ int main(){
 
   vector<double> parForPofB;
   parForPofB.push_back(0.01); //dt
-  parForPofB.push_back(0.5); //dB
+  parForPofB.push_back(2.0); //dB
 
   vector<double> parForPofT;
   parForPofT.push_back(0.0); //phase
   parForPofT.push_back(0.01); //dt
-  parForPofT.push_back(0.5); //dB
+  parForPofT.push_back(2.0); //dB
 
   TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc("/home/l_wojek/analysis/WordsOfWisdom.dat", NFFT);
 
-  parForVortex[0] = 500.0; //app.field
+  parForVortex[0] = 1000.0; //app.field
-  parForVortex[1] = 300.0; //lambda
+  parForVortex[1] = 1000.0; //lambda
-  parForVortex[2] = 3.0; //xi
+  parForVortex[2] = 4.0; //xi
 
   vortexLattice->SetParameters(parForVortex);
   vortexLattice->CalculateGrid();
 
-//   ofstream ofy("testVortex-B.dat");
+  ofstream ofy("testVortex-B.dat");
-//   for (unsigned int j(0); j < NFFT * NFFT; j++) {
+  for (unsigned int j(0); j < NFFT * NFFT; j++) {
-//     ofy << vortexLattice->DataB()[j] << " ";
+    ofy << vortexLattice->DataB()[j] << " ";
-//     if (!((j+1)%(NFFT)))
+    if (!((j+1)%(NFFT)))
-//       ofy << endl;
+      ofy << endl;
-//   }
+  }
-//   ofy.close();
+  ofy.close();
 
   TPofBCalc *PofB = new TPofBCalc(parForPofB);
   PofB->Calculate(vortexLattice, parForPofB);
 
-//   const double *b(PofB->DataB());
+  const double *b(PofB->DataB());
-//   const double *pb(PofB->DataPB());
+  double *pb(PofB->DataPB());
-//   unsigned int s(PofB->GetPBSize());
+  unsigned int s(PofB->GetPBSize());
-// 
+
     double test(0.0);
-// 
+
-//   ofstream ofx("testVortex.dat");
+  ofstream ofx("testVortex.dat");
-//   for (unsigned int i(0); i < s; i++) {
+  for (unsigned int i(0); i < s; i++) {
-//     ofx << b[i] << " " << pb[i] << endl;
+    ofx << b[i] << " " << pb[i] << endl;
-//     test+=pb[i];
+    test+=pb[i];
-//   }
+  }
-//   ofx.close();
+  ofx.close();
   
-//  cout << test << endl;
+ cout << test << endl;
 
   TPofTCalc poft(PofB, "/home/l_wojek/analysis/WordsOfWisdom.dat", parForPofT);
   
@@ -69,12 +69,12 @@ int main(){
   
   
 
-//     ofstream of8("testVortex-Pt.dat");
+     ofstream of8("testVortex-Pt.dat");
      for (double i(0.); i<12.0; i+=0.003) {
        test = poft.Eval(i);
-//       of8 << i << " " << poft.Eval(i) << endl;
+       of8 << i << " " << poft.Eval(i) << endl;
      }
-//     of8.close();
+     of8.close();
 
 //   parForVortex[0] = 500.0; //app.field
 //   parForVortex[1] = 100.0; //lambda