/**********************************************************************************
 * Copyright 2010 Christoph Pflaum 
 * 		Department Informatik Lehrstuhl 10 - Systemsimulation
 *		Friedrich-Alexander Universität Erlangen-Nürnberg
 * 
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 **********************************************************************************/



#include <string>
#include <stdexcept>
#include "../../../program/source/mympi.h"  // von 3D UGBlocks
#include "../abbrevi2D.h"
//#include "../../../linAlgExptemp/source/linAlg.h"
//#include "../parameter.h"
#include "../math_lib/math_lib.h"
#include "../grid/elements2D.h"
#include "../grid/parti2D.h"
#include "../grid/ug2D.h"
#include "../grid/blockgrid2D.h"
#include "../grid/marker2D.h"
//#include "parallel.h"
#include "extemp2D.h"
#include "variable2D.h"
#include "variable2D_cc.h"
#include "update2D.h"
#include "interpolTwoD.h"
#include "../extemp/co_fu2D.h"

#include "assert.h"

/////////////////////////////////////////////////////////////
// 1. Interpolate from  blockgrid to rectangular blockgrid
/////////////////////////////////////////////////////////////


bool contained_in_tri(D2vector lam) {
  if(lam.x < -0.1)                 return false;
  if(lam.y < -0.1)                 return false;
  if(lam.x + lam.y > 1.1)  return false;
  return true;
}

bool new_lam_worse(D2vector lam_old, D2vector lam_new) {
  if(MIN(lam_new) < MIN(lam_old) &&  MIN(lam_old) < -0.1) return true;
  if(MAX(lam_new) > MAX(lam_old) &&  MAX(lam_old) >  1.1) return true;
  return false;
}

Interpolate_on_structured_2Dgrid::Interpolate_on_structured_2Dgrid(int nx_, int ny_,
							           D2vector pWS_, D2vector pEN_,
							           Blockgrid2D& blockgrid_) {
  assert(nx_ > 1);
  assert(ny_ > 1);
  
  nx = nx_;
  ny = ny_;  
  blockgrid = &blockgrid_;  
    
  ug = blockgrid->Give_unstructured_grid();
    
  pWS = pWS_;
  pEN = pEN_;
  
  initialize();
}


Interpolate_on_structured_2Dgrid::Interpolate_on_structured_2Dgrid(Blockgrid2D* blockgrid_, double factorIncrease) {
  blockgrid = blockgrid_;
  ug = blockgrid->Give_unstructured_grid();  
    
  D2vector boxWS;
  D2vector boxEN;
  bool set = false;
    
  int Nmax =0;
  
  for(int id_quad=0;id_quad<ug->Give_number_rectangle();++id_quad) {
      int Nx = blockgrid->Give_Nx_rectangle(id_quad);
      int Ny = blockgrid->Give_Ny_rectangle(id_quad);

      for(int j=0;j<Ny;++j)
	  for(int i=0;i<Nx;++i) {
            // corner points of general hex-cell
	    D2vector cWS = blockgrid->Give_coord_rectangle(id_quad,i,  j  );
	    D2vector cES = blockgrid->Give_coord_rectangle(id_quad,i+1,j  );
	    D2vector cWN = blockgrid->Give_coord_rectangle(id_quad,i,  j+1);
	    D2vector cEN = blockgrid->Give_coord_rectangle(id_quad,i+1,j+1);


            // bounding box calculation 
	    if(set) {
	       boxWS.x = MIN(MIN(MIN(cWS.x,cES.x),MIN(cWN.x,cEN.x)),boxWS.x);
	       boxWS.y = MIN(MIN(MIN(cWS.y,cES.y),MIN(cWN.y,cEN.y)),boxWS.y);
	    
	       boxEN.x = MAX(MAX(MAX(cWS.x,cES.x),MAX(cWN.x,cEN.x)),boxEN.x);
	       boxEN.y = MAX(MAX(MAX(cWS.y,cES.y),MAX(cWN.y,cEN.y)),boxEN.y);
	    }
	    else {
	       boxWS.x = MIN(MIN(cWS.x,cES.x),MIN(cWN.x,cEN.x));
	       boxWS.y = MIN(MIN(cWS.y,cES.y),MIN(cWN.y,cEN.y));
	    
	       boxEN.x = MAX(MAX(cWS.x,cES.x),MAX(cWN.x,cEN.x));
	       boxEN.y = MAX(MAX(cWS.y,cES.y),MAX(cWN.y,cEN.y));
	       set = true;
	    }
	    
	    int N = MAX(blockgrid->Give_Nx_rectangle(id_quad),blockgrid->Give_Ny_rectangle(id_quad));
	    if(Nmax < N) Nmax = N;
	  }
  }
  
  nx = Nmax*factorIncrease;
  ny = Nmax*factorIncrease;  
  
  
  pWS = boxWS;
  pEN = boxEN;
  
  initialize();
}

