diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp
index 78ab4e4aeabd81d0a48138ff43c378abc6703f9c..2f1bebed7e11e410391831aee73164668fd3eb05 100644
--- a/App/src/IsGISAXSTools.cpp
+++ b/App/src/IsGISAXSTools.cpp
@@ -141,7 +141,7 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
     }
 
     // creation of 2D with variable bin size
-//    std::cout << "XXX " << (int)haxises[0].nbins << " "  << (int)haxises[1].nbins;
+//    std::cout << "XXXYYY " << (int)haxises[0].nbins << " "  << (int)haxises[1].nbins;
 //    for(size_t i=0; i<haxises[0].xbins.size(); ++i) {
 //        std::cout << i << " axis0:" << haxises[0].xbins[i] << std::endl;
 //    }
@@ -172,6 +172,77 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s
 
 
 
+
+/* ************************************************************************* */
+// create TH2D from OutputData
+/* ************************************************************************* */
+//TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const std::string &histo_name)
+//{
+//    assert(&output);
+//    if (output.getNdimensions() !=2) throw( "IsGISAXSTools::getOutputDataTH2D() -> Warning! Expected number of dimensiobs is 2.");
+
+//    std::vector<AxisStructure > haxises;
+//    haxises.resize(output.getNdimensions());
+
+//    // we assume variable bin size and prepare [nbins+1] array of left edges of each bin plus right edge of the last bin
+//    for(size_t i_axis=0; i_axis<output.getNdimensions(); ++i_axis) {
+//        const IAxis *axis = output.getAxis(i_axis);
+//        if( !axis ) throw("IsGISAXSTools::getOutputDataTH123D() -> Error! Can't cast axis");
+//        double dx(0);
+//        haxises[i_axis].nbins = axis->getSize();
+//        haxises[i_axis].name = axis->getName();
+//        if( axis->getSize() == 0) {
+//            throw LogicErrorException("IsGISAXSTools::getOutputDataTH123D() -> Error! Axis with zero size.");
+//        }else if( axis->getSize() == 1 ) {
+//            // only one bin, let's invent fake bin size
+//            dx = 0.1*(*axis)[0];
+//            haxises[i_axis].xbins.push_back((*axis)[0]-dx/2.);
+//            haxises[i_axis].xbins.push_back((*axis)[0]+dx/2.);
+//        }else {
+//            for(size_t i_bin=0; i_bin<axis->getSize(); ++i_bin) {
+//                if(i_bin == 0) {
+//                    dx = (*axis)[1]-(*axis)[0];
+//                } else {
+//                    dx = (*axis)[i_bin] - (*axis)[i_bin-1];
+//                }
+//                haxises[i_axis].xbins.push_back( (*axis)[i_bin] - dx/2.);
+//            }
+//        haxises[i_axis].xbins.push_back((*axis)[axis->getSize()-1] + dx/2.); // right bin edge of last bin, so for 100 bins size of vector will be 101
+//        }
+//    }
+
+//    // creation of 2D with variable bin size
+//    std::cout << "XXXYYY " << (int)haxises[0].nbins << " "  << (int)haxises[1].nbins;
+//    for(size_t i=0; i<haxises[0].xbins.size(); ++i) {
+//        std::cout << i << " axis0:" << haxises[0].xbins[i] << std::endl;
+//    }
+//    for(size_t i=0; i<haxises[1].xbins.size(); ++i) {
+//        std::cout << i << " axis1:" << haxises[1].xbins[i] << std::endl;
+//    }
+
+//    TH2D *hist2 = new TH2D(histo_name.c_str(), histo_name.c_str(), (int)haxises[0].nbins, &haxises[0].xbins[0], (int)haxises[1].nbins, &haxises[1].xbins[0]);
+//    hist2->GetXaxis()->SetTitle( haxises[0].name.c_str() );
+//    hist2->GetYaxis()->SetTitle( haxises[1].name.c_str() );
+
+//    OutputData<double>::const_iterator it = output.begin();
+//    while (it != output.end())
+//    {
+//        double x = output.getValueOfAxis( haxises[0].name.c_str(), it.getIndex() );
+//        double y = output.getValueOfAxis( haxises[1].name.c_str(), it.getIndex() );
+//        double value = *it++;
+//        //std::cout << "OOO " << x << " " << y << std::endl;
+//        hist2->Fill(x, y, value);
+//    }
+//    hist2->SetContour(50);
+//    hist2->SetStats(0);
+//    hist2->GetYaxis()->SetTitleOffset(1.1);
+
+//    gStyle->SetPalette(1);
+//    gStyle->SetOptStat(0);
+//    return hist2;
+//}
+
+
 /* ************************************************************************* */
 // create TH1D, TH2D or TH3D from OutputData
 /* ************************************************************************* */