Interpolate_on_structured_2Dgrid::~Interpolate_on_structured_2Dgrid() {
  delete[] ids_quad;
  delete[] ids_i;
  delete[] ids_j;

  delete[] typ_tri;

  delete[] lambda;
}

void Interpolate_on_structured_2Dgrid::initialize() {
  int typ;
    
  int ilmin, jlmin;
  int ilmax, jlmax;

  double factor = 0.1;
  //  double factor = 0.00001;

  D2vector lam;
  


  if(nx>1)
    hx = (pEN.x - pWS.x) / (nx-1);
  else  
    hx = 1.0;
  if(ny>1)
    hy = (pEN.y - pWS.y) / (ny-1);
  else  
    hy = 1.0;

  int num_total = nx * ny;

  D2vector cWS, cES;
  D2vector cWN, cEN;


  D2vector boxWS, boxEN;

  D2vector ploc;

  ids_quad = new int[num_total];

  ids_i = new int[num_total];
  ids_j = new int[num_total];

  typ_tri = new int[num_total];

  lambda = new D2vector[num_total];

  for(int i=0;i<num_total;++i) ids_quad[i] = -1;

  for(int id_quad=0;id_quad<ug->Give_number_rectangle();++id_quad) {
      int Nx = blockgrid->Give_Nx_rectangle(id_quad);
      int Ny = blockgrid->Give_Ny_rectangle(id_quad);

      for(int j=0;j<Ny;++j)
	  for(int i=0;i<Nx;++i) {
            // corner points of general hex-cell
	    cWS = blockgrid->Give_coord_rectangle(id_quad,i,  j  );
	    cES = blockgrid->Give_coord_rectangle(id_quad,i+1,j  );
	    cWN = blockgrid->Give_coord_rectangle(id_quad,i,  j+1);
	    cEN = blockgrid->Give_coord_rectangle(id_quad,i+1,j+1);


            // bounding box calculation 
	    boxWS.x = MIN(MIN(cWS.x,cES.x),MIN(cWN.x,cEN.x)) - factor *hx;
	    boxWS.y = MIN(MIN(cWS.y,cES.y),MIN(cWN.y,cEN.y)) - factor *hy;
	    
	    boxEN.x = MAX(MAX(cWS.x,cES.x),MAX(cWN.x,cEN.x)) + factor *hx;
	    boxEN.y = MAX(MAX(cWS.y,cES.y),MAX(cWN.y,cEN.y)) + factor *hy;
	    
	    // calculation of indices of a collection of cells of structured grid which contains bounding box
	    ilmin = Ganzzahliger_Anteil((boxWS.x - pWS.x) / hx);
	    jlmin = Ganzzahliger_Anteil((boxWS.y - pWS.y) / hy);
	    

	    ilmax = Ganzzahliger_Anteil((boxEN.x - pWS.x) / hx);
	    jlmax = Ganzzahliger_Anteil((boxEN.y - pWS.y) / hy);

	    /*
	    cout << " indices: "
		 << " ilmin: " << ilmin 
		 << " jlmin: " << jlmin 
		 << " klmin: " << klmin 
		 << " ilmax: " << ilmax 
		 << " jlmax: " << jlmax 
		 << " klmax: " << klmax 
		 << " boxWSD.x: " << boxWSD.x
		 << " cWSD.x: " << cWSD.x
		 << " Nx: " <<  Nx
		 << endl;
	    */
	    /*	    
bool now;
if(boxWSD.z < 0 && boxENT.z > 0.0 && boxWSD.y < 0.5 && boxENT.y > 0.5 && boxWSD.x < 1.0 && boxENT.x > 1.0 ) {
  cout << "\n \n WSD   : ";  boxWSD.Print();
  cout << "\n ENT  : ";  boxENT.Print();
  now = true;
 }
 else now = false;

 
 cout << " tt: " << boxWSD.x << " " << pWSD.x << " " << hx << endl;
 cout << (boxWSD.x - pWSD.x) << endl;

 cout << " z: " << 0.1 << " g: " << Ganzzahliger_Anteil(0.1) << endl;
 cout << " z: " << -0.1 << " g: " << Ganzzahliger_Anteil(-0.1) << endl;

 cout << " z: " << 5.1 << " g: " << Ganzzahliger_Anteil(5.1) << endl;
 cout << " z: " << -5.1 << " g: " << Ganzzahliger_Anteil(-5.1) << endl;
	    */

	    if(ilmin<0) ilmin=0;
	    if(jlmin<0) jlmin=0;

	    for(int il = ilmin; (il <= ilmax) && (il < nx);++il)
	      for(int jl = jlmin; (jl <= jlmax) && (jl < ny);++jl) {

		  ploc = D2vector(il * hx, jl * hy) + pWS;

		  //		  cout << "HI" << endl;

		  typ = -1;

		  lam = lambda_of_p_in_tri(ploc,cWN,cWS,cES);
		  if(contained_in_tri(lam)) typ=0;
		  else {
		    lam = lambda_of_p_in_tri(ploc,cEN,cWN,cES);
		    if(contained_in_tri(lam)) typ=1;
		  }

		  /*		  
		  cout << " typ " << typ << id_hex 
		       << " il: " << il 
		       << " jl: " << jl 
		       << " kl: " << kl 
		       << endl; 
		  */

		  if(typ!=-1) {
		    int ind_global;
		    ind_global = il+nx*(jl);      
		    bool stop;
		    stop=false;

		    if(ids_quad[ind_global]!=-1) {
		      stop=new_lam_worse(lambda[ind_global],lam);
		    }

		    if(stop==false) {
		      ids_quad[ind_global] = id_quad;
		      ids_i[ind_global] = i;
		      ids_j[ind_global] = j;
		      
		      typ_tri[ind_global] = typ;
		      
		      lambda[ind_global] = lam;
		    }
		    //go_on = false;
		  }

		  /*
		  cout << " out "
		       << " ilmin: " << ilmin
		       << " ilmax: " << ilmax
		       << " jlmin: " << jlmin
		       << " jlmax: " << jlmax
		       << " klmin: " << klmin
		       << " klmax: " << klmax;
		  cout << "\n   "; cWSD.Print();
		  cout << "\n   "; cESD.Print();
		  cout << "\n   "; cWND.Print();
		  cout << "\n   "; cEND.Print();
		  cout << "\n   "; cWST.Print();
		  cout << "\n   "; cEST.Print();
		  cout << "\n   "; cWNT.Print();
		  cout << "\n   "; cENT.Print();
		  cout << "\n   p: "; ploc.Print();

		  cout << "\n   : ";  boxWSD.Print();
		  cout << "\n   : ";  boxENT.Print();

		  cout << endl;
		  */
		}
	  }
  }


  for(int i=0;i<num_total;++i) {
    if(ids_quad[i]==-1) {
      // wir nehmen default value!!
      /*
      cout << i 
	   << " Error: Interpolate_on_structured_2Dgrid: I cannot interpolate all data!"
	   << endl;
      ids_hex[i] = 0;
      */
    }
    else {
      //cout << i << " Interpolate_on_structured_2Dgrid: o.k.!" << endl;
    }
  }
}