diff --git a/App/src/TestMesoCrystal2.cpp b/App/src/TestMesoCrystal2.cpp
index f51bcc52b899d266032075106bf0d49333991578..8c87b1deeadf9db6f60af76e22af20a1beb47dad 100644
--- a/App/src/TestMesoCrystal2.cpp
+++ b/App/src/TestMesoCrystal2.cpp
@@ -107,7 +107,7 @@ void TestMesoCrystal2::draw_results()
     m_simulation->runSimulation();
     m_simulation->normalize();
 
-    IsGISAXSTools::drawOutputDataComparisonResults(*m_simulation->getOutputData(), *m_real_data, "found", "founf params", 100, 1e6, 100);
+    IsGISAXSTools::drawOutputDataComparisonResults(*m_simulation->getOutputData(), *m_real_data, "found", "found params", 100, 1e6, 100);
 
     TCanvas *c1 = new TCanvas("meso_real_data","meso_real_data",1024, 768);
     c1->cd(); gPad->SetLogz();  gPad->SetRightMargin(0.12); gPad->SetLeftMargin(0.125);
@@ -125,6 +125,11 @@ void TestMesoCrystal2::draw_results()
     hist_simu->GetYaxis()->SetTitleOffset(1.35);
     hist_simu->DrawCopy("CONT4 Z");
 
+//    OutputDataIOFactory::writeOutputData(*m_simulation->getOutputData(), "meso_simul.txt");
+//    OutputDataIOFactory::writeOutputData(*m_real_data, "meso_real.txt");
+//    OutputData<double> *tmp_real = OutputDataIOFactory::getOutputData("meso_real.txt");
+//    OutputData<double> *tmp_simul = OutputDataIOFactory::getOutputData("meso_simul.txt");
+//    IsGISAXSTools::drawOutputDataComparisonResults(*tmp_real, *tmp_simul, "tmp", "bla bla", 100, 1e6);
 }
 
 
@@ -174,7 +179,6 @@ void TestMesoCrystal2::run_fit()
 void TestMesoCrystal2::initializeRealData()
 {
     delete m_real_data;
-    // reading data file
     //std::string file_name = Utils::FileSystem::GetHomePath()+"Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_qyqz.txt.gz";
     std::string file_name = Utils::FileSystem::GetHomePath()+"Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_phitheta.txt.gz";
 
@@ -496,7 +500,7 @@ void TestMesoCrystal2::initializeSimulation(const OutputData<double> *output_dat
         // if there is output_data as input parameter, build detector using output_data axises
         const IAxis *axis0 = output_data->getAxis(0);
         const IAxis *axis1 = output_data->getAxis(1);
-        std::cout << "Axis!!! " << axis0->getSize() << " " << axis0->getMin() << " " << axis0->getMax() << " " << axis1->getSize() << " " << axis1->getMin() << " " << axis1->getMax() << std::endl;
+        //std::cout << "Axis!!! " << axis0->getSize() << " " << axis0->getMin() << " " << axis0->getMax() << " " << axis1->getSize() << " " << axis1->getMin() << " " << axis1->getMax() << std::endl;
         //m_simulation->setDetectorParameters(axis0->getSize(), axis0->getMin(), axis0->getMax(), axis1->getSize(), axis1->getMin(), axis1->getMax());
         m_simulation->setDetectorParameters(*m_real_data);
     }
diff --git a/Core/Tools/inc/OutputDataWriteStrategy.h b/Core/Tools/inc/OutputDataWriteStrategy.h
index d0750e0cb1767a60395ace2276ee902d382101be..d1998b9167e828feb98efd48e898df1bc5d9b588 100644
--- a/Core/Tools/inc/OutputDataWriteStrategy.h
+++ b/Core/Tools/inc/OutputDataWriteStrategy.h
@@ -41,4 +41,15 @@ public:
     virtual void writeOutputData(const OutputData<double> &data, std::ostream &output_stream);
 };
 