void Interpolate_on_structured_2Dgrid::interpolateAbs(Variable2D<std::complex<double> >& u, 
                                                      double* data, double defaultInterpolation) {
    int i,j, id_quad, typ;
    int ind_global;
    
    
    for(int id=0;id<ug->Give_number_rectangle();++id)
        u.template Update<rectangleEl>(id);
    
    for(int js = 0; js < ny;++js)
        for(int is = 0; is < nx;++is) {
            ind_global = is+nx*js;
            i = ids_i[ind_global];
            j = ids_j[ind_global];
            id_quad = ids_quad[ind_global];
            
            if(id_quad < 0) data[ind_global] = defaultInterpolation;
            else {
                Give_corner_data_of_quad<std::complex<double> > du(u, id_quad, i,j);
                
                typ = typ_tri[ind_global];
                
                if(typ==0) data[ind_global] 
                      = absExpression(interpolate_in_tri(lambda[ind_global],du.WN(),du.WS(),du.ES()));
                if(typ==1) data[ind_global] 
                      = absExpression(interpolate_in_tri(lambda[ind_global],du.EN(),du.WN(),du.ES()));
            }
        }
}



///////////////////////////////////////////////////////////////////////////////////
//  ColorMap
///////////////////////////////////////////////////////////////////////////////////

      
uint ColorMap::colormap[ColormapSize];  
          
void ColorMap::makeColorMap() {
   for(int i = 0; i < ColormapSize; ++i) {
       colormap[i] = rgbFromWaveLength(380.0 + (i * 400.0 / ColormapSize));
   }  
}


QImage ColorMap::getImageColorMap(QSize resultSize) {
        makeColorMap();
        int widthI = resultSize.width();
        int heightI = resultSize.height();
	
        QImage image(resultSize, QImage::Format_RGB32);
  
        for(int y = 0; y < heightI; ++y) {
            uint *scanLine = reinterpret_cast<uint *>(image.scanLine(y));

            for(int x = 0; x < widthI; ++x) {
	        double coord_y = ((double)(heightI-1-y))/(heightI-1);
       		int valueIntU =	coord_y * ColormapSize;
		*scanLine++ = colormap[valueIntU % ColormapSize];
	    }
	}
        return image;
}

uint ColorMap::rgbFromWaveLength(double wave) {
    double r = 0.0;
    double g = 0.0;
    double b = 0.0;qRgb(int(r * 255), int(g * 255), int(b * 255));

    if (wave >= 380.0 && wave <= 440.0) {
        r = -1.0 * (wave - 440.0) / (440.0 - 380.0);
        b = 1.0;
    } else if (wave >= 440.0 && wave <= 490.0) {
        g = (wave - 440.0) / (490.0 - 440.0);
        b = 1.0;
    } else if (wave >= 490.0 && wave <= 510.0) {
        g = 1.0;
        b = -1.0 * (wave - 510.0) / (510.0 - 490.0);
    } else if (wave >= 510.0 && wave <= 580.0) {
        r = (wave - 510.0) / (580.0 - 510.0);
        g = 1.0;
    } else if (wave >= 580.0 && wave <= 645.0) {
        r = 1.0;
        g = -1.0 * (wave - 645.0) / (645.0 - 580.0);
    } else if (wave >= 645.0 && wave <= 780.0) {
        r = 1.0;
    }


    double s = 1.0;
    if (wave > 700.0)
        s = 0.3 + 0.7 * (780.0 - wave) / (780.0 - 700.0);
    else if (wave <  420.0)
        s = 0.3 + 0.7 * (wave - 380.0) / (420.0 - 380.0);

    r = std::pow(r * s, 0.8);
    g = std::pow(g * s, 0.8);
    b = std::pow(b * s, 0.8);
    return qRgb(int(r * 255), int(g * 255), int(b * 255));
}

   


///////////////////////////////////////////////////////////////////////////////////
//  QImageMaker
///////////////////////////////////////////////////////////////////////////////////
	  
	  /*
QImageMaker::QImageMaker(Variable2D<double>& u, int thicknessContourLines_) :
QImageMakerAbstract(u.Give_blockgrid(),thicknessContourLines) {
   maxU = Maximum(u);
   minU = Minimum(u);    
      
   makeColorMap();
   
   interpolate(u,data,0.0);
}	  
	*/
      
	/*
QImageMakerAbs::QImageMakerAbs(Variable2D< std::complex<double> >& u, int thicknessContourLines_) :
QImageMakerAbstract(u.Give_blockgrid(),thicknessContourLines) {
   maxU = L_infty(u);
   minU = 0.0;   
    
//   Function2d1<double,std::complex<double> > absFC(absComplexHere);    
//   maxU = Maximum(absFC(u)); //L_infty(u);
//   minU = Minimum(absFC(u));
//   cout << " Maximum: " << maxU << " Minimum: " << minU << endl;   
//   Variable2D<double> u_test(*u.Give_blockgrid());
//   u_test = absFC(u);

   makeColorMap();
   
   interpolateAbs(u,data,0.0);
}		
	*/
    
    /*
QImageMakerAbstract::QImageMakerAbstract(Blockgrid2D* bg_, int thicknessContourLines_) : Interpolate_on_structured_2Dgrid(bg_) {
   thicknessContourLines = thicknessContourLines_;
   //Variable2D<double>* u = &u_; 
   //ug = u->Give_Ug();
   bg = bg_;
    
   valueBackground = qRgb(255, 255, 255);
   contourColor    = qRgb(0,   0,   0  );
   valueEmpty      = qRgb(248, 248, 255);

   data = new double[nx*ny];
}
*/
    