+
+//! Strategy to write OutputData to ascii files
+//! 1d array for x-axis, 1d array for y-axis, 2d array for data
+
+class OutputDataWriteStreamV1 : public IOutputDataWriteStrategy
+{
+public:
+    virtual void writeOutputData(const OutputData<double> &data, std::ostream &output_stream);
+};
+
+
 #endif // OUTPUTDATAWRITESTRATEGY_H
diff --git a/Core/Tools/src/OutputDataIOFactory.cpp b/Core/Tools/src/OutputDataIOFactory.cpp
index 6d0235613c313432b9be13092b14e1c4af18d4b2..0afada18e75fcfa885b88cbeb2f00ce359d2feed 100644
--- a/Core/Tools/src/OutputDataIOFactory.cpp
+++ b/Core/Tools/src/OutputDataIOFactory.cpp
@@ -63,6 +63,8 @@ OutputDataIOFactory::OutputDataWriter_t OutputDataIOFactory::getWriter(const std
     IOutputDataWriteStrategy *write_strategy(0);
     if( Utils::FileSystem::GetFileExtension(file_name) == ".ima") {
         write_strategy = new OutputDataWriteStreamIMA();
+    }else if(Utils::FileSystem::GetFileExtension(file_name) == ".txt") {
+        write_strategy = new OutputDataWriteStreamV1();
     } else {
         throw LogicErrorException("OutputDataIOFactory::getWriter() -> Error. Don't know how to write file '" + file_name+std::string("'"));
     }
diff --git a/Core/Tools/src/OutputDataReadStrategy.cpp b/Core/Tools/src/OutputDataReadStrategy.cpp
index a134df05d6ddc9448117eb3015063c668032b7f6..88099df4af45784d0bef8b1a29a9d76f9296ba45 100644
--- a/Core/Tools/src/OutputDataReadStrategy.cpp
+++ b/Core/Tools/src/OutputDataReadStrategy.cpp
@@ -69,8 +69,10 @@ OutputData<double > *OutputDataReadStreamIMA::readOutputData(std::istream &input
     }
 
     // creating new OutputData and filling it with values from buffer_2d
-    int y_size = (int)buff_2d.size();
-    int x_size = buff_2d.size() ? (int)buff_2d[0].size() : 0;
+//    int y_size = (int)buff_2d.size();
+//    int x_size = buff_2d.size() ? (int)buff_2d[0].size() : 0;
+    int x_size = (int)buff_2d.size();
+    int y_size = buff_2d.size() ? (int)buff_2d[0].size() : 0;
     OutputData<double> *p_result = new OutputData<double>;
     p_result->addAxis(NDetector2d::PHI_AXIS_NAME, x_size, 0.0, double(x_size));
     p_result->addAxis(NDetector2d::ALPHA_AXIS_NAME, y_size, 0.0, double(y_size));
@@ -102,7 +104,7 @@ OutputData<double > *OutputDataReadStreamV1::readOutputData(std::istream &input_
 {
     std::string sline;
     vdouble1d_t buff_xaxis, buff_yaxis;
-    vdouble2d_t buff_data; // [y][x]
+    vdouble2d_t buff_data; // [x][y]
 
     while( std::getline(input_stream, sline) )
     {
@@ -119,14 +121,14 @@ OutputData<double > *OutputDataReadStreamV1::readOutputData(std::istream &input_
         }
     }
 
-    // check consistency of y dimension and data buffer
-    if( buff_data.size() != buff_yaxis.size()) {
-        throw LogicErrorException("OutputDataReadASCII::readOutputData() -> Error. Unconsistent y-size.");
-    }
     // check consistency of x dimension and data buffer
-    for(size_t i = 0; i<buff_yaxis.size(); ++i) {
-        if( buff_data[i].size() != buff_xaxis.size()) {
-            throw LogicErrorException("OutputDataReadASCII::readOutputData() -> Error. Unconsistent x-size.");
+    if( buff_data.size() != buff_xaxis.size()) {
+        throw LogicErrorException("OutputDataReadASCII::readOutputData() -> Error. Unconsistent x-size.");
+    }
+    // check consistency of y dimension and data buffer
+    for(size_t i = 0; i<buff_xaxis.size(); ++i) {
+        if( buff_data[i].size() != buff_yaxis.size()) {
+            throw LogicErrorException("OutputDataReadASCII::readOutputData() -> Error. Unconsistent y-size.");
         }
     }
 
@@ -145,7 +147,7 @@ OutputData<double > *OutputDataReadStreamV1::readOutputData(std::istream &input_
     {
         size_t index_x = p_result->toCoordinates(it.getIndex())[0];
         size_t index_y = p_result->toCoordinates(it.getIndex())[1];
-        *it = buff_data[index_y][index_x];
+        *it = buff_data[index_x][index_y];
         ++it;
     }
     return p_result;
diff --git a/Core/Tools/src/OutputDataWriteStrategy.cpp b/Core/Tools/src/OutputDataWriteStrategy.cpp
index 88b4cf1602b6b8932b35883bc7847eef0d7a04da..f83245cf2823f08ba1bf4739f162c08b29f173e5 100644
--- a/Core/Tools/src/OutputDataWriteStrategy.cpp
+++ b/Core/Tools/src/OutputDataWriteStrategy.cpp
@@ -18,8 +18,9 @@
 #include <fstream>
 #include <sstream>
 #include <iomanip>
+#include <cassert>
 
-//! read data from ASCII file (2D assumed) and fill newly created OutputData with it
+// write OutputData to IsGisaxs *.ima files
 
 void OutputDataWriteStreamIMA::writeOutputData(const OutputData<double> &data, std::ostream &output_stream)
 {
@@ -31,3 +32,41 @@ void OutputDataWriteStreamIMA::writeOutputData(const OutputData<double> &data, s
         if(it.getIndex()%row_length==0) output_stream << std::endl;
     }
 }
+
+
+// write OutputData to ascii files
+// '#' - comments (not treated)
+// ' ' - delimeter between double's
+// Each line is 1D array (or 1D slice of 2D array)
+// Records:
+// line representing x bins [nX]
+// line representing y bins [nY]
+// [nX] lines of [nY] size representing data itself
+void OutputDataWriteStreamV1::writeOutputData(const OutputData<double> &data, std::ostream &output_stream)
+{
+    assert(data.getNdimensions() == 2);
+
+    output_stream << "# 2D scattering data" << std::endl;
+    output_stream << "# shape " << data.getAxis(0)->getSize() << " " << data.getAxis(1)->getSize() << std::endl;
+
+    // writing x and y axis
+
+    for(size_t i_dim=0; i_dim<data.getNdimensions(); ++i_dim) {
+        (i_dim == 0 ? output_stream << "# 0-axis [phi]" : output_stream << "# 1-axis [theta]");
+        output_stream << std::endl;
+        for(size_t i_bin=0; i_bin<data.getAxis(i_dim)->getSize(); ++i_bin) {
+            output_stream << std::scientific << std::setprecision(m_precision) << (*data.getAxis(i_dim))[i_bin] << "    ";
+        }
+        output_stream << std::endl;
+    }
+
+    output_stream << "# data" << std::endl;
+    size_t row_length = data.getAxis(1)->getSize();
+    OutputData<double>::const_iterator it = data.begin();
+    while(it != data.end()) {
+        double z_value = *it++;
+        output_stream << std::scientific << std::setprecision(m_precision) << z_value << "    ";
+        if(it.getIndex()%row_length==0) output_stream << std::endl;
+    }
+}
+
diff --git a/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_phitheta.txt.gz b/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_phitheta.txt.gz
index d958ba767f85b064f875ee94f5b36d8e88421bc2..6e64191372bd072ace7c6746fb565b4d6ec976c5 100644
Binary files a/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_phitheta.txt.gz and b/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_phitheta.txt.gz differ
diff --git a/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_qyqz.txt.gz b/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_qyqz.txt.gz
deleted file mode 100644
index d8863f1d2a718141af2c7037bf5a55771087c8dc..0000000000000000000000000000000000000000
Binary files a/Examples/MesoCrystals/ex02_fitspheres/004_230_P144_im_full_qyqz.txt.gz and /dev/null differ
diff --git a/Examples/MesoCrystals/ex02_fitspheres/edf2ascii.py b/Examples/MesoCrystals/ex02_fitspheres/edf2ascii.py
index f2544004ceae6fd0ce26736e01a47295074425ad..82b982a76ee70d1abc1472e5008a052c04a52042 100644
--- a/Examples/MesoCrystals/ex02_fitspheres/edf2ascii.py
+++ b/Examples/MesoCrystals/ex02_fitspheres/edf2ascii.py
@@ -20,24 +20,61 @@ CENTER_Y = 942.0
 # 1D array with y-values (position of rows) as single line
 # 2D data array where single line represent second dimension [nrows][ncols]
 #-----------------------------------------------------------------------------
-def save_text_file(filename, cols, rows, data):
+#def save_text_file(filename, cols, rows, data):
+  #if ".gz" in filename:
+    #f = gzip.open(filename, 'wb')
+  #else:
+    #f = file(filename, 'w')
+  #f.write("# 2D scattering data\n")
+  ## writing header containing the shape of array
+  #header = "# shape"
+  #for ndim in range(0,len(data.shape)):
+      #header += " " + str(data.shape[ndim])
+  #f.write(header+"\n")
+  ## writing x-axis values
+  #f.write("# x-axis \n")
+  #savetxt(f, cols, fmt='%1.10e', delimiter=" ", newline=" ")
+  #f.write("\n")
+  ## writing y-axis values
+  #f.write("# y-axis \n")
+  #savetxt(f, rows, fmt='%1.10e', delimiter=" ", newline=" ")
+  #f.write("\n")
+
+  ## writing data array
+  #f.write("# data\n")
+  #savetxt(f, data, fmt='%1.10e', delimiter=" ", newline="\n")
+  #f.close()
+
+
+#-----------------------------------------------------------------------------
+# here we've got data[theta][phi]
+# We have to invert to BornAgain convention [phi][theta]
+# saving 3 numpy arrays in text file
+# 1D array with x-values (position of columns) as single line
+# 1D array with y-values (position of rows) as single line
+# 2D data array where single line represent second dimension [nrows][ncols]
+#-----------------------------------------------------------------------------
+def save_text_file2(filename, phi, theta, data):
   if ".gz" in filename:
     f = gzip.open(filename, 'wb')
   else:
     f = file(filename, 'w')
   f.write("# 2D scattering data\n")
   # writing header containing the shape of array
+  # !!!
+  data = data.transpose()
+
   header = "# shape"
   for ndim in range(0,len(data.shape)):
       header += " " + str(data.shape[ndim])
   f.write(header+"\n")
   # writing x-axis values
-  f.write("# x-axis \n")
-  savetxt(f, cols, fmt='%1.10e', delimiter=" ", newline=" ")
+  f.write("# 0-axis [phi] \n")
+  savetxt(f, phi, fmt='%1.10e', delimiter=" ", newline=" ")
   f.write("\n")
   # writing y-axis values
-  f.write("# y-axis \n")
-  savetxt(f, rows, fmt='%1.10e', delimiter=" ", newline=" ")
+  f.write("# 1-axis [theta] \n")
+  savetxt(f, theta, fmt='%1.10e', delimiter=" ", newline=" ")
   f.write("\n")
 
   # writing data array
@@ -46,6 +83,8 @@ def save_text_file(filename, cols, rows, data):
   f.close()
 
 
+
+
 #-----------------------------------------------------------------------------
 # modifying 3 arrays (rows, cols and data) to fit in the range rmin, rmax, cmin, cmax
 # input:
@@ -144,6 +183,11 @@ def make_file(filename):
   Q_z=Q_z[ ::-1]
   data=data[ ::-1,:]
 
+  # The grid orientation follows the MATLAB convention: an array C with shape (nrows, ncolumns) is plotted 
+  # with the column number as X and the row number as Y, increasing up; hence it is plotted the way the 
+  # array would be printed, except that the Y axis is reversed.
+  # e.g. here we have data[Qz][Qy] (data[theta][phi])
+
   # ---
   # to save png file
   gca().set_aspect('equal')
@@ -160,14 +204,25 @@ def make_file(filename):
 
   # data .vs. qy,qz
   newqy, newqz, newdata = modify_arrays(Q_y, Q_z, data, 0.0, 0.3, 0.0, 0.3)
-  save_text_file(filename.rsplit('.',1)[0]+'_qyqz.txt.gz', newqy, newqz, newdata)
+  #save_text_file(filename.rsplit('.',1)[0]+'_qyqz.txt.gz', newqy, newqz, newdata)
 
   # data .vs. phi,theta
   phi = phi/180.*pi
   tth=(tth-0.32)/180.*pi
   tth=tth[ ::-1]
   newphi, newtheta, newdata2 = modify_arrays(phi, tth, data, 0.02, 0.06, 0.0, 0.04)
-  save_text_file(filename.rsplit('.',1)[0]+'_phitheta.txt.gz', newphi, newtheta, newdata2)
+  #save_text_file(filename.rsplit('.',1)[0]+'_phitheta.txt.gz', newphi, newtheta, newdata2)
+  
+  gca().set_aspect('equal')
+  pcolormesh(newphi,newtheta, newdata2, norm=norm, cmap='gist_ncar')
+  xlabel('$Q_y [\AA^{-1}]$ (phi)')
+  ylabel('$Q_z [\AA^{-1}]$ (theta)')
+  colorbar()
+  title(file_title)
+  savefig('tmp.png',dpi=300)
+  close()
+  
+  save_text_file2(filename.rsplit('.',1)[0]+'_phitheta.txt.gz', newphi, newtheta, newdata2)
 
 
 #-----------------------------------------------------------------------------
diff --git a/Tests/FunctionalTests/TestPyCore/README b/Tests/FunctionalTests/TestPyCore/README
index 4c6ed7d646002e9c44c9ae1872e04ac5a1ff63ee..eb6f0922834274f439283b3a367aebf996ea4243 100644
--- a/Tests/FunctionalTests/TestPyCore/README
+++ b/Tests/FunctionalTests/TestPyCore/README
@@ -1,6 +1,6 @@
 Collection of functional tests (Python)
 
-Collection contains set of functional tests for BornAgainCore library.
+Collection contains functional tests for BornAgainCore library.
 Each test defines simple geometry, runs simulation and then compare results of the simulation with reference data.
 
 To build all tests