QImageMaker::QImageMaker(Variable2D<double>& u, int thicknessContourLines_) :
             Interpolate_on_structured_2Dgrid(u.Give_blockgrid()) {
   bg = u.Give_blockgrid();
   thicknessContourLines = thicknessContourLines_;   

   maxU = Maximum(u);
   minU = Minimum(u);    
      
   makeColorMap();
   
   valueBackground = qRgb(255, 255, 255);
   contourColor    = qRgb(0,   0,   0  );
   valueEmpty      = qRgb(248, 248, 255);

   data = new double[nx*ny];   
   
   interpolate(u,data,0.0);
}  

QImageMaker::QImageMaker(Variable2D<std::complex<double> >& u, int thicknessContourLines_) :
            Interpolate_on_structured_2Dgrid(u.Give_blockgrid()) {
   bg = u.Give_blockgrid();
   thicknessContourLines = thicknessContourLines_;   
   
   maxU = L_infty(u);
   minU = 0.0;         

   makeColorMap();
   
   valueBackground = qRgb(255, 255, 255);
   contourColor    = qRgb(0,   0,   0  );
   valueEmpty      = qRgb(248, 248, 255);

   data = new double[nx*ny];   
   
   interpolateAbs(u,data,0.0);
}  

//QImageMakerAbstract::~QImageMakerAbstract() {
QImageMaker::~QImageMaker() {
  delete[] data; 
}

QImage QImageMaker::makeImage(QSize resultSize, unsigned int numberContours, 
//QImage QImageMakerAbstract::makeImage(QSize resultSize, unsigned int numberContours, 
                                      bool showData, bool origionalSize) {
       //double scaleFactor = this->scaleFactor;
       //double centerX = this->centerX;
       //double centerY = this->centerY;
	
 
	std::vector<int> contourVals;
	contourVals.resize(numberContours);
	
	for(int i=0;i<numberContours;++i) {
	    contourVals.at(i) = ColormapSize * (i+1)/(numberContours+2) ;
	}
	
	
	/*
	resultSize = QSize(400,550);
        centerX = -0.637011;
        centerY = -0.0395159;
        scaleFactor = 0.00403897;
*/
	
        int widthI = resultSize.width();
        int heightI = resultSize.height();

	if(origionalSize) {
	   double ratioGiven = (double)widthI / (double)heightI;
           double ratioHere  = (pEN.x - pWS.x) / (pEN.y - pWS.y);
        
	   if(ratioHere >= ratioGiven) { //Das Bild ist zu breit
	      resultSize.setHeight(widthI / ratioHere);
	      heightI = resultSize.height();
  	   }
	   else {
	      resultSize.setWidth(heightI * ratioHere);
	      widthI = resultSize.width();	  
	   }
	}
	
	double Ihx = (pEN.x - pWS.x) / (widthI-1);
	double Ihy = (pEN.y - pWS.y) / (heightI-1);
	
	
        QImage image(resultSize, QImage::Format_RGB32);
  
        for(int y = 0; y < heightI; ++y) {
            uint *scanLine = reinterpret_cast<uint *>(image.scanLine(y));

            for(int x = 0; x < widthI; ++x) {
	      
	        double coord_x = x*Ihx;
	        double coord_y = y*Ihy;
		
                if(bg->inStrechedDomain(pWS + D2vector(coord_x,coord_y))==false) {  
//                   *scanLine++ = colormap[valueBackground % ColormapSize];
                   *scanLine++ = valueBackground;		  
		}
		else {
                   int i = coord_x / hx;   
                   int j = coord_y / hy;   
    
                   if(i>=nx-1) i=nx-2;  
		   if(j>=ny-1) j=ny-2;
  
                //cout << "i: " << i << " j: " << j << endl;
    
    
                   double uWS = data[i    +nx*(j)];
                   double uES = data[(i+1)+nx*(j)];
                   double uWN = data[i    +nx*(j+1)];
                   double uEN = data[(i+1)+nx*(j+1)];
  
		   bool WSin = bg->inStrechedDomain(pWS + D2vector(i*hx,j*hy));
		   bool ESin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,j*hy));
		   bool WNin = bg->inStrechedDomain(pWS + D2vector(i*hx,(j+1)*hy));
		   bool ENin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,(j+1)*hy));
		   if(WSin==false || ESin==false || WNin==false || ENin==false) {
		      double valueOk;
		      if(WSin) valueOk=uWS;
		      if(WNin) valueOk=uWN;
		      if(ESin) valueOk=uES;
		      if(ENin) valueOk=uEN;
                      if(WSin==false) uWS = valueOk;
                      if(WNin==false) uWN = valueOk;
                      if(ESin==false) uES = valueOk;
                      if(ENin==false) uEN = valueOk;
		   }
		   
                // assert( (i+1)+nx*(j+1) < nx*ny);
  
                   double locX = coord_x / hx - i;
                   double locY = coord_y / hy - j;
  
                   double valueU = uWS * (1.0 - locX) * (1.0 - locY) +
                                   uES *        locX  * (1.0 - locY) +
                                   uWN * (1.0 - locX) *        locY  +
                                   uEN *        locX  *        locY;
       		   int valueIntU =	(valueU - minU) / (maxU - minU) * ColormapSize;
		
		   bool notFoundContour = true;
		   if(numberContours>0) {
		      for(int i=0;i<numberContours && notFoundContour;++i) {
//		          if(valueIntU == contourVals.at(i)) {
		          if(ABS(valueIntU-contourVals.at(i)) <= thicknessContourLines) {			
//			     *scanLine++ = colormap[contourColor % ColormapSize];
			     *scanLine++ = contourColor;
			     notFoundContour = false;
			  }
		      }
		   }
		   if(notFoundContour) {
		      if(showData) *scanLine++ = colormap[valueIntU % ColormapSize];
		      else         *scanLine++ = valueEmpty;
		   }
		}
				
  	        //*scanLine++ = colormap[rgbFromU(valueU) % ColormapSize];
	        /*
                if(x > 0 && y > 0)
                   *scanLine++ = colormap[6 % ColormapSize];
                else {
                   *scanLine++ = qRgb(0, 0, 0);
		}
		*/
	    }
	}
        return image;
}


///////////////////////////////////////////////////////////////////////////////////
//  VectorOfArrowsCreator
///////////////////////////////////////////////////////////////////////////////////

VectorOfArrowsCreator::VectorOfArrowsCreator(Variable2D<double>& ux, Variable2D<double>& uy) 
                : Interpolate_on_structured_2Dgrid(ux.Give_blockgrid()) {
  //ug = u->Give_Ug();
  bg = ux.Give_blockgrid();
  
  valueBackground = qRgb(255, 255, 255);
  contourColor    = qRgb(0,   0,   0  );
  valueEmpty      = qRgb(248, 248, 255);
  
  makeColorMap();
  
  dataX = new double[nx*ny];
  dataY = new double[nx*ny];
   
  interpolate(ux,dataX,0.0);
  interpolate(uy,dataY,0.0);   
}

VectorOfArrowsCreator::~VectorOfArrowsCreator() {
  delete[] dataX; 
  delete[] dataY;   
}
 
void VectorOfArrowsCreator::createVectorList(QVector<QPointF>& startP, QVector<QPointF>& endP,
                   QSize resultSize, unsigned int numberArrowX, unsigned int numberArrowY, 
		   bool origionalSize) {
      
   startP.resize(0);
   endP.resize(0);   
      
   int widthI = resultSize.width();
   int heightI = resultSize.height();
   
   if(origionalSize) {
	   double ratioGiven = (double)widthI / (double)heightI;
           double ratioHere  = (pEN.x - pWS.x) / (pEN.y - pWS.y);
        
	   if(ratioHere >= ratioGiven) { //Das Bild ist zu breit
	      resultSize.setHeight(widthI / ratioHere);
	      heightI = resultSize.height();
  	   }
	   else {
	      resultSize.setWidth(heightI * ratioHere);
	      widthI = resultSize.width();	  
	   }
    }
   
    double Phx = widthI  / numberArrowX;
    double Phy = heightI / numberArrowY;
   
    double Ph = MIN(Phx,Phy) * 0.5;
    
    double Ihx = (pEN.x - pWS.x) / numberArrowX;
    double Ihy = (pEN.y - pWS.y) / numberArrowY;   
   
    for(int y = 0; y < numberArrowX; ++y) {
        for(int x = 0; x < numberArrowY; ++x) {  
	    double coord_x = (x+0.5)*Ihx;
	    double coord_y = (y+0.5)*Ihy;
	    
            if(bg->inStrechedDomain(pWS + D2vector(coord_x,coord_y))) {
                   int i = coord_x / hx;   
                   int j = coord_y / hy;   
    
                   if(i>=nx-1) i=nx-2;  
		   if(j>=ny-1) j=ny-2;
  
                   double uWSx = dataX[i    +nx*(j)];
                   double uESx = dataX[(i+1)+nx*(j)];
                   double uWNx = dataX[i    +nx*(j+1)];
                   double uENx = dataX[(i+1)+nx*(j+1)];

                   double uWSy = dataY[i    +nx*(j)];
                   double uESy = dataY[(i+1)+nx*(j)];
                   double uWNy = dataY[i    +nx*(j+1)];
                   double uENy = dataY[(i+1)+nx*(j+1)];
		   
		   bool WSin = bg->inStrechedDomain(pWS + D2vector(i*hx,j*hy));
		   bool ESin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,j*hy));
		   bool WNin = bg->inStrechedDomain(pWS + D2vector(i*hx,(j+1)*hy));
		   bool ENin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,(j+1)*hy));
		   if(WSin==false || ESin==false || WNin==false || ENin==false) {
		      double valueOkX, valueOkY;
		      if(WSin) { valueOkX=uWSx;  valueOkY=uWSy;  }
		      if(WNin) { valueOkX=uWNx;  valueOkY=uWNy;  }
		      if(ESin) { valueOkX=uESx;  valueOkY=uESy;  }
		      if(ENin) { valueOkX=uENx;  valueOkY=uENy;  }
                      if(WSin==false) { uWSx = valueOkX;   uWSy = valueOkY; }
                      if(WNin==false) { uWNx = valueOkX;   uWNy = valueOkY; } 
                      if(ESin==false) { uESx = valueOkX;   uESy = valueOkY; }
                      if(ENin==false) { uENx = valueOkX;   uENy = valueOkY; }
		   }
		   	        
                   double locX = coord_x / hx - i;
                   double locY = coord_y / hy - j;
  
                   double valueUx = uWSx * (1.0 - locX) * (1.0 - locY) +
                                    uESx *        locX  * (1.0 - locY) +
                                    uWNx * (1.0 - locX) *        locY  +
                                    uENx *        locX  *        locY;
	      
                   double valueUy = uWSy * (1.0 - locX) * (1.0 - locY) +
                                    uESy *        locX  * (1.0 - locY) +
                                    uWNy * (1.0 - locX) *        locY  +
                                    uENy *        locX  *        locY;
                   double norm = sqrt(valueUx*valueUx + valueUy*valueUy);	   
		   
		   if(norm>1.0e-10) {
		      valueUx = valueUx / norm * Ph;
		      valueUy = valueUy / norm * Ph;
		      
		      QPointF A((x+0.5)*Phx, (y+0.5)*Phy);
		      QPointF B((x+0.5)*Phx + valueUx, (y+0.5)*Phy + valueUy);
		      
		      startP.append(A);
		      endP.append(B);		      
		   }
	    }
	}
    }
}
 

VectorOfArrowsCreatorComplex::VectorOfArrowsCreatorComplex(Variable2D<std::complex<double> >& ux, 
                                                           Variable2D<std::complex<double> >& uy) 
                : Interpolate_on_structured_2Dgrid(ux.Give_blockgrid()) {
  //ug = u->Give_Ug();
  bg = ux.Give_blockgrid();
  
  valueBackground = qRgb(255, 255, 255);
  contourColor    = qRgb(0,   0,   0  );
  valueEmpty      = qRgb(248, 248, 255);
  
  makeColorMap();
  
  dataX = new std::complex<double>[nx*ny];
  dataY = new std::complex<double>[nx*ny];
   
  interpolate(ux,dataX,std::complex<double>(0.0,0.0));
  interpolate(uy,dataY,std::complex<double>(0.0,0.0));   
}

VectorOfArrowsCreatorComplex::~VectorOfArrowsCreatorComplex() {
  delete[] dataX; 
  delete[] dataY;   
}
 
void VectorOfArrowsCreatorComplex::createVectorList(QVector<QPointF>& startP, QVector<QPointF>& endP,
                   QSize resultSize, unsigned int numberArrowX, unsigned int numberArrowY, 
		   bool origionalSize) {
      
   startP.resize(0);
   endP.resize(0);   
      
   int widthI = resultSize.width();
   int heightI = resultSize.height();
   
   if(origionalSize) {
	   double ratioGiven = (double)widthI / (double)heightI;
           double ratioHere  = (pEN.x - pWS.x) / (pEN.y - pWS.y);
        
	   if(ratioHere >= ratioGiven) { //Das Bild ist zu breit
	      resultSize.setHeight(widthI / ratioHere);
	      heightI = resultSize.height();
  	   }
	   else {
	      resultSize.setWidth(heightI * ratioHere);
	      widthI = resultSize.width();	  
	   }
    }
   
    double Phx = widthI  / numberArrowX;
    double Phy = heightI / numberArrowY;
   
    double Ph = MIN(Phx,Phy) * 0.45;
    
    double Ihx = (pEN.x - pWS.x) / numberArrowX;
    double Ihy = (pEN.y - pWS.y) / numberArrowY;   
   
    for(int y = 0; y < numberArrowX; ++y) {
        for(int x = 0; x < numberArrowY; ++x) {  
	    double coord_x = (x+0.5)*Ihx;
	    double coord_y = (y+0.5)*Ihy;
	    
            if(bg->inStrechedDomain(pWS + D2vector(coord_x,coord_y))) {
                   int i = coord_x / hx;   
                   int j = coord_y / hy;   
    
                   if(i>=nx-1) i=nx-2;  
		   if(j>=ny-1) j=ny-2;
  
                   std::complex<double> uWSx = dataX[i    +nx*(j)];
                   std::complex<double> uESx = dataX[(i+1)+nx*(j)];
                   std::complex<double> uWNx = dataX[i    +nx*(j+1)];
                   std::complex<double> uENx = dataX[(i+1)+nx*(j+1)];

                   std::complex<double> uWSy = dataY[i    +nx*(j)];
                   std::complex<double> uESy = dataY[(i+1)+nx*(j)];
                   std::complex<double> uWNy = dataY[i    +nx*(j+1)];
                   std::complex<double> uENy = dataY[(i+1)+nx*(j+1)];
		   
		   bool WSin = bg->inStrechedDomain(pWS + D2vector(i*hx,j*hy));
		   bool ESin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,j*hy));
		   bool WNin = bg->inStrechedDomain(pWS + D2vector(i*hx,(j+1)*hy));
		   bool ENin = bg->inStrechedDomain(pWS + D2vector((i+1)*hx,(j+1)*hy));
		   if(WSin==false || ESin==false || WNin==false || ENin==false) {
		      std::complex<double> valueOkX, valueOkY;
		      if(WSin) { valueOkX=uWSx;  valueOkY=uWSy;  }
		      if(WNin) { valueOkX=uWNx;  valueOkY=uWNy;  }
		      if(ESin) { valueOkX=uESx;  valueOkY=uESy;  }
		      if(ENin) { valueOkX=uENx;  valueOkY=uENy;  }
                      if(WSin==false) { uWSx = valueOkX;   uWSy = valueOkY; }
                      if(WNin==false) { uWNx = valueOkX;   uWNy = valueOkY; } 
                      if(ESin==false) { uESx = valueOkX;   uESy = valueOkY; }
                      if(ENin==false) { uENx = valueOkX;   uENy = valueOkY; }
		   }
		   	        
                   double locX = coord_x / hx - i;
                   double locY = coord_y / hy - j;
  
                   std::complex<double> valueUx = uWSx * (1.0 - locX) * (1.0 - locY) +
                                    uESx *        locX  * (1.0 - locY) +
                                    uWNx * (1.0 - locX) *        locY  +
                                    uENx *        locX  *        locY;
	      
                   std::complex<double> valueUy = uWSy * (1.0 - locX) * (1.0 - locY) +
                                    uESy *        locX  * (1.0 - locY) +
                                    uWNy * (1.0 - locX) *        locY  +
                                    uENy *        locX  *        locY;
                                    
                   double signX = 1.0;
                   double signY = 1.0;
//test NOW		   
                   if(valueUx.real() < 0.0) signX = -1.0;
                   if(valueUy.real() < 0.0) signY = -1.0;
                   
                   double absUx = absExpression(valueUx);
                   double absUy = absExpression(valueUy);
                   double norm  = sqrt(absUx*absUx + absUy*absUy);
                   double cosAngle = signX * absUx / norm;
                   double sinAngle = signY * absUy / norm;
                   
		   if(norm>1.0e-10) {
		      cosAngle = cosAngle * Ph;
		      sinAngle = sinAngle * Ph;
		      
		      QPointF A((x+0.5)*Phx, (y+0.5)*Phy);
		      QPointF B((x+0.5)*Phx + cosAngle, (y+0.5)*Phy + sinAngle);
		      
		      startP.append(A);
		      endP.append(B);		      
		   }
	    }
	}
    }
}