/****************************************************************************
 **
 ** Copyright (C) 2011 Christian B. Huebschle & George M. Sheldrick
 ** All rights reserved.
 ** Contact: chuebsch@moliso.de
 **
 ** This file is part of the ShelXle
 **
 ** This file may be used under the terms of the GNU Lesser
 ** General Public License version 2.1 as published by the Free Software
 ** Foundation and appearing in the file COPYING included in the
 ** packaging of this file.  Please review the following information to
 ** ensure the GNU Lesser General Public License version 2.1 requirements
 ** will be met: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html.
 **
 **
 ****************************************************************************/
#include "fourxle.h"
#include <QtCore>
#include <QtGui>
#include "deprecation.h"
#include "itsme.h"
//üäöß 
FourXle::FourXle(Molecule *mole_, ChGL *chgl_,QToolBar *toolView, double resol, double wght){
  //  printf("FourXle\n");
  ownfcf6_=NULL;
  mole =mole_;
  f000=0.0;
  overridef000=-1.0;
  n1=n2=n3=n4=n5=0;
  voxelstr.clear();
  chgl=chgl_;
  MYFCF=false;
  oldatomsize=-1;
  chgl->habzutun=false;
  map_radius=5.0;
  maptrunc = 2;
  mapcontrol=NULL;
  resDir=QDir::current();
  // for (int i=0; i<37; i++) chgl->foubas[i]=0;
  /*
     QVector< QVector< float > > chgl->fchgl->fVertexes;
     QVector< QVector< float > > chgl->fNormals;
     */
  chgl->fVertexes.resize(28);
  chgl->fNormals.resize(28);
  HKLMX=200;
  datfo_fc=datfo=NULL;
  nodex=nodey=nodez=NULL;noGoMap=NULL;
  nodeX=nodeY=nodeZ=NULL;
  urs=V3(0,0,0);
  nr=0;
  nc=0;
  sigma[0]=sigma[1]=iso[0]=iso[1]=0;
  chgl->lintrans=0.8;
  chgl->linwidth=0.5;
  rr=resol;
  rw=wght; 
  chgl->foact=toolView->addAction(QIcon(":/xle-icons/fomaps.svg"),"toggle Fo map",chgl,SLOT(updateGL()));
  chgl->foact->setWhatsThis("<h1>Toggle Fo map <img src=\":/xle-icons/fomaps.svg\"></h1> Toggle F-observed map on/off.");
  connect(chgl->foact,SIGNAL(destroyed(QObject*)),this,SLOT(foactDestroyed()));
  chgl->foact->setCheckable(true);
  chgl->fofcact=toolView->addAction(QIcon(":/xle-icons/diffmap.svg"),"toggle Fo-Fc map",chgl,SLOT(updateGL()));
  chgl->fofcact->setWhatsThis("<h1>Toggle Fo-Fc map<img src=\":/xle-icons/diffmap.svg\"></h1> Toggle (F-observed)-(F-calculated) map on/off. Fo-Fc map is also known as difference Fourier map.");
  connect(chgl->fofcact,SIGNAL(destroyed(QObject*)),this,SLOT(fofcactDestroyed()));
  chgl->fofcact->setCheckable(true);
  chgl->fofcact->setVisible(false);
  chgl->foact->setVisible(false);
  doMaps = new QCheckBox("Calculate Maps");
  doMaps->setChecked(true);
  doMaps->setWhatsThis("<h1>Calculate maps</h1> Enables Fo and Fo-Fc map calculation. This can only be performed if LIST is set to 6.");
  keepIso = new QCheckBox("Keep map iso values");
  keepIso->setChecked(false);
  keepIso->setVisible(false);
  connect(keepIso,SIGNAL(toggled(bool)), this,SLOT(keepIsoValues(bool )));
  connect(chgl,SIGNAL(destroyed(QObject*)),this,SLOT(chglDestroyed()));
  char *ENUS=setlocale(LC_ALL,"C");
  if (!ENUS)ENUS=setlocale(LC_ALL,"en_US.UTF-8"); 
  {
    buttonBoxMC = new QDialogButtonBox( QDialogButtonBox::Apply | QDialogButtonBox::Cancel);
    applyMC          = buttonBoxMC->button(QDialogButtonBox::Apply);
    applyDF          = buttonBoxMC->addButton("Apply defaults",QDialogButtonBox::ActionRole);
    wl= new QLabel("Factor to downweight weak data");
    pl= new QLabel("Map precision: ");
    dl= new QLabel("Fo-Fc map: ");
    ol= new QLabel("F-observed map: ");
    sl= new QLabel("Map truncation type: ");
    lw= new QLabel("Line width: ");
    ltz= new QLabel("Line transparency: ");
    md = new QDialog();
    md->setWindowTitle("Map Control");
    connect(applyMC , SIGNAL(clicked()), this , SLOT(controlMap()));
    connect(applyMC , SIGNAL(clicked()), this , SLOT(openMapControl()));
    connect(applyDF , SIGNAL(clicked()), this , SLOT(mapDefault()));
    connect(buttonBoxMC, SIGNAL(rejected()), md, SLOT(hide()));

    mapprec = new QDoubleSpinBox(md);
    mapprec->setMinimum(1.1);
    mapprec->setMaximum(8.0);
    mapprec->setValue(rr);
    mapprec->setSingleStep(0.1);

    lineTrans = new QDoubleSpinBox(md);
    lineTrans->setMinimum(0.1);
    lineTrans->setMaximum(1.0);
    lineTrans->setValue(chgl->lintrans);
    lineTrans->setSingleStep(0.1);

    lineWidth = new QDoubleSpinBox(md);
    lineWidth->setMinimum(0.15);
    lineWidth->setMaximum(4.0);
    lineWidth->setValue(chgl->linwidth);
    lineWidth->setSingleStep(0.1);

    weak = new QDoubleSpinBox(md);
    weak->setMinimum(0.1);
    weak->setMaximum(8.0);
    weak->setValue(rw);
    weak->setSingleStep(0.1);

    difmaps = new QDoubleSpinBox(md);
    difmaps->setMinimum(0.025);
    difmaps->setMaximum(9999);
    difmaps->setValue(fabs(iso[1]));
    difmaps->setSingleStep(fabs(iso[1])/10.0); 
    difmaps->setSuffix("e/A^3");

    fomaps = new QDoubleSpinBox();
    fomaps->setMinimum(0.025);
    fomaps->setMaximum(9999);
    fomaps->setValue(iso[0]);
    fomaps->setSingleStep(iso[0]/10.0);
    fomaps->setSuffix("e/A^3");

    cutrange = new QDoubleSpinBox(md);
    cutrange->setMinimum(0.5);
    cutrange->setMaximum(3.0);
    cutrange->setValue(1.41);
    cutrange->setSingleStep(0.1);
    cutrange->setSuffix("A");
    cutrange->setPrefix("cut range ");

    maprad = new QDoubleSpinBox();
    maprad->setValue(map_radius);
    maprad->setMinimum(0);
    maprad->setMaximum(25.0);
    maprad->setSuffix("A");
    maprad->setPrefix("Map radius: ");
    maprad->setSingleStep(0.3);
    mapSchnitt=  new QComboBox();
    mapSchnitt->addItem("one complete unit cell",0);
    mapSchnitt->addItem("sphere around rotation center",1);
    mapSchnitt->addItem("cut range around visible atoms or peaks",2);
    mapSchnitt->setCurrentIndex(maptrunc);
    connect(mapSchnitt,SIGNAL(currentIndexChanged(int)),this,SLOT(schnittart(int)));
    chgl->lighting = new QCheckBox("Lighting");
    chgl->lighting->setChecked(true);
    chgl->niceTrans= new QCheckBox("Undistorted transparency");
    chgl->niceTrans->setChecked(false);
    chgl->fillMap= new QCheckBox("Filled map surface");
    chgl->fillMap->setChecked(false);
    mdl = new QGridLayout();
    mdl->addWidget(mapprec,0,1);
    mdl->addWidget(weak,1,1);
    mdl->addWidget(difmaps,2,1);
    mdl->addWidget(fomaps,3,1);
    mdl->addWidget(maprad,6,2);
    mdl->addWidget(mapSchnitt,6,1);
    mdl->addWidget(cutrange,6,2);
    mdl->addWidget(lineTrans,7,1);
    mdl->addWidget(lineWidth,8,1);
    mdl->addWidget(chgl->lighting,9,1);
    mdl->addWidget(chgl->niceTrans,9,0);
    mdl->addWidget(chgl->fillMap,9,2);
    connect(lineTrans,SIGNAL(valueChanged(double)),this,SLOT(updateLines()));
    connect(lineWidth,SIGNAL(valueChanged(double)),this,SLOT(updateLines()));

    connect(chgl->fillMap,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));
    connect(chgl->niceTrans,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));
    connect(chgl->lighting,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));    
    mdl->addWidget(pl,0,0);
    mdl->addWidget(wl,1,0);
    mdl->addWidget(dl,2,0);
    mdl->addWidget(ol,3,0);
    //mdl->addWidget(rl,5,0);
    mdl->addWidget(sl,6,0);
    mdl->addWidget(ltz,7,0);
    mdl->addWidget(lw,8,0);

    QTabWidget *tabWidget = new QTabWidget(md);
    QWidget *gwmd,*admd;
    gwmd = new QWidget();

    gwmd->setLayout(mdl);
    admd = new QWidget();
    QVBoxLayout *mdmain= new QVBoxLayout;
    tabWidget->addTab(gwmd,"General");
    tabWidget->addTab(admd,"Advanced");

    mole->einstellung->beginGroup("Colors");
    chgl->fopc=mole->einstellung->value("FourFoPlusColor",   QColor("#0000FF")).value<QColor>();
    chgl->fomc=mole->einstellung->value("FourFoMinusColor",  QColor("#FFA600")).value<QColor>();
    chgl->dipc=mole->einstellung->value("FourDiffPlusColor", QColor("#00B200")).value<QColor>();
    chgl->dimc=mole->einstellung->value("FourDiffMinusColor",QColor("#CC0000")).value<QColor>();
    mole->einstellung->endGroup();

    foPlusButton =new QPushButton("Fo positive color");
    foPlusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
          "border: 0px"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg( chgl->fopc.name())
        .arg( chgl->fopc.darker(200).name())
        .arg( chgl->fopc.lighter(100).name())
        .arg("#ffffff"));
    connect(foPlusButton,SIGNAL(pressed()),this,SLOT(colorDLGFOP()));
    defColorButton =new QPushButton("Default color scheme");
    connect(defColorButton,SIGNAL(pressed()),this,SLOT(defaultColors()));
    foMinusButton =new QPushButton("Fo negative color");
    foMinusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
          "border: 0px"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg( chgl->fomc.name())
        .arg( chgl->fomc.darker(200).name())
        .arg( chgl->fomc.lighter(100).name())
        .arg("#ffffff"));
    connect(foMinusButton,SIGNAL(pressed()),this,SLOT(colorDLGFOM()));

    diffMinusButton =new QPushButton("Fo-Fc negative color");
    diffMinusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
          "border: 0px"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg( chgl->dimc.name())
        .arg( chgl->dimc.darker(200).name())
        .arg( chgl->dimc.lighter(100).name())
        .arg("#ffffff"));
    connect(diffMinusButton,SIGNAL(pressed()),this,SLOT(colorDLGDIM()));

    diffPlusButton =new QPushButton("Fo-Fc positive color");
    diffPlusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
          "border: 0px"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg( chgl->dipc.name())
        .arg( chgl->dipc.darker(200).name())
        .arg( chgl->dipc.lighter(100).name())
        .arg("#ffffff"));
    connect(diffPlusButton,SIGNAL(pressed()),this,SLOT(colorDLGDIP()));

    amdl = new QVBoxLayout();
    amdl->addWidget(foPlusButton);
    amdl->addWidget(foMinusButton);
    amdl->addWidget(diffPlusButton);
    amdl->addWidget(diffMinusButton);
    amdl->addWidget(defColorButton);
    jnkButton = new QPushButton("fractal dimension analysis");
    amdl->addWidget(jnkButton);
    admd->setLayout(amdl);
    mdmain->addWidget(tabWidget);
    mdmain->addWidget(buttonBoxMC);

    md->setLayout(mdmain);
    md->setModal (false);
    ownfcf6=new QCheckBox("fcf6 by hkl");
    ownfcf6->setChecked(MYFCF);
    ownfcf6_=new QCheckBox("Compute fcf6 file from structure factors and observed data from hkl file");
    ownfcf6_->setChecked(MYFCF);
    amdl->addWidget(ownfcf6_);
    connect(ownfcf6_,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
    connect(doMaps,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
    mole->einstellung->beginGroup("Window");//mole->einstellung->beginGroup("");
    if (mole->einstellung->contains("Fo-Map")) chgl->foact->setChecked(mole->einstellung->value("Fo-Map").toBool());
    if (mole->einstellung->contains("Fc-Fo-Map")) chgl->fofcact->setChecked(mole->einstellung->value("Fc-Fo-Map").toBool());
    if (mole->einstellung->contains("Map-truncation-type")) maptrunc=mole->einstellung->value("Map-truncation-type").toInt();
    if (mole->einstellung->contains("Map-Line-Transparency")) chgl->lintrans = mole->einstellung->value("Map-Line-Transparency").toDouble();
    if (mole->einstellung->contains("Map-Line-Width")) chgl->linwidth = mole->einstellung->value("Map-Line-Width").toDouble();
    if (mole->einstellung->contains("Map-Resolution")) rr = mole->einstellung->value("Map-Resolution").toDouble();
    if (mole->einstellung->contains("Map-Weight")) rw = mole->einstellung->value("Map-Weight").toDouble();
    if (mole->einstellung->contains("Map-Lighting")) chgl->lighting->setChecked(mole->einstellung->value("Map-Lighting").toBool());
    if (mole->einstellung->contains("Nice-Transparency")) chgl->niceTrans->setChecked(mole->einstellung->value("Nice-Transparency").toBool());
    if (mole->einstellung->contains("Solid-Map")) chgl->fillMap->setChecked(mole->einstellung->value("Solid-Map").toBool());
    if (mole->einstellung->contains("Calculate-Maps")) doMaps->setChecked(mole->einstellung->value("Calculate-Maps").toBool());
    if (mole->einstellung->contains("KeepMapIsoValues")){
      keepIso->setChecked(mole->einstellung->value("KeepMapIsoValues").toBool());
      // keepiso->setChecked(mole->einstellung->value("KeepMapIsoValues").toBool());
    }
    mole->einstellung->endGroup();    
    schnittart(maptrunc);
    //    printf("@!#$  %p\n",mole->einstellung);
    md->hide();
  }
  connect(ownfcf6,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
  isBelo=false;
  //printf("constr. belo?%d\n",isBelo);
  //printf("%s\n",ENUS);
}

FourXle::~FourXle(){
    //printf("Maps will close...\n");
  if (chgl!=NULL) disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
  if (chgl!=NULL) disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
  if (chgl!=NULL) disconnect(chgl,SIGNAL(inimibas()),0,0);
  //printf("%d\n",__LINE__);
  killmaps();
  //printf("%d\n",__LINE__);
  delete doMaps; 
  //printf("Maps closed.\n");
}

void FourXle::mapDefault(){
  map_radius= 5.0;
  maptrunc = 2;
  rr= 2.0;
  rw= 1.0;
  chgl->lintrans= 0.8;
  chgl->linwidth= 0.5;
  iso[1]= sigma[0] * 2.7f;
  iso[0]= sigma[1] * 1.2f;
  chgl->niceTrans->setChecked(false);
  chgl->lighting->setChecked(true);
  chgl->fillMap->setChecked(false);
  weak->setValue(rw);
  mapprec->setValue(rr);
  difmaps->setValue(fabs(iso[1]));
  difmaps->setSingleStep(fabs(iso[1])/10.0);
  lineTrans->setValue(chgl->lintrans);
  lineWidth->setValue(chgl->linwidth);
  fomaps->setValue(iso[0]);
  fomaps->setSingleStep(iso[0]/10.0);
  maprad->setValue(map_radius);
  mapSchnitt->setCurrentIndex(maptrunc);
  controlMap();
}

void FourXle::schnittart(int trunc){
  if (trunc==1) maprad->show();
  else maprad->hide();
  if (trunc==2) cutrange->show();
  else cutrange->hide();
}

void FourXle::keepIsoValues(bool b){
  keepiso->setChecked(b);
  keepIso->setChecked(b);
  keepIso->setVisible(b);
  if (n5>0){
    sigma[0] = sigmaFoFc();
    sigma[1]  = sigmaFo();
    iso[1]= sigma[0] * 2.7f;
    iso[0]= sigma[1] * 1.2f;
    fomaps->setValue(iso[0]);
    difmaps->setValue(fabs(iso[1]));
    controlMap();
  }
}

void FourXle::updateLines(){
  chgl->lintrans=lineTrans->value();
  chgl->linwidth=lineWidth->value();
  chgl->updateGL();
}

void FourXle::killFourXleFirst(){
  //fprintf(stderr,"chgl dies...argh...\n");
  killmaps(); 
  n5=0;
  if (mapcontrol!=NULL) mapcontrol->setVisible(false);
}

void FourXle::doMapsNow(bool b){
  //isBelo=false;
  bool bb=MYFCF;
  bool B = (MYFCF!=ownfcf6->isChecked())||(MYFCF!=ownfcf6_->isChecked());
  MYFCF=(B)?!MYFCF:MYFCF;
  if (bb!=MYFCF) {
    ownfcf6->disconnect();
    ownfcf6_->disconnect();
    ownfcf6->setChecked(MYFCF);
    ownfcf6_->setChecked(MYFCF);
    connect(ownfcf6,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
    connect(ownfcf6_,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
    emit reloadRes();
    return;
  } 
  if (!b) {
    killmaps(); 
    n5=0;
    if (mapcontrol!=NULL) mapcontrol->setVisible(false);
    chgl->updateGL();
    return;
  } else {
    bool hamrsho=((iso[0]==iso[1])&&((iso[0]==0.0f)||(iso[0]==0.03f)));
    if (!loadFouAndPerform(fouName.toLocal8Bit(),hamrsho)){
      emit bigmessage(QString("<font color=red>Error while loading %1!</font><br>").arg(fouName));      
      if (mapcontrol!=NULL) mapcontrol->setVisible(false);
    } else {        
      if (mapcontrol!=NULL) mapcontrol->setVisible(true);
      if (!hamrsho) {
        sigma[0] = sigmaFoFc();
        sigma[1]  = sigmaFo();
        iso[1]= sigma[0] * 2.7f;
        iso[0]= sigma[1] * 1.2f;
        fomaps->setValue(iso[0]);
        difmaps->setValue(fabs(iso[1]));
        inimap();
      }
    }
    chgl->updateGL();
  }

} 

void FourXle::controlMap(){
  maptrunc=mapSchnitt->currentIndex (); 
  if ((rw!=weak->value())||(rr!=mapprec->value())){
    rr =  mapprec->value();  
    rw = weak->value();
    if (!loadFouAndPerform(fouName.toLocal8Bit(),false)){
      emit bigmessage(QString("<font color=red>Error while loading %1!</font><br>").arg(fouName));      
    }
  }    
  chgl->lintrans = lineTrans->value();
  chgl->linwidth = lineWidth->value();
  map_radius = maprad->value();
  iso[0] = fomaps->value();
  iso[1] = -difmaps->value();
  inimap();
  chgl->updateGL();
}

void FourXle::killmaps(){
  //printf("KM belo?%d chgl->%p dataFo-Fc->%p dataFo->%p\n",isBelo,chgl,datfo_fc,datfo);
  isBelo=false;
 // printf("These nodes will be killed %p %p %p\n",nodex,nodey,nodez);
  /*! deletes all fourier maps and frees the memory
  */
  // printf("I kill maps\n");

  if (chgl!=NULL) disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
  if (chgl!=NULL) disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
  if (chgl!=NULL) disconnect(chgl,SIGNAL(inimibas()),0,0);  
  //printf("emit halt\n");
  //emit haltMW();
  //printf("halt emited\n");
  if ((chgl!=NULL)&&(chgl->fofcact!=NULL)) chgl->fofcact->setVisible(false);
  if ((chgl!=NULL)&&(chgl->foact!=NULL)) chgl->foact->setVisible(false);
  if (datfo!=NULL) free(datfo);
  if (datfo_fc!=NULL) free(datfo_fc);
  if (nodex!=NULL) free(nodex);
  nodex=NULL;
  if (nodey!=NULL) free(nodey);
  nodey=NULL;
  if (nodez!=NULL) free(nodez);
  nodez=NULL;
  if (noGoMap!=NULL) free(noGoMap);
  noGoMap=NULL;
  if (nodeX!=NULL) free(nodeX);
  if (nodeY!=NULL) free(nodeY);
  if (nodeZ!=NULL) free(nodeZ);
  nodeX=nodeY=nodeZ=NULL;
  datfo=datfo_fc=NULL;
  if (chgl!=NULL){
     //printf("I kill maps %d! %d!\n",chgl->fVertexes.size(),chgl->fNormals.size());
     int sumsize=0;
     for (int i=0; i<chgl->fVertexes.size(); i++) {sumsize+=chgl->fVertexes[i].size();chgl->fVertexes[i].clear();}
     for (int i=0; i<chgl->fNormals.size(); i++) {sumsize+=chgl->fNormals[i].size();chgl->fNormals[i].clear();}
     //printf("Freed %d normals and vertices.\n",sumsize);
  }//else printf("chgl is already dead!\n");
}

void FourXle::colorDLGFOP(){
  QColor color;
  color=QColorDialog::getColor(chgl->fopc, md);
  if (color.isValid()) {
    chgl->fopc=color;
    foPlusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg(color.name())
        .arg(color.darker(200).name())
        .arg(color.lighter(100).name())
        .arg("#ffffff"));
  }
  mole->einstellung->beginGroup("Colors");
  mole->einstellung->setValue("FourFoPlusColor",   chgl->fopc);
  mole->einstellung->setValue("FourFoMinusColor",  chgl->fomc);
  mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
  mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
  mole->einstellung->endGroup();
  md->update();

}

void FourXle::colorDLGFOM(){
  QColor color;
  color=QColorDialog::getColor(chgl->fopc, md);
  if (color.isValid()) {
    chgl->fomc=color;
    foMinusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg(color.name())
        .arg(color.darker(200).name())
        .arg(color.lighter(100).name())
        .arg("#ffffff"));
  }
  mole->einstellung->beginGroup("Colors");
  mole->einstellung->setValue("FourFoPlusColor",   chgl->fopc);
  mole->einstellung->setValue("FourFoMinusColor",  chgl->fomc);
  mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
  mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
  mole->einstellung->endGroup();
  md->update();

}

void FourXle::colorDLGDIP(){
  QColor color;
  color=QColorDialog::getColor(chgl->fopc, md);
  if (color.isValid()) {
    chgl->dipc=color;

    diffPlusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg(color.name())
        .arg(color.darker(200).name())
        .arg(color.lighter(100).name())
        .arg("#ffffff"));
  }
  mole->einstellung->beginGroup("Colors");
  mole->einstellung->setValue("FourFoPlusColor",   chgl->fopc);
  mole->einstellung->setValue("FourFoMinusColor",  chgl->fomc);
  mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
  mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
  mole->einstellung->endGroup();
  md->update();
}

void FourXle::colorDLGDIM(){
  QColor color;
  color=QColorDialog::getColor(chgl->fopc, md);
  if (color.isValid()) {
    chgl->dimc=color;

    diffMinusButton->setStyleSheet(QString(
          "QPushButton {"
          "font-size:16px;"
          "font-weight:bold;"
          "border: 1px solid #000000;"
          "border-radius: 9px;"
          "color: %4;"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
          "}"
          "QPushButton:hover {"
          "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
          "}"
          "QPushButton:flat {"
          "    border: none; /* no border for a flat push button */"
          "}"
          )
        .arg(color.name())
        .arg(color.darker(200).name())
        .arg(color.lighter(100).name())
        .arg("#ffffff"));
  }
  mole->einstellung->beginGroup("Colors");
  mole->einstellung->setValue("FourFoPlusColor",   chgl->fopc);
  mole->einstellung->setValue("FourFoMinusColor",  chgl->fomc);
  mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
  mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
  mole->einstellung->endGroup();
  md->update();

}

void FourXle::openMapControl(){
  weak->setValue(rw);
  mapprec->setValue(rr);
  difmaps->setValue(fabs(iso[1]));
  difmaps->setSingleStep(fabs(iso[1])/10.0);
  lineTrans->setValue(chgl->lintrans);
  lineWidth->setValue(chgl->linwidth);
  fomaps->setValue(iso[0]);
  fomaps->setSingleStep(iso[0]/10.0);
  maprad->setValue(map_radius);
  mapSchnitt->setCurrentIndex(maptrunc);
  md->show();
}

void FourXle::defaultColors(){
  chgl->fopc = QColor("#0000FF");
  chgl->fomc = QColor("#FFA600");
  chgl->dipc = QColor("#00B200");
  chgl->dimc = QColor("#CC0000");
  foPlusButton->setStyleSheet(QString(
        "QPushButton {"
        "font-size:16px;"
        "font-weight:bold;"
        "border: 1px solid #000000;"
        "border-radius: 9px;"
        "color: %4;"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
        "}"
        "QPushButton:hover {"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
        "border: 0px"
        "}"
        "QPushButton:flat {"
        "    border: none; /* no border for a flat push button */"
        "}"
        )
      .arg( chgl->fopc.name())
      .arg( chgl->fopc.darker(200).name())
      .arg( chgl->fopc.lighter(100).name())
      .arg("#ffffff"));

  foMinusButton->setStyleSheet(QString(
        "QPushButton {"
        "font-size:16px;"
        "font-weight:bold;"
        "border: 1px solid #000000;"
        "border-radius: 9px;"
        "color: %4;"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
        "}"
        "QPushButton:hover {"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
        "border: 0px"
        "}"
        "QPushButton:flat {"
        "    border: none; /* no border for a flat push button */"
        "}"
        )
      .arg( chgl->fomc.name())
      .arg( chgl->fomc.darker(200).name())
      .arg( chgl->fomc.lighter(100).name())
      .arg("#ffffff"));

  diffMinusButton->setStyleSheet(QString(
        "QPushButton {"
        "font-size:16px;"
        "font-weight:bold;"
        "border: 1px solid #000000;"
        "border-radius: 9px;"
        "color: %4;"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
        "}"
        "QPushButton:hover {"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
        "border: 0px"
        "}"
        "QPushButton:flat {"
        "    border: none; /* no border for a flat push button */"
        "}"
        )
      .arg( chgl->dimc.name())
      .arg( chgl->dimc.darker(200).name())
      .arg( chgl->dimc.lighter(100).name())
      .arg("#ffffff"));

  diffPlusButton->setStyleSheet(QString(
        "QPushButton {"
        "font-size:16px;"
        "font-weight:bold;"
        "border: 1px solid #000000;"
        "border-radius: 9px;"
        "color: %4;"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
        "}"
        "QPushButton:hover {"
        "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
        "border: 0px"
        "}"
        "QPushButton:flat {"
        "    border: none; /* no border for a flat push button */"
        "}"
        )
      .arg( chgl->dipc.name())
      .arg( chgl->dipc.darker(200).name())
      .arg( chgl->dipc.lighter(100).name())
      .arg("#ffffff"));
  mole->einstellung->beginGroup("Colors");
  mole->einstellung->setValue("FourFoPlusColor",   chgl->fopc);
  mole->einstellung->setValue("FourFoMinusColor",  chgl->fomc);
  mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
  mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
  mole->einstellung->endGroup();
  md->update();

} 

bool FourXle::loadFouAndPerform(const char filename[],bool neu){
  /*! loads a fcf file prepares the reciprocal data for the fourier transpormation and performs it.
  */
  foti.start();

  const int it[143]= {2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,45,48,50,54,60,64,72,75,80,81,90,96,100,
    108,120,125,128,135,144,150,160,162,180,192,200,216,225,240,243,250,256,270,288,300,320,324,360,375,384,400,405,
    432,450,480,486,500,512,540,576,600,625,640,648,675,720,729,750,768,800,810,864,900,960,972,1000,1024,1080,1125,
    1152,1200,1215,1250,1280,1296,1350,1440,1458,1500,1536,1600,1620,1728,1800,1875,1920,1944,2000,2025,2048,2160,
    2187,2250,2304,2400,2430,2500,2560,2592,2700,2880,2916,3000,3072,3125,3200,3240,3375,3456,3600,3645,3750,3840,
    3888,4000,4050,4096,4320,4374,4500,4608,4800,4860,5000};//!multiples of 2 3 and 5 
  mole->einstellung->beginGroup("Colors");
  chgl->fopc=mole->einstellung->value("FourFoPlusColor",   QColor("#0000FF")).value<QColor>();
  chgl->fomc=mole->einstellung->value("FourFoMinusColor",  QColor("#FFA600")).value<QColor>();
  chgl->dipc=mole->einstellung->value("FourDiffPlusColor", QColor("#00B200")).value<QColor>();
  chgl->dimc=mole->einstellung->value("FourDiffMinusColor",QColor("#CC0000")).value<QColor>();
  mole->einstellung->endGroup();
  //printf("%d\n",__LINE__);
  killmaps();
  //printf("%d {%s}\n",__LINE__,filename);
  int ok;
  if (!doMaps->isChecked()) return false;
  if (strstr(filename,".fcf")==NULL) return false;
  FILE *f;
  //printf("%d\n",__LINE__);
  ok= readHeader(filename);
  //printf("%d\n",__LINE__);
  if (ok) {
    switch (ok){
      case 1: emit bigmessage("<font color=red>Map generation failed. SHELXL LIST code was not 6.</font>");break;
      case 2: emit bigmessage(QString("<font color=red>Map generation failed. File %1 corrupted.</font>").arg(QString::fromLocal8Bit(filename)));break;
      case 3: emit bigmessage(QString("<font color=red>Map generation failed. Cannot open file %1.</font>").arg(QString::fromLocal8Bit(filename)));break;
      case 4: emit bigmessage(QString("<font color=red>Map generation failed. No reflection data in file %1.</font>").arg(QString::fromLocal8Bit(filename)));break;
    }
    return false;
  }
  //printf("%d\n",__LINE__);
  f=fopen(filename,"rb");
  if (f==NULL)return false;
  char line[122]="";
  while (strstr(line,"_refln_phase_calc")==NULL) {
    if (fgets(line,120,f)){;}//just a stupid supression for annoying warnig
  }
  nr=0;
  lr[nr].ih=0;
  lr[nr].ik=0;
  lr[nr].il=0; 
  lr[nr].fo=f000*f000;
  lr[nr].so=0.5f;
  lr[nr].fc=f000;
  lr[nr].phi=0.0f;
  nr=1;
  //printf("f000 %g\n",f000);
  emit bigmessage(QString::fromLocal8Bit(filename));

  while (!feof(f)){
    if (fgets(line,120,f)){;}//just a stupid supression for annoying warnig
    int rdi=
      sscanf(line,"%d %d %d %f %f %f %f",&lr[nr].ih,&lr[nr].ik, &lr[nr].il ,&lr[nr].fo, &lr[nr].so, &lr[nr].fc, &lr[nr].phi);
    if (rdi==7) {
      if ((abs(lr[nr].ih)<HKLMX)&&
          (abs(lr[nr].ik)<HKLMX)&&
          (abs(lr[nr].il)<HKLMX)&&
          ((lr[nr].ih|lr[nr].ik|lr[nr].il)!=0))

        //printf("%4d%4d%4d %g %g\n", lr[nr].ih, lr[nr].ik, lr[nr].il ,lr[nr].fo, lr[nr].so);
        nr++;
    }

    // else printf("?? %d {%s}\n",rdi,line);

    if (nr>=LM) {
      emit bigmessage(QString("<font color=red>to many reflections in fcf file</font>"));
      return false;
    }
  }
  fclose(f);
  if (nr<2) {
    emit bigmessage(QString("<font color=red>Map generation failed. No reflection data in file %1.</font>").arg(QString::fromLocal8Bit(filename)));
    return false;
  }

  for (int i=0;i<nr;i++){
    double u=lr[i].ih,v=lr[i].ik,w=lr[i].il;
    int mh=lr[i].ih,mk=lr[i].ik,ml=lr[i].il;
    double p,q=lr[i].phi/180.0*M_PI;
    lr[i].phi=fmod(4*M_PI+q,2*M_PI);
    for (int k=0; k<ns; k++){
      int nh,nk,nl;
      double t=1.0;
      nh=(int) (u*sy[0][k]+ v*sy[3][k] + w*sy[6][k]);
      nk=(int) (u*sy[1][k]+ v*sy[4][k] + w*sy[7][k]);
      nl=(int) (u*sy[2][k]+ v*sy[5][k] + w*sy[8][k]);
      if((nl<0)||((nl==0)&&(nk<0))||((nl==0)&&(nk==0)&&(nh<0)))
      {nh*=-1;nk*=-1;nl*=-1;t=-1.0;}
      if ((nl<ml)||((nl==ml)&&(nk<mk))||((nl==ml)&&(nk==mk)&&(nh<=mh))) continue;
      mh=nh;mk=nk;ml=nl;
      p=u*sy[9][k]+v*sy[10][k]+w*sy[11][k];
      lr[i].phi=fmod(4*M_PI+t*fmod(q-2*M_PI*p,2*M_PI)-0.01,2*M_PI)+0.01;

    }
    lr[i].ih=mh;
    lr[i].ik=mk;
    lr[i].il=ml;
  }
  sorthkl(nr,lr);
  int n=-1;
  {int i=0;
    while(i<nr){
      double t=0.;
      double u=0.;
      double v=0.;
      double z=0.;
      double p=0.;
      int m;
      int k=i;
      while ((i<nr)&&(lr[i].ih==lr[k].ih)&&(lr[i].ik==lr[k].ik)&&(lr[i].il==lr[k].il)) {
        t=t+1.;
        u+=lr[i].fo;
        v+=1./(lr[i].so*lr[i].so);
        z+=lr[i].fc;
        p=lr[i].phi;
        i++;
      }
      m=n+1;
      lr[m].fo=sqrt(fmax(0.,u/t));
      lr[m].so=sqrt(lr[m].fo*lr[m].fo+sqrt(1./v))-lr[m].fo;
      lr[m].fc=z/t;
      lr[m].phi=p;
      n=m;
      lr[n].ih=lr[k].ih;
      lr[n].ik=lr[k].ik;
      lr[n].il=lr[k].il;
    }
  }
  n++;
  nr=n;
  //  printf("%4d%4d%4d %g %g %g %g\n",lr[0].ih,lr[0].ik,lr[0].il,lr[0].fo,lr[0].so,lr[0].fc,lr[0].phi);
  {
    float DX;
    float DY;
    float DZ;	


    {
      int mh=0, mk=0, ml=0,j;
      for (int n=0; n<nr; n++){
        double u=lr[n].ih,v=lr[n].ik,w=lr[n].il;
        int a,b,c;
        for (int k=0; k<ns;k++){
          a=abs((int)(u*sy[0][k]+v*sy[3][k]+w*sy[6][k]));
          b=abs((int)(u*sy[1][k]+v*sy[4][k]+w*sy[7][k]));
          c=abs((int)(u*sy[2][k]+v*sy[5][k]+w*sy[8][k]));
          mh=(mh<a)?a:mh;
          mk=(mk<b)?b:mk;
          ml=(ml<c)?c:ml;
        }
      }
      j=(int)(rr*mh+.5);
      for (int i=0; it[i]< j; i++)n1=it[i+1];
      j=(int)(rr*mk+.5);
      for (int i=0; it[i]< j; i++)n2=it[i+1];
      j=(int)(rr*ml+.5);
      for (int i=0; (it[i]< j)||((nc)&&(it[i]%2)); i++) n3=it[i+1];
      if (!voxelstr.isEmpty()) {
        if (voxelstr.contains('x')) {
          QStringList lll=voxelstr.split('x');
          if (lll.size()==3){
            n1=lll.at(0).toInt();
            n2=lll.at(1).toInt();
            n3=lll.at(2).toInt();
            //printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n %d %d %d \n",n1,n2,n3);
          }
        }
      }
      n4=n2*n1;
      n5=n3*n4;
      datfo=(float*) malloc(sizeof(float)*n5);
      if (datfo==NULL){
          fprintf(stderr,"Too less memory for F_observed!\n");
          datfo=NULL;
          mapDefault();
          info="";
          emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
          return false;
      }
      datfo_fc=(float*) malloc(sizeof(float)*n5);
      if (datfo_fc==NULL){
          fprintf(stderr,"Too less memory for F_observed - F_calculated!\n");
          free(datfo);
          datfo=NULL;
          datfo_fc=NULL;          
          mapDefault();
          info="";
          emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
          return false;
      }
      DX=1.0f/n1;
      DY=1.0f/n2;
      DZ=1.0f/n3;
    } 
    for (int typ=0; typ<2;typ++){
      double miZ=99999.99,maZ=-99999.99;
      //      mInimum[typ]=99999.99; mAximum[typ]=-99999.99;
      int nbytes,dims[3];
      dims[0]=n3;
      dims[1]=n2;
      dims[2]=n1;
#ifdef FFTW3_H
      B=(fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*n5);
      if (B==NULL)return false; 
      for (int i=0; i<n5; i++){B[i][0]=0;B[i][1]=0;}
#else
      B=(kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes = (sizeof(kiss_fft_cpx)*n5));
      if (B==NULL){
          fprintf(stderr,"Too less memory for FFT map!\n");
          free(datfo);
          free(datfo_fc);
          datfo=NULL;
          datfo_fc=NULL;
          mapDefault();
          info="";
          emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
          return false;
      }
      for (int i=0; i<n5; i++){B[i].r=0;B[i].i=0;}
#endif
      for (int i=0; i<nr;i++){
        float  u,v,w;
        u=lr[i].ih;
        v=lr[i].ik;
        w=lr[i].il;
        float  ss,s=0,t=0,q,p;
        for (int n=0; n<ns;n++){
          int j,k,l;
          j=(int) (u*sy[0][n]+ v*sy[3][n] + w*sy[6][n]);
          k=(int) (u*sy[1][n]+ v*sy[4][n] + w*sy[7][n]);
          l=(int) (u*sy[2][n]+ v*sy[5][n] + w*sy[8][n]);
          if((abs(j-lr[i].ih)+abs(k-lr[i].ik)+abs(l-lr[i].il))==0)s+=1.0;
          if(abs(j+lr[i].ih)+abs(k+lr[i].ik)+abs(l+lr[i].il)==0)t+=1.0;
        }
        if (i==0) {//printf("v%f s%f t%f\n",C[14],s,t);
          s=1;t=0;//f000 
        }
        if(typ==0) ss=(lr[i].fo-lr[i].fc)/(C[14]*(s+t));
        else if(typ==2) //here this happens never
          ss=(lr[i].fc)/(C[14]*(s+t));
        else ss=(lr[i].fo)/(C[14]*(s+t));
        if(lr[i].fc>1.E-6) ss=ss/(1.+rw*pow(lr[i].so/lr[i].fc,4));
        for (int n=0; n<ns;n++){
          int j,k,l,m;
          j=(int) (u*sy[0][n]+ v*sy[3][n] + w*sy[6][n]);
          k=(int) (u*sy[1][n]+ v*sy[4][n] + w*sy[7][n]);
          l=(int) (u*sy[2][n]+ v*sy[5][n] + w*sy[8][n]);
          //          q=(-2*M_PI*(u*sy[9][n]+v*sy[10][n]+w*sy[11][n]))-M_PI*(j*DX+k*DY+l*DZ);
          q=(lr[i].phi-2*M_PI*(u*sy[9][n]+v*sy[10][n]+w*sy[11][n]))-M_PI*(j*DX+k*DY+l*DZ);
          j=(999*n1+j)%n1;
          k=(999*n2+k)%n2;
          l=(999*n3+l)%n3;
          m=j+n1*(k+n2*l);
          p=ss*cosf(q);
          q=ss*sinf(q);
#ifdef FFTW3_H
          B[m][0]=p;
          B[m][1]=q;
#else
          B[m].r=p;
          B[m].i=q;
#endif
          j*=-1;
          if(j<0)j=n1+j;
          k*=-1;
          if(k<0)k=n2+k;
          l*=-1;
          if(l<0)l=n3+l;
          m=j+n1*(k+n2*l);
#ifdef FFTW3_H
          B[m][0]=p;
          B[m][1]=-q;
#else
          B[m].r=p;
          B[m].i=-q;
#endif
        }
      }
#ifdef FFTW3_H
      fwd_plan = fftwf_plan_dft_3d(n3,n2,n1,B,B,FFTW_FORWARD,FFTW_ESTIMATE);
      fftwf_execute(fwd_plan);
      fftwf_destroy_plan(fwd_plan);
#else
      fwd_plan = kiss_fftnd_alloc(dims,3,0,0,0);
      kiss_fftnd( fwd_plan,B,B);
      free(fwd_plan);
#endif
      float t=0;
      double DM=0.,  DS=0., DD  ;
      for (int i=0; i<n5;i++){
#ifdef FFTW3_H
        DD=B[i][0];
#else
        DD=B[i].r;
#endif
        //	maxi[typ]=fmax(maxi[typ],DD); mini[typ]=fmin(mini[typ],DD);
        miZ=fmin(miZ,DD);
        maZ=fmax(maZ,DD);
        DM+=DD;
        DS+=DD*DD;
#ifdef FFTW3_H
        if (typ==1) datfo[i]=B[i][0];
        else if (typ==0) datfo_fc[i]=B[i][0];
#else
        if (typ==1) datfo[i]=B[i].r;
        else if (typ==0) datfo_fc[i]=B[i].r;

#endif
      }
      sigma[typ]=t=sqrt((DS/n5)-((DM/n5)*(DM/n5)));
      free(B);
    }//1
  }//2
  urs=V3(0,0,0);int gt=0;
  for (int i=0; i<mole->showatoms.size();i++)
    if (mole->showatoms.at(i).an>-1) {urs+=mole->showatoms.at(i).frac;gt++;}
  urs*=1.0/gt;
  urs=V3(1,1,1)-1.0*urs;
  mole->frac2kart(urs,urs);

  nodex= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodex==NULL){
      fprintf(stderr,"Too less memory for node X!\n");
      free(datfo);
      free(datfo_fc);
      datfo=NULL;
      datfo_fc=NULL;
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  nodey= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodey==NULL){
      fprintf(stderr,"Too less memory for node Y!\n");
      free(nodex);
      free(datfo);
      free(datfo_fc);
      nodex = NULL;
      datfo = NULL;
      datfo_fc = NULL;
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  nodez= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodez==NULL){
      fprintf(stderr,"Too less memory for node Z!\n");
      free(nodex);
      free(nodey);
      free(datfo);
      free(datfo_fc);
      nodex = NULL;
      nodey = NULL;
      datfo = NULL;
      datfo_fc = NULL;
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  noGoMap= (int*) malloc(sizeof(int)*n5*27);
  if (noGoMap==NULL){
      fprintf(stderr,"Too less memory for noGoMap!\n");

      free(nodex);
      free(nodey);
      free(nodez);
      free(datfo);
      free(datfo_fc);
      nodex= NULL;
      nodey= NULL;
      nodez= NULL;
      datfo= NULL;
      datfo_fc= NULL;
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  for (int no=0;no<(n5*27);no++) noGoMap[no]=0; 
  for (int o=0; o<n5;o++){
    nodex[o].flag=0;
    nodey[o].flag=0;
    nodez[o].flag=0;
  }
  dx=V3(1.0/(n1),0,0);
  dy=V3(0,1.0/(n2),0);
  dz=V3(0,0,1.0/(n3));
  mole->frac2kart(dx,dx); 
  mole->frac2kart(dy,dy); 
  mole->frac2kart(dz,dz);     
  delDA[ 0]=  -n1*dx -n2*dy -n3*dz;//nx ny,nz??
  delDA[ 1]=         -n2*dy -n3*dz;
  delDA[ 2]=   n1*dx -n2*dy -n3*dz;
  delDA[ 3]=  -n1*dx        -n3*dz;
  delDA[ 4]=                -n3*dz;
  delDA[ 5]=   n1*dx        -n3*dz;
  delDA[ 6]=  -n1*dx +n2*dy -n3*dz;
  delDA[ 7]=          n2*dy -n3*dz;
  delDA[ 8]=   n1*dx +n2*dy -n3*dz;
  delDA[ 9]=  -n1*dx -n2*dy       ;
  delDA[10]=         -n2*dy       ;
  delDA[11]=   n1*dx -n2*dy       ;
  delDA[12]=  -n1*dx              ;
  delDA[13]=  V3(0,0,0)           ;
  delDA[14]=   n1*dx              ;
  delDA[15]=  -n1*dx +n2*dy       ;
  delDA[16]=          n2*dy       ;
  delDA[17]=   n1*dx +n2*dy       ;
  delDA[18]=  -n1*dx -n2*dy +n3*dz;
  delDA[19]=         -n2*dy +n3*dz;
  delDA[20]=   n1*dx -n2*dy +n3*dz;
  delDA[21]=  -n1*dx        +n3*dz;
  delDA[22]=                +n3*dz;
  delDA[23]=   n1*dx        +n3*dz;
  delDA[24]=  -n1*dx +n2*dy +n3*dz;
  delDA[25]=          n2*dy +n3*dz;
  delDA[26]=   n1*dx +n2*dy +n3*dz;
  emit bigmessage(QString("Uniq Reflections: %1 <br>Fourier grid dimensions: %2 %3 %4 = %5 grid points.<br>Time for loading and fourier transformation: <b>%6 s</b>.<br>").arg(nr).arg(n1).arg(n2).arg(n3).arg(n5).arg(foti.elapsed()/1000.0,1,'f',1));
  if (neu) gen_surface(neu);
  return true;
}

bool FourXle::loadMap(int NX, int NY, int NZ, float *dat){
  bool wbi=isBelo;
  //printf("%d\n",__LINE__);
  killmaps();
  //printf("%d\n",__LINE__);
  isBelo=wbi;

  n1=NX;
  n2=NY;
  n3=NZ;
  n4=n2*n1;
  n5=n3*n4;
  datfo=(float*) malloc(sizeof(float)*n5);
  if (datfo==NULL){
      fprintf(stderr,"Too less memory for F_observed - F_calculated!\n");
      free(datfo_fc);
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  for (int i=0; i<n5; i++)datfo[i]=0.0f;
  datfo_fc=dat;

  sigma[0] = sigmaFoFc();
  sigma[1] = sigmaFo();
  iso[1]= sigma[0] * 2.5f;
  iso[0]= sigma[1] * 1.2f;
  fomaps->setValue(iso[0]);
  difmaps->setValue(iso[1]);

  urs=V3(0,0,0);int gt=0;
  for (int i=0; i<mole->showatoms.size();i++)
    if (mole->showatoms.at(i).an>-1) {urs+=mole->showatoms.at(i).frac;gt++;}
  urs*=1.0/gt;
  urs=V3(1,1,1)-1.0*urs;
  mole->frac2kart(urs,urs);
  nodex= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodex==NULL){
      fprintf(stderr,"Too less memory for node X!\n");
      free(datfo);
      free(datfo_fc);
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  nodey= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodey==NULL){
      fprintf(stderr,"Too less memory for node Y!\n");
      free(datfo);
      free(datfo_fc);
      free(nodex);
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  nodez= (FNode*)malloc(sizeof(FNode)*n5);
  if (nodez==NULL){
      fprintf(stderr,"Too less memory for node Z!\n");
      free(datfo);
      free(datfo_fc);
      free(nodex);
      free(nodey);
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  noGoMap= (int*) malloc(sizeof(int)*n5*27);
  if (noGoMap==NULL){
      fprintf(stderr,"Too less memory for noGoMap!\n");
      free(datfo);
      free(datfo_fc);
      free(nodex);
      free(nodey);
      free(nodez);
      mapDefault();
      info="";
      emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
      return false;
  }
  for (int no=0;no<(n5*27);no++) noGoMap[no]=0;
  for (int o=0; o<n5;o++){
    nodex[o].flag=0;
    nodey[o].flag=0;
    nodez[o].flag=0;
  }
  dx = V3(1.0 / (n1), 0, 0);
  dy = V3(0, 1.0 / (n2), 0);
  dz = V3(0, 0, 1.0 / (n3));
  mole->frac2kart(dx,dx);
  mole->frac2kart(dy,dy);
  mole->frac2kart(dz,dz);
  delDA[ 0]=  -n1*dx -n2*dy -n3*dz;//nx ny,nz??
  delDA[ 1]=         -n2*dy -n3*dz;
  delDA[ 2]=   n1*dx -n2*dy -n3*dz;
  delDA[ 3]=  -n1*dx        -n3*dz;
  delDA[ 4]=                -n3*dz;
  delDA[ 5]=   n1*dx        -n3*dz;
  delDA[ 6]=  -n1*dx +n2*dy -n3*dz;
  delDA[ 7]=          n2*dy -n3*dz;
  delDA[ 8]=   n1*dx +n2*dy -n3*dz;
  delDA[ 9]=  -n1*dx -n2*dy       ;
  delDA[10]=         -n2*dy       ;
  delDA[11]=   n1*dx -n2*dy       ;
  delDA[12]=  -n1*dx              ;
  delDA[13]=  V3(0,0,0)           ;
  delDA[14]=   n1*dx              ;
  delDA[15]=  -n1*dx +n2*dy       ;
  delDA[16]=          n2*dy       ;
  delDA[17]=   n1*dx +n2*dy       ;
  delDA[18]=  -n1*dx -n2*dy +n3*dz;
  delDA[19]=         -n2*dy +n3*dz;
  delDA[20]=   n1*dx -n2*dy +n3*dz;
  delDA[21]=  -n1*dx        +n3*dz;
  delDA[22]=                +n3*dz;
  delDA[23]=   n1*dx        +n3*dz;
  delDA[24]=  -n1*dx +n2*dy +n3*dz;
  delDA[25]=          n2*dy +n3*dz;
  delDA[26]=   n1*dx +n2*dy +n3*dz;
  rr =  mapprec->value();
  rw = weak->value();
  //gen_surface_silent();
  return true;

}


void FourXle::trimm(char s[]){
  /*! a trimm function for c-strings.
  */
  char sc[409];
  int j=0;
  int len=strlen(s);
  strncpy(sc,s,400);
  for (int i=0; i<len; i++) if ((sc[i]!='\'')&&(!isspace(sc[i]))) s[j++]=toupper(sc[i]);
  s[j]='\0';
}

void FourXle::deletes(char *s, int count){
  /*! deletes count characters at the begining of s.
  */
  if ((s==NULL)||(count <1)||((size_t)count>strlen(s))) return;
  for (int i=0; i<count; i++) s[i]=' ';
  trimm(s);
}

void FourXle::exportMaps(int na, const char filename[], const char atomlist[]){
  if ((datfo==NULL)||(datfo_fc==NULL)) return;
  FILE *fo,*fof1;//,*f1f2;
  char foname[4096];
  char fof1name[4096];
  //  char f1f2name[4096];
  int len=strlen(filename);
  float factor=0.1481847095290449f;//a0**3
  double a0=0.52917720859;  //bohr=0.5291775108
  int i=0;
  //
  //FO MAP
  //
  V3 dx1=V3(1.0/(n1),0,0);
  V3 dy1=V3(0,1.0/(n2),0);
  V3 dz1=V3(0,0,1.0/(n3));
  mole->frac2kart(dx1,dx1); 
  mole->frac2kart(dy1,dy1); 
  mole->frac2kart(dz1,dz1);     

  strncpy(foname,filename,len-4);
  foname[len-4]='\0';
  strcat(foname,"_fo_densitymap.cube");
  fo=fopen(foname,"w");
  if (fo==NULL) return ;
  double fomax=-9e37,fomin=9e37,em=0,ep=0,et=0;
  for (int i=0; i<n5; i++ ){
    if (datfo[i]>0.0) ep+=datfo[i];
    else em+=datfo[i];
    et+=datfo[i];
    fomax=fmax(datfo[i],fomax);
    fomin=fmin(datfo[i],fomin);
  }
  //printf("Fo %g %g %g %g %g\n",fomin,fomax,ep/n5*C[14],em/n5*C[14],et/n5*C[14] );
  fprintf(fo,"F observed map written by ShelXle\nDensity obtained from Fo in fcf with phases of model transformed using fft\n");
  fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",na,
      (dx1.x+dy1.x+dz1.x)*0.5/a0,
      (dx1.y+dy1.y+dz1.y)*0.5/a0,
      (dx1.z+dy1.z+dz1.z)*0.5/a0
      //0.0,0.0,0.0
      );
  fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",n1,dx1.x/a0,dx1.y/a0,dx1.z/a0);
  fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",n2,dy1.x/a0,dy1.y/a0,dy1.z/a0);
  fprintf(fo,"%5d%12.6f%12.6f%12.6f"  ,n3,dz1.x/a0,dz1.y/a0,dz1.z/a0);//no newline here because it is in atomlist
  fprintf(fo,"%s",atomlist);
  for (int xi=0;xi<n1;xi++)
    for (int yi=0;yi<n2;yi++)
      for (int zi=0;zi<n3;zi++)
      {
        //fprintf(fo,"%s%13.5E",(((i%6)==0)?"\n":""),datfo[dex(xi,yi,zi)]*factor);
        fprintf(fo,"%s",QString("%1%2").arg(((i%6)==0)?"\n":"").arg(datfo[dex(xi,yi,zi)]*factor,13,'E',5).toStdString().c_str());
        i++;
      }
  fclose(fo);
  //
  //FO-F1 MAP
  //

  strncpy(fof1name,filename,len-4);
  fof1name[len-4]='\0';
  strcat(fof1name,"_fo-fc_densitymap.cube");
  fof1=fopen(fof1name,"w");
  if (fof1==NULL) return ;
  fomax=-9e37;
  fomin= 9e37;
  em=0;
  ep=0;
  et=0;
  for (int i=0; i<n5; i++ ){
    if (datfo_fc[i]>0.0) ep+=datfo_fc[i];
    else em+=datfo_fc[i];
    et+=datfo_fc[i];
    fomax=fmax(datfo_fc[i],fomax);
    fomin=fmin(datfo_fc[i],fomin);
  }
  printf("Fo-Fc %g %g %g %g %g\n",fomin,fomax,ep/n5*C[14],em/n5*C[14],et/n5*C[14] );
  fprintf(fof1,"F obs-fc map written by ShelXle\nDifferene density Fo-fc in fcf file transformed using fft\n");
  fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",na,
      (dx1.x+dy1.x+dz1.x)*0.5/a0,
      (dx1.y+dy1.y+dz1.y)*0.5/a0,
      (dx1.z+dy1.z+dz1.z)*0.5/a0
      //0.0,0.0,(dx.z+dy.z+dz.z)/a0
      );
  fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",n1,dx1.x/a0,dx1.y/a0,dx1.z/a0);
  fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",n2,dy1.x/a0,dy1.y/a0,dy1.z/a0);
  fprintf(fof1,"%5d%12.6f%12.6f%12.6f"  ,n3,dz1.x/a0,dz1.y/a0,dz1.z/a0);//no newline here because it is in atomlist
  fprintf(fof1,"%s",atomlist);
  for (int xi=0;xi<n1;xi++)
    for (int yi=0;yi<n2;yi++)
      for (int zi=0;zi<n3;zi++)
      {
        //fprintf(fof1,"%s%13.5E",(((i%6)==0)?"\n":""),datfo_fc[dex(xi,yi,zi)]*factor);
        fprintf(fof1,"%s",QString("%1%2").arg(((i%6)==0)?"\n":"").arg(datfo_fc[dex(xi,yi,zi)]*factor,13,'E',5).toStdString().c_str());
        i++;
      }
  fclose(fof1);
  }
  // */
//}

int FourXle::readHeader(const char *filename){
  /*! reads the header of an fcf file
  */
  FILE *f=NULL;
  char line[122],*dum=NULL;
  //size_t zlen=120;
  int ok=0;
  int i;double T,V;
  f=fopen(filename,"r");
  //printf("%d\n",__LINE__);
  if (f==NULL) return 3;
  //printf("%d\n",__LINE__);
  ns=0;
  sy[0][ns]=1.0;
  sy[1][ns]=0.0;
  sy[2][ns]=0.0;

  sy[3][ns]=0.0;
  sy[4][ns]=1.0;
  sy[5][ns]=0.0;

  sy[6][ns]=0.0;
  sy[7][ns]=0.0;
  sy[8][ns]=1.0;

  sy[9][ns]=0.0;
  sy[10][ns]=0.0;
  sy[11][ns]=0.0;
  ns=1;
  int listcode=0;
  do{
    dum=fgets(line,120,f);
    //printf("%d\n",__LINE__);
    if (dum==NULL){fclose(f);return 2;};
    if (feof(f)){fclose(f);return 2;};
    while (dum[0]==' ') dum++;
    if (!strncmp(dum,"_shelx_title",12)) {
      sscanf(line,"_shelx_title %[^\r\n]",titl);
      trimm(titl);
    }
    if (!strncmp(dum,"_exptl_crystal_F_000",20)){
      sscanf(line,"_exptl_crystal_F_000 %f",&f000);
      if (overridef000!=-1.0) f000=overridef000;
      //printf("F000 = %g %g %d %s\n",f000,overridef000,(overridef000!=-1.0f),line);
    }
    if (!strncmp(dum,"_shelx_refln_list_code",22)) {
      sscanf(line,"_shelx_refln_list_code %d",&listcode);
      //qDebug()<<listcode;
      if (listcode!=6) {fclose(f);return 1;}
    }
    if (!strncmp(dum,"_cell_length_a",14)) {
      sscanf(line,"_cell_length_a %lf",&C[0]);
    } 
    if (!strncmp(dum,"_cell_length_b",14)) {
      sscanf(line,"_cell_length_b %lf",&C[1]);
    } 
    if (!strncmp(dum,"_cell_length_c",14)) {
      sscanf(line,"_cell_length_c %lf",&C[2]);
    } 
    if (!strncmp(dum,"_cell_angle_alpha",17)) {
      sscanf(line,"_cell_angle_alpha %lf",&C[3]);
    } 
    if (!strncmp(dum,"_cell_angle_beta",16)) {
      sscanf(line,"_cell_angle_beta %lf",&C[4]);
    } 
    if (!strncmp(dum,"_cell_angle_gamma",17)) {
      sscanf(line,"_cell_angle_gamma %lf",&C[5]);
      for (i=0;i<3;i++){
        if (C[i]<0.1) return 2;
        T=.0174533*C[i+3];
        if (T<0.001) return 2;
        D[i]=sin(T);
        D[i+3]=cos(T);
        C[i+6]=(D[i]/(C[i]*C[i]));
      } 
      V=1.-D[3]*D[3]-D[4]*D[4]-D[5]*D[5]+2.*D[3]*D[4]*D[5];
      C[6]/=V;
      C[7]/=V;
      C[8]/=V;
      C[9]= 2.*sqrt(C[7]*C[8])*(D[4]*D[5]-D[3])/(D[2]*D[2]);
      C[10]=2.*sqrt(C[6]*C[8])*(D[3]*D[5]-D[4])/(D[0]*D[2]);
      C[11]=2.*sqrt(C[6]*C[7])*(D[3]*D[4]-D[5])/(D[0]*D[1]);
      C[14]=C[1]*C[2]*C[0]*sqrt(V);
      D[6]=C[1]*C[2]*D[0]/C[14];
      D[7]=C[0]*C[2]*D[1]/C[14];
      D[8]=C[0]*C[1]*D[2]/C[14];
    }
    if ((!strncmp(dum,"_symmetry_equiv_pos_as_xyz",26))||(!strncmp(dum,"_space_group_symop_operation_xyz",32))){
      //      char s1[50],s2[50],s3[50];
      //      char *kill=NULL,*nom=NULL,*div=NULL ;
      dum=fgets(line,120,f);
      trimm(line);
      while (strchr(line,'Y')) {
        QString sc=QString(line).toUpper().remove("SYMM").trimmed();
        sc=sc.remove("'");
        sc=sc.remove(" ");
        QStringList axe=sc.split(",");
        QStringList bruch;
        if (axe.size()!=3) return false;
        double _sx[3],_sy[3],_sz[3],t[3];
        for (int i=0; i<3; i++){
          _sx[i]=0;_sy[i]=0;_sz[i]=0;t[i]=0;
          if (axe.at(i).contains("-X")) {_sx[i]=-1.0;axe[i].remove("-X");}
          else if (axe.at(i).contains("X")) {_sx[i]=1.0;axe[i].remove("X");}
          if (axe.at(i).contains("-Y")) {_sy[i]=-1.0;axe[i].remove("-Y");}
          else if (axe.at(i).contains("Y")) {_sy[i]=1.0;axe[i].remove("Y");}
          if (axe.at(i).contains("-Z")) {_sz[i]=-1.0;axe[i].remove("-Z");}
          else if (axe.at(i).contains("Z")) {_sz[i]=1.0;axe[i].remove("Z");}
          if (axe.at(i).endsWith("+")) axe[i].remove("+");
          if (axe.at(i).contains("/")) {
            bruch=axe.at(i).split("/");
            if (bruch.size()==2) t[i]=bruch.at(0).toDouble() / bruch.at(1).toDouble();
          }
          else if (!axe.at(i).isEmpty()) t[i]=axe.at(i).toDouble();
        }
        sy[0][ns]=_sx[0];
        sy[1][ns]=_sy[0];
        sy[2][ns]=_sz[0];

        sy[3][ns]=_sx[1];
        sy[4][ns]=_sy[1];
        sy[5][ns]=_sz[1];

        sy[6][ns]=_sx[2];
        sy[7][ns]=_sy[2];
        sy[8][ns]=_sz[2];


        sy[9][ns]=t[0];
        sy[10][ns]=t[1];
        sy[11][ns]=t[2];


        strcpy(line,"");
        dum=fgets(line,120,f);
        trimm(line);
        ns++;

      }
    }
    if (!strncmp(dum,"_refln_phase_calc",17)) ok=1;
  }while((!ok)&&(!feof(f)));

  if (listcode!=6) return 1;
  for (int i=0; i<ns;i++){
    for (int n=i+1; n<ns;n++){
      int u=0,v=0;
      for (int j=0; j<9; j++){
        u+=abs(sy[j][n]-sy[j][i]);
        v+=abs(sy[j][n]+sy[j][i]);
      }
      if (fmin(u,v)>0.01) continue;
      for (int j=0; j<12; j++){
        sy[j][n]=sy[j][ns-1];
      }
      ns--;
    }
  }
  fclose(f);

  return 0;
}

void FourXle::sorthkl(int nr, Rec r[]){
  /*! sorts the reflection list
  */
  Rec *hilf= (Rec*) malloc(sizeof(Rec)*nr);
  if (hilf==NULL)return ; 
  int i,j,k,nj,ni,spalte;int index[4096];
  for (spalte=0; spalte<3; spalte++){
    j=-999999;
    k=999999;
    switch (spalte) {
      case 0: for (i=0; i<nr; i++){ j=(j<r[i].ih)?r[i].ih:j; k=(k>r[i].ih)?r[i].ih:k; } break;
      case 1: for (i=0; i<nr; i++){ j=(j<r[i].ik)?r[i].ik:j; k=(k>r[i].ik)?r[i].ik:k; } break;
      case 2: for (i=0; i<nr; i++){ j=(j<r[i].il)?r[i].il:j; k=(k>r[i].il)?r[i].il:k; } break;
    }
    nj=-k;
    ni=(nj+j+1);
    for (i=0; i<=ni; i++) index[i]=0;
    for (i=0; i<nr; i++){
      switch (spalte){
        case 0: j=r[i].ih+nj; break;
        case 1: j=r[i].ik+nj; break;
        case 2: j=r[i].il+nj; break;
      }
      index[j]++;/*brauch ich das? -->JA!*/
      hilf[i].ih=r[i].ih;
      hilf[i].ik=r[i].ik;
      hilf[i].il=r[i].il;
      hilf[i].fo=r[i].fo;
      hilf[i].so=r[i].so;
      hilf[i].fc=r[i].fc;
      hilf[i].phi=r[i].phi;
    }/*/4*/
    j=0;
    for (i=0; i<ni; i++){
      k=j;
      j+=index[i];
      index[i]=k;     
    }/*/5*/
    for (i=0; i<nr;i++){
      switch (spalte) {
        case 0: j=hilf[i].ih +nj;break;
        case 1: j=hilf[i].ik +nj;break;
        case 2: j=hilf[i].il +nj;break;
      }     
      index[j]++;   
      j=index[j]-1;
      r[j].ih=hilf[i].ih;
      r[j].ik=hilf[i].ik;
      r[j].il=hilf[i].il;
      r[j].fo=hilf[i].fo;
      r[j].so=hilf[i].so;
      r[j].fc=hilf[i].fc;
      r[j].phi=hilf[i].phi;
    }/*/6*/
  }/*/spalten*/
  free(hilf);
}

void FourXle::bewegt(V3 nm){/*!moves the rotation center to nm*/
 // printf("moved\n");
  V3 v , alturs=urs;
  mole->kart2frac(nm,v);
  urs=V3(1,1,1)-1.0*v;
  mole->frac2kart(urs,urs);
  if ((chgl->objCnt==acnt)&&(maptrunc==2)) return; 
  if (urs==alturs) return;
  gen_surface(false);
  chgl->pause=false;
  chgl->updateGL();
}

void FourXle::inimap(){/*! reinitializes the display lists for screenshots*/
  gen_surface(false);

}

void FourXle::gen_surface(bool neu,int imin,int imax){
  GenMapWorker *w = new GenMapWorker(chgl, mole, this, neu, imin, imax);
  chgl->habzutun=true;
  //connect(this,SIGNAL(haltMW()),w,SLOT(terminate()));
  connect(w, SIGNAL(finished()),this, SLOT(reportInfo()));
  connect(w, SIGNAL(finished()), chgl, SLOT(updateGL()));
  connect(w, SIGNAL(finished()), w, SLOT(deleteLater()));
  connect(w, SIGNAL(finished()),chgl, SLOT(fertig()));
  w->start();
}
//GenMapWorker



void GenMapWorker::gen_surface(bool neu,int imin,int imax){
  if (fxle->noGoMap==NULL) return;
  /*! creates iso surfaces for fo-fc- fo-fc+ fo+ and fo- maps if neu then the iso values are calculated fro the sigma value of the map.
   *
   */
  if (!chgl->neutrons) imax=qMin(3,imax);
  V3 v=V3(0,0,0);
  mole->kart2frac(chgl->altemitte,v);
  fxle->urs=V3(1,1,1)-1.0*v;
  mole->frac2kart(fxle->urs,fxle->urs);
  //chgl->pause=true;
  chgl->fVertexes.resize(28);
  chgl->fNormals.resize(28);
  fxle->suti.start();
  fxle->disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
  fxle->disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
  fxle->disconnect(chgl,SIGNAL(inimibas()),0,0);
  //  if ((neu)||(mole->showatoms.size()!=oldatomsize)||(chgl->unhide->isVisible())){
  for (int no=0;no<(fxle->n5*27);no++) fxle->noGoMap[no]=0;
  mole->kart2frac(fxle->dx,fxle->dxc);
  mole->kart2frac(fxle->dy,fxle->dyc);
  mole->kart2frac(fxle->dz,fxle->dzc);
  int incx = (int) (fxle->cutrange->value()/fxle->dx.x);
  int incy = (int) (fxle->cutrange->value()/fxle->dy.y);
  int incz = (int) (fxle->cutrange->value()/fxle->dz.z);
  int incmin=qMin(incx,qMin(incy,incz));
  incmin*=incmin;
  for (int g=0; g<mole->showatoms.size();g++){
    if (mole->showatoms.at(g).hidden) continue;
    mole->kart2frac(mole->showatoms.at(g).pos,fxle->oc);
    int ax=(int)((fxle->oc.x)/fxle->dxc.x-0.499), bx=(int)((fxle->oc.y)/fxle->dyc.y-0.499), cx=(int)((fxle->oc.z)/fxle->dzc.z-0.499);
    for (int aa=ax-incx; aa<ax+incx; aa++){
      for (int bb=bx-incy; bb<bx+incy; bb++){
        for (int cc=cx-incz; cc<cx+incz; cc++){
          if (incmin<((aa-ax)*(aa-ax)+(bb-bx)*(bb-bx)+(cc-cx)*(cc-cx))) continue;
          int dEx=fxle->dex3(aa,bb,cc);
          if ((dEx>0)&&(dEx<(27*fxle->n5-1)))
            fxle->noGoMap[dEx]=1;
        }
      }
    }
  }
  //  oldatomsize=mole->showatoms.size();
  //  }
  fxle->tri=0;
  float FXO=0.1f, FXD=0.01f;
  QString fxvaluesfile = fxle->resDir.absolutePath()+"_fourxle.fixv";
  //  qDebug()<<fxvaluesfile ;
  if (fxle->keepIso->isChecked()){
    if (QFile::exists(fxvaluesfile)){
      QFile fxfi(fxvaluesfile);
      fxfi.open(QIODevice::ReadOnly|QIODevice::Text);
      QStringList lll = QString(fxfi.readAll()).split(' ');
      if (lll.size()==2){
        FXO=lll.at(0).toFloat();
        FXD=lll.at(1).toFloat();
      }
      fxfi.close();
    }else{
      QFile fxfi(fxvaluesfile);
      fxfi.open(QIODevice::WriteOnly|QIODevice::Text);
      FXO = fxle->sigma[1];
      FXD = fxle->sigma[0];
      fxfi.write(QString("%1 %2").arg(FXO,0,'g').arg(FXD,0,'g').toLocal8Bit());
      fxfi.close();
    } 
  }else{
    QFile::remove(fxvaluesfile);
  }
  printf("sigma0 %f sigma1 %f \n",fxle->sigma[0]*2.7f,fxle->sigma[1]*1.2f);
  for (int fac=imin; fac<imax; fac++){
    switch (fac){
      case 0:                   
        if (neu) fxle->iso[1]=(fxle->keepIso->isChecked())?FXD*2.7f:fxle->sigma[0]*2.7f;
        else fxle->iso[1]=fabs(fxle->iso[1]);
        fxle->mtyp=1;
        break;
      case 1:
        if (neu) fxle->iso[1]=(fxle->keepIso->isChecked())?-FXD*2.7f:-fxle->sigma[0]*2.7f;
        else fxle->iso[1]=-fabs(fxle->iso[1]);
        fxle->mtyp=1;
        break;
      case 2:                   
        if (neu) fxle->iso[0]=(fxle->keepIso->isChecked())?FXO*1.2f:fxle->sigma[1]*1.2f;
        //   printf("blau %g %d\n",fxle->iso[0],neu);
        fxle->mtyp=0;
        break;
      case 3:                   
        if (neu) fxle->iso[0]=(fxle->keepIso->isChecked())?-FXO*1.2f:-fxle->sigma[1]*1.2f;
        else fxle->iso[0]=-fabs(fxle->iso[0]);
        //   printf("orange %g %d\n",fxle->iso[0],neu);
        fxle->mtyp=0;
        break;
    }
    {
      int ix,iy,iz;
      for (iz=0; iz<fxle->n3;iz++){
        for (iy=0; iy<fxle->n2;iy++){
          for (ix=0; ix<fxle->n1;ix++){
            fxle->CalcVertex(ix,iy,iz);
          }
        }
      }
    }

    int mper=(chgl->niceTrans->isChecked())?6:1;
    for (int pers=0;pers<mper;pers++){        

        if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
        if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
        if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
      int facnpers=fac+pers*4;
      chgl->fVertexes[facnpers].clear();
      chgl->fNormals[facnpers].clear();

      int h,k,l;
      int step=0;
      switch (pers) {
        case 0:
          //fprintf(stderr, "makeElement case 0 n1=%d n2=%d n3=%d N4=%d fac+pers=%d nodes %p %p %p\n",fxle->n1,fxle->n2,fxle->n3,fxle->n4,facnpers,fxle->nodex,fxle->nodey,fxle->nodez);
          for ( h=0; h<fxle->n1;h++){
            for ( k=0; k<fxle->n2;k++){
              for ( l=0; l<fxle->n3;l++){
                if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          //printf("makeElement made %d Vertexes %p %p %p\n", chgl->fVertexes[facnpers].size(), fxle->nodex, fxle->nodey, fxle->nodez);
          break;
        case 1:
          //fprintf(stderr, "makeElement case 1\n");
          for ( h=fxle->n1; h>=0;h--){
            for ( k=fxle->n2; k>=0;k--){
              for ( l=fxle->n3; l>=0;l--){
                if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          break;
        case 2:
          //fprintf(stderr, "makeElement case 2\n");
          for ( k=0; k<fxle->n2;k++){
            for ( h=0; h<fxle->n1;h++){
              for ( l=0; l<fxle->n3;l++){
                  if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                  if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                  if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          break;
        case 3:
          //fprintf(stderr, "makeElement case 3\n");
          for ( k=fxle->n2; k>=0;k--){
            for ( h=fxle->n1; h>=0;h--){
              for ( l=fxle->n3; l>=0;l--){
                  if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                  if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                  if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          break;
        case 4:
          //fprintf(stderr, "makeElement case 4\n");
          for ( l=0; l<fxle->n3;l++){
            for ( h=0; h<fxle->n1;h++){
              for ( k=0; k<fxle->n2;k++){
                  if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                  if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                  if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          break;
        case 5:
          //fprintf(stderr, "makeElement case 5\n");
          for ( l=fxle->n3; l>=0;l--){
            for ( h=fxle->n1; h>=0;h--){
              for ( k=fxle->n2; k>=0;k--){
                  if (fxle->nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
                  if (fxle->nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
                  if (fxle->nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
                fxle->MakeElement(h,k,l,facnpers);
              }
            }
            step++;step%=mper;
          }
          break;
      }
    }
  }
  fxle->iso[0]=fabs(fxle->iso[0]);

  fxle->info=QString("<b>%13 Map:</b><font color=\"%12\">%1 e&Aring;<sup>-3</sup></font>"
      "<font color=\"%11\"> %2 e&Aring;<sup>-3</sup> (&sigma; = %7) </font><br><font color=grey> Hint:  [%3 Scroll (up or down)] to change. </font><br>"
      "<b>Fo-Map:</b><font color=\"%10\">%4  e&Aring;<sup>-3</sup> (&sigma; = %8)</font><br><font color=grey> Hint:  [%5 Scroll (up or down)] to change. </font><br>Time for creating map surfaces <b>%6 s</b>. %9 Triangles drawn.<br>")
    .arg(fxle->iso[1],6,'g',2)
    .arg(-fxle->iso[1],6,'g',2)
    .arg(QKeySequence(Qt::ControlModifier).toString(QKeySequence::NativeText))
    .arg(fxle->iso[0],6,'g',2)
    .arg(QKeySequence(Qt::ShiftModifier).toString(QKeySequence::NativeText))
    .arg(fxle->suti.elapsed()/1000.0,1,'f',1)
    .arg(fxle->sigma[0],6,'g',2)
    .arg(fxle->sigma[1],6,'g',2)
    .arg(fxle->tri)
    .arg(chgl->fopc.name())
    .arg(chgl->dipc.name())
    .arg(chgl->dimc.name())
    .arg((chgl->dipc.name() == QString("#008080"))?"BODD":"Fo-Fc1")//quick and dirty check for BODD
    ;  
  //chgl->pause=false;
  fxle->connect(chgl,SIGNAL(inimibas()),fxle,SLOT(inimap()));
  fxle->connect(chgl,SIGNAL(neuemitte(V3)),fxle, SLOT(bewegt(V3)));
  fxle->connect(chgl,SIGNAL(diffscroll(int ,int )),fxle,SLOT(change_iso(int ,int )));
  chgl->fofcact->setVisible(!chgl->fVertexes[0].isEmpty()||!chgl->fVertexes[1].isEmpty());
  chgl->foact->setVisible(!chgl->fVertexes[2].isEmpty());
  //printf("wer ist empty %d %d %d %d\n",!chgl->fVertexes[0].isEmpty(),!chgl->fVertexes[1].isEmpty(),!chgl->fVertexes[2].isEmpty(),!chgl->fVertexes[0].isEmpty()||!chgl->fVertexes[1].isEmpty());
  //  printf("Fomap ist da? %d Fo-Fcmap ist da? %d \n",chgl->fVertexes[2].isEmpty(),chgl->fVertexes[0].isEmpty());

  fxle->acnt=chgl->objCnt;
}


void FourXle::reportInfo(){
  emit bigmessage(info);
}

void FourXle::change_iso(int numsteps,int diff){
  /*! canges the iso value of the fo or fo-fc maps and redraws them
  */
  //printf("isochaged\n");
  iso[diff]=fabs(iso[diff]);
  iso[diff]+=iso[diff]*numsteps/10.0f;
  int mi=0,ma=5;
  switch (diff) {
    case 0: 
      mi=2;ma=4;
      break;
    case 1:
      mi=0; ma=2;
      break;

  }    
  gen_surface(false,mi,ma);
  chgl->updateGL();
}

void FourXle::foactDestroyed() {
  if (chgl!=NULL) chgl->foact = NULL;
}

void FourXle::fofcactDestroyed() {
  if (chgl!=NULL)chgl->fofcact = NULL;
}

void FourXle::chglDestroyed() {
  //fprintf(stderr,"chgl is dead omg ...argh...\n");

  chgl = NULL;
}

void FourXle::CalcVertex( int ix, int iy, int iz) {
  V3 mdz=(0.5*dx)+(0.5*dy)+(0.5*dz);
  V3 o,fl,m2u=V3(0.5,0.5,0.5);
  mole->frac2kart(m2u,m2u);
  double vo=0, vx=0,vy=0,vz=0;
  int idx=dex(ix,iy,iz);
  nodex[idx].flag=0;
  nodey[idx].flag=0;
  nodez[idx].flag=0;
  if (mtyp==0){//*datfo,*datfo_fc,*datf1_f2
    vo = datfo[ idx]   -iso[mtyp];
    vx = datfo[ dex(ix+1,iy,iz)]   -iso[mtyp];
    vy = datfo[ dex(ix,iy+1,iz)]   -iso[mtyp];
    vz = datfo[ dex(ix,iy,iz+1)]   -iso[mtyp];
  }else {
    vo = datfo_fc[idx]   -iso[mtyp];
    vx = datfo_fc[dex(ix+1,iy,iz)]   -iso[mtyp];
    vy = datfo_fc[dex(ix,iy+1,iz)]   -iso[mtyp];
    vz = datfo_fc[dex(ix,iy,iz+1)]   -iso[mtyp];
  }
  V3 nor=V3(0,0,0);//Normalize((vx-vo)*dx+(vy-vo)*dy+(vz-vo)*dz);
  if (Intersect(vo,vx)) {
    if (maptrunc==2) o=dx*((vo/(vo-vx))+ix) + dy*iy + dz*iz+m2u;
    else  o=dx*((vo/(vo-vx))+ix) + dy*iy + dz*iz+urs;

    mole->kart2frac(o,o);
    o+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.x),floor(o.y),floor(o.z));
    o+=-1.0*fl;
    o+=V3(0.5,0.5,0.5);
    mole->frac2kart(o,o);
    if (maptrunc!=2) o+=-1.0*urs;
    else o+=-1.0*m2u;
    o+=mdz;
    //    orte.append(o);
    nodex[idx].vertex=o;
    nodex[idx].normal=nor;
    nodex[idx].flag=1;
  }
  if (Intersect(vo,vy)) {
    if (maptrunc==2)o=dx*ix + dy*((vo/(vo-vy))+iy) + dz*iz+m2u;
    else o=dx*ix + dy*((vo/(vo-vy))+iy) + dz*iz+urs;
    mole->kart2frac(o,o);
    o+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.x),floor(o.y),floor(o.z));
    o+=-1.0*fl;
    o+=V3(0.5,0.5,0.5);
    mole->frac2kart(o,o);
    if (maptrunc!=2) o+=-1.0*urs;
    else o+=-1.0*m2u;
    o+=mdz;
    //    orte.append(o);
    nodey[idx].vertex=o;
    nodey[idx].normal=nor;
    nodey[idx].flag=1;
  }
  if (Intersect(vo,vz)) {
    if (maptrunc==2)o=dx*ix + dy*iy + dz*((vo/(vo-vz))+iz)+m2u;
    else o=dx*ix + dy*iy + dz*((vo/(vo-vz))+iz)+urs;
    mole->kart2frac(o,o);
    o+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.x),floor(o.y),floor(o.z));
    o+=-1.0*fl;
    o+=V3(0.5,0.5,0.5);
    mole->frac2kart(o,o);
    if (maptrunc!=2) o+=-1.0*urs;
    else o+=-1.0*m2u;
    o+=mdz;
    //    orte.append(o);
    nodez[idx].vertex=o;
    nodez[idx].normal=nor;
    nodez[idx].flag=1;
  }
}

void FourXle::makeFaces(int n, int faps, FNode poly[] ){
  if (n<3) return;  //less then 3 vertices  -> nothing to do
  V3 mid_ver = V3(0,0,0);
  V3 mid_nor2 = V3(0,0,0);
  for(int i=0; i<n; i++ ){
    if ((i>0)&&(i<n-1)) mid_nor2 += Normalize((poly[i].vertex-poly[0].vertex)%(poly[i+1].vertex-poly[0].vertex));
    mid_ver += poly[i].vertex;
  }
  mid_nor2=Normalize(mid_nor2);
  for(int i=0; i<n; i++ ){
    poly[i].normal = mid_nor2;
  }
  mid_ver *= (1.0/n);
  //mid_nor2= Normalize(mid_nor2);
  V3 mit=V3(1,1,1);
  mole->frac2kart(mit,mit);
  for (int w=0; w<27; w++){
    if (maptrunc==0){w=13;}
    else if ((maptrunc==1)&&(Distance(mit-urs,mid_ver+delDA[w])>(map_radius*map_radius))) continue;
    else if (maptrunc==2) {
      mole->kart2frac(mid_ver+delDA[w],oc);
      int  ax = static_cast<int>((oc.x)/dxc.x-0.499),
           bx = static_cast<int>((oc.y)/dyc.y-0.499),
           cx = static_cast<int>((oc.z)/dzc.z-0.499);
      int DEX=dex3(ax,bx,cx);
      if (!noGoMap[DEX]) continue;
    }    
    for (int k=1; k<n-1;k++){
      tri++;

      chgl->fVertexes[faps].append(poly[k].vertex.x+delDA[w].x);
      chgl->fVertexes[faps].append(poly[k].vertex.y+delDA[w].y);
      chgl->fVertexes[faps].append(poly[k].vertex.z+delDA[w].z);

      chgl->fVertexes[faps].append(poly[(k+1)].vertex.x+delDA[w].x);
      chgl->fVertexes[faps].append(poly[(k+1)].vertex.y+delDA[w].y);
      chgl->fVertexes[faps].append(poly[(k+1)].vertex.z+delDA[w].z);

      chgl->fVertexes[faps].append(poly[0].vertex.x+delDA[w].x);
      chgl->fVertexes[faps].append(poly[0].vertex.y+delDA[w].y);
      chgl->fVertexes[faps].append(poly[0].vertex.z+delDA[w].z);


      chgl->fNormals[faps].append(poly[k].normal.x);
      chgl->fNormals[faps].append(poly[k].normal.y);
      chgl->fNormals[faps].append(poly[k].normal.z);

      chgl->fNormals[faps].append(poly[(k+1)].normal.x);
      chgl->fNormals[faps].append(poly[(k+1)].normal.y);
      chgl->fNormals[faps].append(poly[(k+1)].normal.z);

      chgl->fNormals[faps].append(poly[0].normal.x);
      chgl->fNormals[faps].append(poly[0].normal.y);
      chgl->fNormals[faps].append(poly[0].normal.z);
    }
    if (maptrunc==0) w=27;
  }
}

int FourXle::IndexSelected( FNode& node0, FNode& node1, FNode& node2, FNode& node3 ) {
  if( node1 && node2 && node3 ){
    double d1 = Distance( node0.vertex, node1.vertex) +  Distance( node3.vertex, node2.vertex) ;
    double d2 = Distance( node0.vertex, node2.vertex) +  Distance( node3.vertex, node1.vertex );
    if( d1 > d2 ) return 2; else return 1;
  }else{
    if(      node1 )   return 1;
    else if( node2 )   return 2;
    else if( node3 )   return 3;
  }
  return 0;

}

void FourXle::MakeElement(int ix, int iy, int iz , int faps) {
    QMutexLocker locker(&mutex);
  int conn[12][2][4] = {
    {{ 0, 1, 7, 6}, { 0, 2, 8, 3}},  //  0
    {{ 1, 2, 5, 4}, { 1, 0, 6, 7}},  //  1
    {{ 2, 0, 3, 8}, { 2, 1, 4, 5}},  //  2
    {{ 3, 8, 2, 0}, { 3, 4,10, 9}},  //  3
    {{ 4, 3, 9,10}, { 4, 5, 2, 1}},  //  4
    {{ 5, 4, 1, 2}, { 5, 6, 9,11}},  //  5
    {{ 6, 5,11, 9}, { 6, 7, 1, 0}},  //  6
    {{ 7, 6, 0, 1}, { 7, 8,11,10}},  //  7
    {{ 8, 7,10,11}, { 8, 3, 0, 2}},  //  8
    {{ 9,10, 4, 3}, { 9,11, 5, 6}},  //  9
    {{10,11, 8, 7}, {10, 9, 3, 4}},  // 10
    {{11, 9, 6, 5}, {11,10, 7, 8}}   // 11
  };
  FNode node[12];
  FNode polygon[12];
  /*
  if (nodex==NULL) {fprintf(stderr,"NODE X is NULL!\n");return;}
  if (nodey==NULL) {fprintf(stderr,"NODE Y is NULL!\n");return;}
  if (nodez==NULL) {fprintf(stderr,"NODE Z is NULL!\n");return;}
  if (
         ((ix+iy*s1+iz*n4)                  %n5)<0||
         ((ix+iy*s1+iz*n4)                  %n5)<0||
         ((ix+iy*s1+iz*n4)                  %n5)<0||
         ((ix+iy*s1+((iz+1)%n3)*n4)         %n5)<0||
         ((ix+iy*s1+((iz+1)%n3)*n4)         %n5)<0||
         ((ix+((iy+1)%n2)*s1+iz*n4)         %n5)<0||
         ((ix+((iy+1)%n2)*s1+iz*n4)         %n5)<0||
         ((((1+ix)%n1)+iy*s1+iz*n4)         %n5)<0||
         ((((1+ix)%n1)+iy*s1+iz*n4)         %n5)<0||
         ((ix+((iy+1)%n2)*s1+((iz+1)%n3)*n4)%n5)<0||
         ((((ix+1)%n1)+iy*s1+((iz+1)%n3)*n4)%n5)<0||
         ((((ix+1)%n1)+((iy+1)%n2)*s1+iz*n4)%n5)<0
         )printf("EEEEK!\n");
         */
  node[ 0] = nodex[(ix+iy*n1+iz*n4)                  %n5];        // 000x
  node[ 1] = nodey[(ix+iy*n1+iz*n4)                  %n5];        // 000y
  node[ 2] = nodez[(ix+iy*n1+iz*n4)                  %n5];        // 000z
  node[ 3] = nodex[(ix+iy*n1+((iz+1)%n3)*n4)         %n5];    // 001y
  node[ 4] = nodey[(ix+iy*n1+((iz+1)%n3)*n4)         %n5];    // 001z
  node[ 5] = nodez[(ix+((iy+1)%n2)*n1+iz*n4)         %n5];    // 010x
  node[ 6] = nodex[(ix+((iy+1)%n2)*n1+iz*n4)         %n5];    // 010y
  node[ 7] = nodey[(((1+ix)%n1)+iy*n1+iz*n4)         %n5];      // 100y
  node[ 8] = nodez[(((1+ix)%n1)+iy*n1+iz*n4)         %n5];      // 100z
  node[ 9] = nodex[(ix+((iy+1)%n2)*n1+((iz+1)%n3)*n4)%n5];// 011x
  node[10] = nodey[(((ix+1)%n1)+iy*n1+((iz+1)%n3)*n4)%n5];// 101y
  node[11] = nodez[(((ix+1)%n1)+((iy+1)%n2)*n1+iz*n4)%n5];// 110z
  if ((static_cast<char>(node[0])+node[1]+node[2]+node[3]+node[4]+node[5]
        +node[6]+node[7]+node[8]+node[9]+node[10]+node[11])==0) return;
  for( int is=0; is<12; is++ ) {
    if( !node[is] ) continue;

    int n=0, i=is, m=0;//,ai=i;
    double dis;
    dis=0.0;
    do {
      polygon[n++]= node[i];
      int sol = IndexSelected(
          node[conn[i][m][0]],
          node[conn[i][m][1]],
          node[conn[i][m][2]],
          node[conn[i][m][3]]);
      //ai=i;
      i = conn[i][m][sol];
      if( sol == 2 ) m ^= 1;
      dis+=Distance(polygon[0].vertex,node[i].vertex);
      node[i].flag= 0;      
    }while( (i!=is)&&(n<11) );
    if (n>=3) {
      if (dis<5){
        //for (int polni=0; polni<n; polni++){polygon[polni].normal+=Normalize((polygon[(polni+1)%n].vertex-polygon[polni].vertex)%(polygon[(polni+n-1)%n].vertex-polygon[polni].vertex));}
        makeFaces( n, faps, polygon );
      }//*
      else {
        int axe=0;
        double delx=0,dely=0,delz=0;
        double mind=100000000;
        V3 minp=V3(10000,10000,10000),lihiun=V3(-10000,-10000,-10000);
        int minii=0;
        for (int polni=1; polni<=n; polni++){
          delx+=fabs(polygon[polni-1].vertex.x - polygon[polni%n].vertex.x);
          dely+=fabs(polygon[polni-1].vertex.y - polygon[polni%n].vertex.y);
          delz+=fabs(polygon[polni-1].vertex.z - polygon[polni%n].vertex.z);
          if (Distance(polygon[polni%n].vertex,lihiun)<mind) {
            mind=Distance(polygon[polni%n].vertex,minp);
            minii=polni%n;
          }
        }
        minp=polygon[minii].vertex;
        axe|=(delx>1)?1:0;
        axe|=(dely>1)?2:0;
        axe|=(delz>1)?4:0;
        for (int polni=0; polni<n; polni++){
          V3 neo=polygon[polni].vertex;
          double lang =Distance(minp,neo);
          if ((lang>Distance(minp,neo+dx*n1))&&(axe&1)) neo+=dx*n1;
          else if ((lang>Distance(minp,neo-dx*n1))&&(axe&1)) neo+=-n1*dx;
          lang =Distance(minp,neo);
          if ((lang>Distance(minp,neo+dy*n2))&&(axe&2)) neo+=n2*dy;
          else if ((lang>Distance(minp,neo-dy*n2))&&(axe&2)) neo+=-n2*dy;
          lang =Distance(minp,neo);
          if ((lang>Distance(minp,neo+n3*dz))&&(axe&4)) neo+=n3*dz;
          else if ((lang>Distance(minp,neo-n3*dz))&&(axe&4)) neo+=-n3*dz;
          polygon[polni].vertex=neo;
        }
        dis=0;
        for (int polni=1; polni<=n; polni++){
          dis+=Distance(polygon[polni-1].vertex , polygon[polni%n].vertex);
        }
        if (dis<5)
          makeFaces( n, faps, polygon );
      }
      // */
    }
  }

}

float FourXle::sigmaFo(){
  float DM=0.0f;
  float DS=0.0f;
  float f=0.0f;
  for (int i=0; i<n5; i++){
    f=datfo[i];
    DM+=f;
    DS+=f*f;
  }
  return sqrt((DS/n5)-((DM/n5)*(DM/n5)));
}

float FourXle::sigmaFoFc(){
  float DM=0.0f;
  float DS=0.0f;
  float f=0.0f;
  for (int i=0; i<n5; i++){
    f=datfo_fc[i];
    DM+=f;
    DS+=f*f;
  }
  return sqrt((DS/n5)-((DM/n5)*(DM/n5)));
}

void FourXle::jnk(){
  if(!n5)return;
  float mini=9e37f,maxi=-9e37f;
  float f,r,fstep;
  double DM,DS,sigma,w,e_net,e_gross;
  r=powf((3*(n1-1)*(n2-1)*(n3-1)),1.0f/3.0f);
  //  printf("%d %d %d %d %f\n",n1,n2,n3,n5,r);
  float *datfo_fcstp=(float*) malloc(sizeof(float)*n5);
  if (datfo_fcstp==NULL)return ; 
  float invstep=100.0f;
  float step=1.0f/invstep;
  float rhomind2=99999.0f,rhomaxd2=-99999.0f;
  QMap<float,int> hash;
  QMap<float,float> hashf;
  DM=0.0;
  DS=0.0;
  w=0.0;
  for (int i=0; i<n5; i++){
    f=datfo_fc[i];
    mini=qMin(f,mini);
    maxi=qMax(f,maxi);
    DM+=f;
    DS+=f*f;
    w+=fabs(f);
    datfo_fcstp[i]=floorf(f*invstep);
  }// */
  e_net=(DM/n5)*C[14];//C[14] ist das Volumen der UC
  e_gross=w/(2*n5)*C[14];
  sigma=sqrt((DS/n5)-((DM/n5)*(DM/n5)));
  //printf("%g %g %g",e_net,e_gross,sigma);
  for (int zi=0;zi<n3;zi++)
    for (int yi=0;yi<n2;yi++)
      for (int xi=0;xi<(n1-1);xi++){
        if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi+1,yi,zi)]){
          fstep=datfo_fcstp[dex(xi,yi,zi)];
          int ze=fstep-datfo_fcstp[dex(xi+1,yi,zi)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi+1,yi,zi)]){
          fstep=datfo_fcstp[dex(xi+1,yi,zi)];
          int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }
      }
  for (int xi=0;xi<n1;xi++)
    for (int zi=0;zi<n3;zi++)
      for (int yi=0;yi<(n2-1);yi++){
        if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi,yi+1,zi)]){
          fstep=datfo_fcstp[dex(xi,yi,zi)];
          int ze=fstep-datfo_fcstp[dex(xi,yi+1,zi)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi,yi+1,zi)]){
          fstep=datfo_fcstp[dex(xi,yi+1,zi)];
          int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }
      }

  for (int yi=0;yi<n2;yi++)
    for (int xi=0;xi<n1;xi++)
      for (int zi=0;zi<(n3-1);zi++){
        if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi,yi,zi+1)]){
          fstep=datfo_fcstp[dex(xi,yi,zi)];
          int ze=fstep-datfo_fcstp[dex(xi,yi,zi+1)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi,yi,zi+1)]){
          fstep=datfo_fcstp[dex(xi,yi,zi+1)];
          int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
          for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
        }
      }
  free(datfo_fcstp);
  //float range=fmaxf(fabsf(mini),fabsf(maxi));
  QDialog *jnkdlg = new QDialog();
  QGraphicsScene*scene= new QGraphicsScene(-30,-50,550,586);
  scene->setBackgroundBrush(QBrush(QColor("#e9f7d6")));
  scene->clear ();
  QGraphicsItem *itm;
  for (int i=0; i<21;i++){
    itm=scene->addLine(i*25,0,i*25,500,(i%5)?QPen(QColor("#cbdbbb"),0):QPen(QColor("#959d9d"),0));
    itm->setData(0,-1);
  }
  for (int i=0; i<31;i++){
    itm=scene->addLine(0,i*16.66666666666667,500,i*16.66666666666667,(i%5)?QPen(QColor("#cbdbbb"),0):QPen(QColor("#959d9d"),0));
    itm->setData(0,-1);
  }
  itm=scene->addLine(250,0,250,500,QPen(QColor("#000000"),0));
  itm=scene->addLine(0,83.33333333333333,500,83.33333333333333,QPen(QColor("#000000"),0));
  itm=scene->addLine(0,250,500,250,QPen(QColor("#000000"),0));
  itm=scene->addLine(0,416.6666666666667,500,416.6666666666667,QPen(QColor("#000000"),0));
  QGraphicsTextItem *txt = scene->addText("3");
  txt->moveBy(-15,73.33333333333333);
  txt = scene->addText("2");
  txt->moveBy(-15,240);
  txt = scene->addText("1");
  txt->moveBy(-15,406.6666666666667);
  txt = scene->addText("-1.0");
  txt->moveBy(-15,500);
  txt = scene->addText("0");
  txt->moveBy(241,500);
  txt = scene->addText("1.0");
  txt->moveBy(485,500);


  QMapIterator<float, int> i(hash);
  while (i.hasNext()) {
    i.next();
    f=hashf[i.key()]=logf(i.value())/logf(r);
    scene->addEllipse(250+(i.key())*250,500-((f-0.5)/3.)*500,4,4,QPen(Qt::NoPen),QBrush(QColor("#0907e6")));

    if (f>2){
      rhomind2=qMin(rhomind2,i.key());
      rhomaxd2=qMax(rhomaxd2,i.key());
    }
    // printf("jnk: %g %g %g %g %g %g\n",i.key(),f,rhomind2,rhomaxd2,250+(i.key())*250,500-((f-0.5)/3.)*500);
  }
  // printf("\ndf(0)= %g %d %g %g\n",hashf.value(0.0f),hash.value(0.0f),rhomind2,rhomaxd2);

  float m,b;
  m=(hashf.value(rhomind2)-hashf.value(rhomind2-step))/step;
  b=hashf.value(rhomind2)-m*rhomind2;
  rhomind2=(2.0f-b)/m;
  m=(hashf.value(rhomaxd2)-hashf.value(rhomaxd2-step))/step;
  b=hashf.value(rhomaxd2)-m*(rhomaxd2);
  rhomaxd2=(2.0f-b)/m;
  // printf("df(0)= %g %d %g %g\n",hashf.value(0.0f),hash.value(0.0f),rhomind2,rhomaxd2);
  txt = scene->addText(QString("df(0) = %1").arg((double)hashf.value(0.0f),7,'f',2),QFont("Courier",10));
  txt->moveBy(328,2);
  txt = scene->addText(QString("min(d=2) = %1 eA^-3").arg((double)rhomind2,7,'f',3),QFont("Courier",10));
  txt->moveBy(305,18);
  txt = scene->addText(QString("max(d=2) = %1 eA^-3").arg((double)rhomaxd2,7,'f',3),QFont("Courier",10));
  txt->moveBy(305,34);
  txt = scene->addText(QString("min      = %1 eA^-3").arg((double)mini,7,'f',3),QFont("Courier",10));
  txt->moveBy(305,50);
  txt = scene->addText(QString("max      = %1 eA^-3").arg((double)maxi,7,'f',3),QFont("Courier",10));
  txt->moveBy(305,66);
  txt = scene->addText(QString("e_net    = %1 e").arg((double)e_net,8,'g',3),QFont("Courier",10));
  txt->moveBy(295,82);
  txt = scene->addText(QString("e_gross  = %1 e").arg((double)e_gross,8,'g',3),QFont("Courier",10));
  txt->moveBy(295,98);
  txt = scene->addText(QString("         = %1 eA^-3").arg((double)sigma,8,'f',3),QFont("Courier",10));
  txt->moveBy(295,116);
  txt = scene->addText(QString("nx ny nz  %1 x %2 x %3 = %4").arg(n1).arg(n2).arg(n3).arg(n5),QFont("Courier",9));
  txt->moveBy(30,2);
  txt = scene->addText(QString("fractal dimension vs. residual density"),QFont("Helvetica",14,QFont::Bold));
  txt->moveBy(100,-40);
  txt = scene->addText(QString("df"),QFont("Helvetica",12,QFont::Bold));
  txt->moveBy(-20,20);
  txt = scene->addText(QString("0"),QFont("Helvetica",12,QFont::Bold));
  txt->moveBy(460,500);

#ifdef __APPLE__
  txt = scene->addText(QString::fromUtf8("𝛒"),QFont("Helvetica",16,QFont::Bold));
  txt->moveBy(450,490);
  txt = scene->addText(QString::fromUtf8("𝛒"),QFont("Helvetica",14));
  txt->moveBy(295,8);
  txt = scene->addText(QString::fromUtf8("𝛒"),QFont("Helvetica",14));
  txt->moveBy(295,24);
  txt = scene->addText(QString::fromUtf8("𝛒"),QFont("Helvetica",14));
  txt->moveBy(295,40);
  txt = scene->addText(QString::fromUtf8("𝛒"),QFont("Helvetica",14));
  txt->moveBy(295,56);
  txt = scene->addText(QString::fromUtf8("𝛔"),QFont("Helvetica",14));
  txt->moveBy(295,106);
#else
  txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",16,QFont::Bold));
  txt->moveBy(450,490);
  txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
  txt->moveBy(295,8);
  txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
  txt->moveBy(295,24);
  txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
  txt->moveBy(295,40);
  txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
  txt->moveBy(295,56);
  txt = scene->addText(QString::fromUtf8("σ"),QFont("Helvetica",14));
  txt->moveBy(295,106);
#endif  // */
  txt = scene->addText(QString("Please cite as: 'K. Meindl, J. Henn, Acta Cryst., 2008, A64, 404-418.'"),QFont("Helvetica",9));
  txt->moveBy(100,516);
  QGraphicsView *view = new QGraphicsView(scene,jnkdlg);

  QVBoxLayout *lt = new QVBoxLayout(jnkdlg);
  lt->addWidget(view);
  jnkdlg->show();


}


inline int FourXle::dex(int x,int y, int z){
  /*! dex is used to adress elemennts of a one dimensional array by three indizes like it is a 3 dimensional array
   * @param x,y,z tree dimensional indices 
   * if x is < 0 or > n1 it is not a problem because % is used to clamp it.
   * if y is < 0 or > n2 it is not a problem because % is used to clamp it.
   * if z is < 0 or > n3 it is not a problem because % is used to clamp it.
   * \returns index of an 1 dimensional array 
   */
  x=(x+n1)%n1;
  y=(y+n2)%n2;
  z=(z+n3)%n3;
  return x+n1*(y+n2*z);
}

inline int FourXle::dex3(int x,int y, int z){
  int n31=3*n1,n32=3*n2,n33=3*n3;
  x=(x+n31)%n31;
  y=(y+n32)%n32;
  z=(z+n33)%n33;
  return x+n31*(y+n32*z);
}

/*=========================P=L=A=N=E===========================================/
 *
 *                         Contours lines on a plane 
 *
 *=============================================================================*/

int FourXle::IndexSelectedP( Node& node0, Node& node1, Node& node2, Node& node3 ) {
  if( node1 && node2 && node3 ){
    GLfloat d1 = Distance( Planorte.at(node0.index).vertex, Planorte.at(node1.index).vertex ) + 
      Distance( Planorte.at(node3.index).vertex, Planorte.at(node2.index).vertex );
    GLfloat d2 = Distance( Planorte.at(node0.index).vertex, Planorte.at(node2.index).vertex ) + 
      Distance( Planorte.at(node3.index).vertex, Planorte.at(node1.index).vertex );
    if( d1 > d2 ) return 2; else return 1;
  }else{
    if(      node1 )   return 1;
    else if( node2 )   return 2;
    else if( node3 )   return 3;
  }
  return 0;
}

int cutTriangle(GLfloat va, GLfloat vb, GLfloat vc){
  if ((va<=0.0f)&&(vb<=0.0f)&&(vc<=0.0f)) return 0;
  if ((va>0.0f)&&(vb>0.0f)&&(vc>0.0f)) return 0;

  if ((va<=0.0f)&&(vb>0.0f)&&(vc>0.0f)) return 1; //A
  if ((va>0.0f)&&(vb<=0.0f)&&(vc<=0.0f)) return 1;

  if ((va>0.0f)&&(vb<=0.0f)&&(vc>0.0f)) return 2;
  if ((va<=0.0f)&&(vb>0.0f)&&(vc<=0.0f)) return 2;

  if ((va>0.0f)&&(vb>0.0f)&&(vc<=0.0f)) return 3;
  if ((va<=0.0f)&&(vb<=0.0f)&&(vc>0.0f)) return 3;

  return 4;
}

void FourXle::findContour(QList<V3> &lines, GLfloat value){
  GLfloat Va,Vb,Vc,mix1,mix2;
  V3 sta,end;
  double scope=cScopeBx->value()*cScopeBx->value();
  for (int i=0; i<Planpgns.size();i++){
    int n=Planpgns.at(i).n;
    if (n==3){
      Va=Planorte.at(Planpgns.at(i).ii[0]).color-value;
      Vb=Planorte.at(Planpgns.at(i).ii[1]).color-value;
      Vc=Planorte.at(Planpgns.at(i).ii[2]).color-value;
      switch (cutTriangle(Va, Vb, Vc)){
        case 1: 
          mix1=(Va/(Va-Vb));
          mix2=(Va/(Va-Vc));
          sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
          end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
          if (Distance(sta,end)>1.0) break;
          if (scope > 0.0){
               if (Distance(sta,aufpunkt)> scope) break;
               if (Distance(end,aufpunkt)> scope) break;
          }
          sta+=pnormal*value*heightScale->value();
          end+=pnormal*value*heightScale->value();
          lines.append(sta);
          lines.append(end);
          break;
        case 2: 
          mix1=(Va/(Va-Vb));
          mix2=(Vb/(Vb-Vc));
          sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
          end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[1]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
          if (Distance(sta,end)>1.0) break;
          if (cScopeBx->value() > 0.0){
              if (Distance(sta,aufpunkt)> scope) break;
              if (Distance(end,aufpunkt)> scope) break;
         }
          sta+=pnormal*value*heightScale->value();
          end+=pnormal*value*heightScale->value();
          lines.append(sta);
          lines.append(end);
          break;
        case 3: 
          mix1=(Vc/(Vc-Vb));
          mix2=(Va/(Va-Vc));
          sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[2]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
          end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
          if (Distance(sta,end)>1.0) break;
          if (cScopeBx->value() > 0.0){
               if (Distance(sta,aufpunkt)> scope) break;
               if (Distance(end,aufpunkt)> scope) break;
          }
          sta+=pnormal*value*heightScale->value();
          end+=pnormal*value*heightScale->value();
          lines.append(sta);
          lines.append(end);
          break;
      }
    }else{
      V3 mid=V3(0.0,0.0,0.0);
      for (int j=0; j<n;j++) mid+=Planorte.at(Planpgns.at(i).ii[j]).vertex;
      mid*=1.0/n;
      GLfloat mixColor=0.0f;
      for (int j=0; j<n;j++) mixColor+=Planorte.at(Planpgns.at(i).ii[j]).color;
      mixColor*=1.0f/n;
      Va=mixColor-value;
      for (int k=0;k<n;k++){
        int kp=(k+1)%n;
        Vb=Planorte.at(Planpgns.at(i).ii[k]).color-value;
        Vc=Planorte.at(Planpgns.at(i).ii[kp]).color-value;
        switch (cutTriangle(Va, Vb, Vc)){
          case 1: 
            mix1=(Va/(Va-Vb));
            mix2=(Va/(Va-Vc));
            sta=(1.0f-mix1)*mid + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
            end=(1.0f-mix2)*mid + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
            if (Distance(sta,end)>1.0) break;
            if (cScopeBx->value() > 0.0){
                 if (Distance(sta,aufpunkt)> scope) break;
                 if (Distance(end,aufpunkt)> scope) break;
            }
            sta+=pnormal*value*heightScale->value();
            end+=pnormal*value*heightScale->value();
            lines.append(sta);
            lines.append(end);
            break;
          case 2: 
            mix1=(Va/(Va-Vb));
            mix1=(Va/(Va-Vb));
            mix2=(Vb/(Vb-Vc));
            sta=(1.0f-mix1)*mid + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
            end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[k]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
            if (Distance(sta,end)>1.0) break;
            if (cScopeBx->value() > 0.0){
                 if (Distance(sta,aufpunkt)> scope) break;
                 if (Distance(end,aufpunkt)> scope) break;
            }
            sta+=pnormal*value*heightScale->value();
            end+=pnormal*value*heightScale->value();
            lines.append(sta);
            lines.append(end);
            break;
          case 3: 
            mix1=(Vc/(Vc-Vb));
            mix2=(Va/(Va-Vc));
            sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[kp]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
            end=(1.0f-mix2)*mid + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
            if (Distance(sta,end)>1.0) break;
            if (cScopeBx->value() > 0.0){
                 if (Distance(sta,aufpunkt)> scope) break;
                 if (Distance(end,aufpunkt)> scope) break;
            }
            sta+=pnormal*value*heightScale->value();
            end+=pnormal*value*heightScale->value();
            lines.append(sta);
            lines.append(end);
            break;
          case 4:
            printf("@@~~~~ %g %g %g\n",Va,Vb,Vc);
        }
      }
    }
  }
}

void FourXle::CalcPlaneVertex( int ix, int iy, int iz , V3 n, V3 ap, V3 off){
  GLfloat vo, vx=0, vy=0, vz=0, mix;
  Ort o;  
  //V3 xpo=V3(1,0,0);
  //mole->frac2kart(xpo,xpo);
  V3 mdz=0.25*dx+0.25*dy+0.25*dz;
  vo = n*(dx*ix + dy*iy + dz*iz-ap); //plane n dot (x - x1)
  vx = n*(dx*(ix+1) + dy*iy + dz*iz - ap);
  vy = n*(dx*ix + dy*(iy+1) + dz*iz - ap);
  vz = n*(dx*ix + dy*iy + dz*(iz+1) - ap);
  //V3 fl,mdz=(0.5*dx)+(0.5*dy)+(0.5*dz);
  if( ix != n1-1 && Intersect(vo,vx) ){
    mix=(vo/(vo-vx));
    o.vertex = dx*(mix+ix) + dy*iy + dz*iz +off;//+urs;
/*
    mole->kart2frac(o.vertex,o.vertex);
    o.vertex+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
    o.vertex+=-1.0*fl;
    o.vertex+=V3(0.5,0.5,0.5);
    mole->frac2kart(o.vertex,o.vertex);
    o.vertex+=-1.0*urs;
// */
    o.vertex+=mdz;
    //o.vertex+=test3;
    //o.vertex+=-1.0*xpo;
    o.normal = n;
    if (sourceMap->currentIndex()==0)
        o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[((1+ix)%n1)+(iy%n2)*n1+(iz%n3)*n4];
    else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[((1+ix)%n1)+(iy%n2)*n1+(iz%n3)*n4];
    Planorte.append(o);
    nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index=Planorte.size()-1;
    nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
  }else{
    nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag   = 0;
  }
  if( iy != n2-1 && Intersect(vo,vy) ){
    mix=(vo/(vo-vy));
    o.vertex = dx*ix + dy*(mix+iy) + dz*iz+off;//+urs;
/*

    mole->kart2frac(o.vertex,o.vertex);
    o.vertex+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
    o.vertex+=-1.0*fl;
    o.vertex+=V3(0.5,0.5,0.5);
    mole->frac2kart(o.vertex,o.vertex);
    o.vertex+=-1.0*urs;
// */
    o.vertex+=mdz;
    //o.vertex+=test3;

    //o.vertex+=-1.0*xpo;
    o.normal = n;
    if (sourceMap->currentIndex()==0)
        o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+ mix*datfo_fc[(ix%n1)+((iy+1)%n2)*n1+(iz%n3)*n4];
    else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+ mix*datfo_fc[(ix%n1)+((iy+1)%n2)*n1+(iz%n3)*n4];
    Planorte.append(o);
    nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index=Planorte.size()-1;
    nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
  }else{
    nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=0;
  }
  if( iz != n3-1 && Intersect(vo,vz) ){
    mix=(vo/(vo-vz));
    o.vertex = dx*ix + dy*iy + (mix+iz)*dz+off;//+urs;
/*
    mole->kart2frac(o.vertex,o.vertex);
    o.vertex+=V3(-0.5,-0.5,-0.5);
    fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
    o.vertex+=-1.0*fl;
    o.vertex+=V3(0.5,0.5,0.5);
    mole->frac2kart(o.vertex,o.vertex);
    o.vertex+=-1.0*urs;
//   */

    o.vertex+=mdz;
   // o.vertex+=test3;
    //o.vertex+=-1.0*xpo;
    o.normal = n;
    if (sourceMap->currentIndex()==0)
        o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[(ix%n1)+(iy%n2)*n1+((iz+1)%n3)*n4];
    else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[(ix%n1)+(iy%n2)*n1+((iz+1)%n3)*n4];
    Planorte.append(o);
    nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index = Planorte.size()-1;
    nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
  }else{
    nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag  = 0;
  }
}

void FourXle::makeFacesp(int nn, Node poly[] ){
  //   /* 
  int n=0;
  int ly[13];
  for (int j=0;j<nn;j++)//zu nahe beieinander liegende verts ignorieren
    if (Distance(Planorte.at(poly[j].index).vertex,Planorte.at(poly[(j+1)%nn].index).vertex)>0.00)
      ly[n++]=poly[j].index;
    else{
      Planorte[poly[j].index].vertex=
        Planorte[poly[(j+1)%nn].index].vertex=
        (Planorte.at(poly[j].index).vertex+Planorte.at(poly[(j+1)%nn].index).vertex)*0.5f;
      Planorte[poly[j].index].color=
        Planorte[poly[(j+1)%nn].index].color=
        (Planorte.at(poly[j].index).color+Planorte.at(poly[(j+1)%nn].index).color)*0.5f;
      ly[n++]=poly[(j+1)%nn].index;
      j++;
    }
  GLfloat mid_col = 0.0;
  V3 mid_ver = V3(0,0,0);
  V3 mid_nor = V3(0,0,0);
  for(int i=0; i<n; i++ ){
    mid_ver += Planorte.at(ly[i]).vertex;
    mid_nor += Planorte.at(ly[i]).normal;
    mid_col += Planorte.at(ly[i]).color;
  }
  mid_ver *= (1.0/n);
  Normalize(mid_nor);
  mid_col *= (1.0/n);
  if (n<3) return;  //weniger als 3 verts -> nichts zu tun
  Polygn neupoly;
  neupoly.vertex=mid_ver;
  neupoly.n=n;
  //if (Planorte.at(poly[0].index).direct*iso_level<0.0){//drehsinn der polygone abfragen* /
    for(int i=0; i<n; i++ ){
      neupoly.ii[i]=ly[i];
    }
  //}
  /*else{
    for(int i=n-1,k=0; i>=0; i-- ){
      neupoly.ii[i]=ly[k++];
    }
  }//else*/
  Planpgns.append(neupoly);
  //  */
}

void FourXle::MakeElementp( int ix, int iy, int iz ,int s1, int s2) {//das ist der Teil des japanischen Programms den ich nicht verstehe.
  //Hauptsache fuktioniert.
  //    /*
  static int conn[12][2][4] = {           //char->int wegen warning g++ >3.0
    {{ 0, 1, 7, 6}, { 0, 2, 8, 3}},  //  0
    {{ 1, 2, 5, 4}, { 1, 0, 6, 7}},  //  1
    {{ 2, 0, 3, 8}, { 2, 1, 4, 5}},  //  2
    {{ 3, 8, 2, 0}, { 3, 4,10, 9}},  //  3
    {{ 4, 3, 9,10}, { 4, 5, 2, 1}},  //  4
    {{ 5, 4, 1, 2}, { 5, 6, 9,11}},  //  5
    {{ 6, 5,11, 9}, { 6, 7, 1, 0}},  //  6
    {{ 7, 6, 0, 1}, { 7, 8,11,10}},  //  7
    {{ 8, 7,10,11}, { 8, 3, 0, 2}},  //  8
    {{ 9,10, 4, 3}, { 9,11, 5, 6}},  //  9
    {{10,11, 8, 7}, {10, 9, 3, 4}},  // 10
    {{11, 9, 6, 5}, {11,10, 7, 8}}   // 11
  };
  static Node node[12];
  static Node polygon[12];

  node[ 0] = nodeX[(ix+iy*s1+iz*s2)                  %n5];// 000x
  node[ 1] = nodeY[(ix+iy*s1+iz*s2)                  %n5];// 000y
  node[ 2] = nodeZ[(ix+iy*s1+iz*s2)                  %n5];// 000z
  node[ 3] = nodeX[(ix+iy*s1+((iz+1)%n3)*s2)         %n5];// 001y
  node[ 4] = nodeY[(ix+iy*s1+((iz+1)%n3)*s2)         %n5];// 001z
  node[ 5] = nodeZ[(ix+((iy+1)%n2)*s1+iz*s2)         %n5];// 010x
  node[ 6] = nodeX[(ix+((iy+1)%n2)*s1+iz*s2)         %n5];// 010y
  node[ 7] = nodeY[(((1+ix)%n1)+iy*s1+iz*s2)         %n5];// 100y
  node[ 8] = nodeZ[(((1+ix)%n1)+iy*s1+iz*s2)         %n5];// 100z
  node[ 9] = nodeX[(ix+((iy+1)%n2)*s1+((iz+1)%n3)*s2)%n5];// 011x
  node[10] = nodeY[(((ix+1)%n1)+iy*s1+((iz+1)%n3)*s2)%n5];// 101y
  node[11] = nodeZ[(((ix+1)%n1)+((iy+1)%n2)*s1+iz*s2)%n5];// 110z

  if (((char)node[0]+node[1]+node[2]+node[3]+node[4]+node[5]
        +node[6]+node[7]+node[8]+node[9]+node[10]+node[11])==0) return;
  for( int is=0; is<12; is++ ) {
    if( !node[is] ) continue;

    int n=0, i=is, m=0;//,ai=i;
    do {
      polygon[n++]= node[i];
      int sol = IndexSelectedP(node[conn[i][m][0]],node[conn[i][m][1]],
          node[conn[i][m][2]],node[conn[i][m][3]]);
      //   ai=i;
      i = conn[i][m][sol];

      if( sol == 2 ) m ^= 1;
      node[i].flag = 0;
    }while( (i!=is)&&(n<11) );
      {
        makeFacesp( n, polygon );
      }
  }
  //  */
}

void FourXle::makePlane(QList<V3> &lines,int a1, int a2, int a3) {
  /*
  farbe[0][0]=0.6;    
  farbe[0][1]=0.0;    
  farbe[0][2]=0.0;    
  farbe[0][3]=1.0;    

  farbe[1][0]=0.6;    
  farbe[1][1]=0.0;    
  farbe[1][2]=0.0;    
  farbe[1][3]=1.0;    

  farbe[2][0]=0.6;    
  farbe[2][1]=0.0;    
  farbe[2][2]=0.0;    
  farbe[2][3]=1.0;    

  farbe[3][0]=0.0;    
  farbe[3][1]=0.0;    
  farbe[3][2]=0.0;    
  farbe[3][3]=1.0;    

  farbe[4][0]=0.0;    
  farbe[4][1]=0.0;    
  farbe[4][2]=1.0;    
  farbe[4][3]=1.0;

  farbe[5][0]=0.0;    
  farbe[5][1]=0.0;    
  farbe[5][2]=1.0;    
  farbe[5][3]=1.0;

  farbe[6][0]=0.0;    
  farbe[6][1]=0.0;    
  farbe[6][2]=1.0;    
  farbe[6][3]=1.0;

  Farben=7;
  */

//  fprintf(stderr, "makePlane %d %d %d %d %d\n",n1,n2,n3,n4,n5);
  //orig=V3(0,0,0);
//  n4=n1*n2;
//  data.clear();

  if (( nodeX =(Node*)malloc(sizeof(Node)*n5))==NULL) {
    qDebug()<<"Less Memory(X) ";
    exit(1);  
  } 
  if (( nodeY =(Node*)malloc(sizeof(Node)*n5))==NULL) {
    qDebug()<<"Less Memory(Y)!!";
    exit(1); 
  }
  if (( nodeZ =(Node*)malloc(sizeof(Node)*n5))==NULL) {
    qDebug()<<"Less Memory(Z)!!";
    exit(1); 
  }
//  printf("UIUI %p %p %p\n",nodeX,nodeY,nodeZ);
  for (int ix=0; ix<n5; ix++){
    nodeX[ix].index=0;
    nodeX[ix].flag=0;
    nodeY[ix].index=0;
    nodeY[ix].flag=0;
    nodeZ[ix].index=0;
    nodez[ix].flag=0;
  }
//  printf("__\n");
  Planorte.clear();
  Planpgns.clear();
  lines.clear();
  chgl->contval.clear(); 

  //test3= ((n1-1)/-2.0) *  dx + ((n2-1)/-2.0) * dy + ((n3-1)/-2.0) * dz + orig;
  //test3=orig;
//  printf("Atomindices: %d %d %d < %d\n",a1,a2,a3,mole->showatoms.size());
  if ((a1>=mole->showatoms.size())||(a2>=mole->showatoms.size())||(a3>=mole->showatoms.size())) return;
  if ((a1<0)||(a2<0)||(a3<0)) return;
  if  (mole->showatoms[a1].pos==mole->showatoms[a2].pos){
    mole->frac2kart(mole->showatoms[a1].frac,mole->showatoms[a1].pos);
    mole->frac2kart(mole->showatoms[a2].frac,mole->showatoms[a2].pos);
    mole->frac2kart(mole->showatoms[a3].frac,mole->showatoms[a3].pos);
  }
  V3 a1v=V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z);
  V3 a2v=V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z);
  V3 a3v=V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z);
  pnormal = Normalize((a2v - a1v) % (a3v - a1v));
  switch (centerIsOn->currentIndex()){
    case 0: aufpunkt = V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z);break;
    case 1: aufpunkt = V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z);break;
    case 2: aufpunkt = V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z);break;
    case 3: {
              aufpunkt = (1.0/3.0)*
               (V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z)+
                V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z)+
                V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z));
              break;
            }
    
  }
//  aufpunkt = V3(mole->showatoms.at(a1).kart.x,mole->showatoms.at(a1).kart.y,mole->showatoms.at(a1).kart.z);
//  printf("%g  %g %g %g   %g %g %g\n",Norm(pnormal),pnormal.x,pnormal.y,pnormal.z,aufpunkt.x,aufpunkt.y,aufpunkt.z);
  double angle = winkel(pnormal, V3(0,0,1));
  V3 ax= Normalize(pnormal % V3(0,0,1));
  double cp = cos(angle), cp1 = 1.0 - cos(angle), sp = sin(angle);
  double X=ax.x,Y=ax.y,Z=ax.z;
  Matrix ROT = Matrix(
         cp+cp1*X*X, -Z*sp+cp1*X*Y,  Y*sp+cp1*X*Z, 
       Z*sp+cp1*X*Y,    cp+cp1*Y*Y, -X*sp+cp1*Y*Z, 
      -Y*sp+cp1*X*Z,  X*sp+cp1*Y*Z,    cp+cp1*Z*Z);
  V3 bv=ROT*Normalize(V3(a2v.x-a1v.x,a2v.y-a1v.y,a2v.z-a1v.z));
  //printf("BV %g %g %g\n",bv.x,bv.y,bv.z);
  double bww=winkel(bv,V3( 1,0,0));
  cp=cos(bww);
  sp=sin(bww);
  //printf("%g %g %g\n",bww*180.0/M_PI,cp,sp);
  Matrix R2 = Matrix(
      cp, sp, 0,
     -sp, cp, 0,
       0,  0, 1);//*/
  bv=R2*bv;
  bww=winkel(bv,V3( 1,0,0));
  if (fabs(bww)>1) {
  R2 = Matrix(
      cp,-sp, 0,
      sp, cp, 0,
       0,  0, 1);//*/
  
  } 
  //printf("BV %g %g %g %g\n",bv.x,bv.y,bv.z,bww);
  V3 vt = ROT * V3(pnormal.x, pnormal.y, pnormal.z);
  vt = R2*vt;
  //printf("ang %g %g %g \n%12.6f %12.6f %12.6f \n%g %g\n%g\n",angle/M_PI*180.0, pnormal*ax,ax*V3(0,0,1), vt.x, vt.y, vt.z,determinant(ROT),Norm(ax),winkel(a1v,V3(0,1,0)));
//  fprintf(stderr, "makePlane\n");
  int qii=0;V3 ap;
  int aa=-1,ee=2;
  if (!extendUnit->isChecked()){aa++;ee--;}
  for (int iix=aa; iix<ee; iix++){
     for (int iiy=aa; iiy<ee; iiy++){
        for (int iiz=aa; iiz<ee; iiz++){
            V3 propo=V3(iix,iiy,iiz)+V3(0.5,0.5,0.5);
            V3 pp=V3(iix,iiy,iiz);
            mole->frac2kart(propo,propo);
            mole->frac2kart(pp,pp);
            double abst=pnormal*(propo-aufpunkt);
            V3 da=propo-abst*pnormal;
            mole->kart2frac(da,da);
            qii++;
           // printf("PLANE %2d %2d %2d %12g %12g %12g  %g  %12g %12g %12g\n",iix,iiy,iiz,da.x,da.y,da.z,abst,da.x-iix,da.y-iiy,da.z-iiz);
            ap=da-V3(iix,iiy,iiz);
            mole->frac2kart(ap,ap);
            //printf("H%d' 1 %8.5f %8.5f %8.5f 11.0 0.05 %5.3f \n",qii, da.x, da.y, da.z,abst);
  //          if ((da.x-iix<1.0)&&(da.x-iix>0.0)&&(da.y-iiy<1.0)&&(da.y-iiy>0.0)&&(da.z-iiz<1.0)&&(da.z-iiz>0.0)){
            for( int ix=0; ix<=n1; ix++ ){
              for( int iy=0; iy<=n2; iy++ ){
                for( int iz=0; iz<=n3; iz++ ){
                  CalcPlaneVertex(ix,iy,iz,pnormal,ap,pp);//
                  //CalcPlaneVertex(ix,iy,iz,pnormal,aufpunkt,V3(0,0,0));
                }
              }
            }
            //printf("%3d %3d %3d Planorte.size = %d\n",iix,iiy,iiz,Planorte.size());
            //fprintf(stderr, "makePlane\n");
            for( int ix=0; ix<n1-1; ix++ )
              for( int iy=0; iy<n2-1; iy++ )
                for( int iz=0; iz<n3-1; iz++ )
                  MakeElementp(ix,iy,iz,n1,n4);
//}
        }
     }
  }
  double minP=1e37,maxP=-1e38;
  for (int i=0; i<Planorte.size(); i++) {
    minP=fmin(Planorte.at(i).color,minP);
    maxP=fmax(Planorte.at(i).color,maxP);
  }  
  //fprintf(stderr, "makePlane Planpgns.size = %d Planorte.size = %d minP = %g maxP = %g\n",Planpgns.size(),Planorte.size(),minP,maxP);
  min= 1e37;
  max=-1e37;
  QStringList numbers = QString(contourValueEdit->text()).split(" ",skipEmptyParts);
  int z=0;
 // printf("ok\n");
  if (!numbers.isEmpty()) chgl->contval[0]=numbers.at(0).toFloat();
  for (int i = 0; i < numbers.size(); i++){
    float contour = numbers.at(i).toFloat();
    findContour(lines,contour);
    chgl->contval[lines.size()]=contour;
    if ((lines.size()-z)>0){
      min=fmin(contour,min);
      max=fmax(contour,max);
    }
    z=lines.size();
  }
  chgl->contMax=max;
  chgl->contMin=min;
 // printf("ok %f %f\n",min,max);
  //legende();
  free(nodeX);
  nodeX=NULL;
  free(nodeY);
  nodeY=NULL;
  free(nodeZ);
  nodeZ=NULL;
  Planpgns.clear();
  Planorte.clear();
  //printf("ok\n");
  QList<V3> lin2d;
  int aspect1=3,aspect2=4;//todo make aspect ratio changable
  switch (aspectRatios->currentIndex()) {
    case 0: aspect1=3,aspect2=4; break;
    case 1: aspect1=9,aspect2=16; break;
    case 2: aspect1=1,aspect2=1; break;
  }
  int width=1024*aspect1/aspect2;
  double miX =  9.e30;
  double maX = -9.e30;
  double miY =  9.e30;
  double maY = -9.e30;
  V3 v1;//,v2=V3(99,99,99);
  for (int i=0; i<lines.size(); i++){
    v1=R2*(ROT * V3(lines.at(i).x-aufpunkt.x,lines.at(i).y-aufpunkt.y,lines.at(i).z-aufpunkt.z));
    //if (Distance(v1,v2)>0.0) 
    lin2d.append(v1);
    miX=fmin(v1.x,miX);
    maX=fmax(v1.x,maX);
    miY=fmin(v1.y,miY);
    maY=fmax(v1.y,maY);


    //v2=v1;
  }
  //printf("!%d == %d! %f < X < %f, %f < Y < %f \n",lin2d.size(),lines.size(),miX,maX,miY,maY);
  if (cScopeBx->value() > 0.0){
    double XX=cScopeBx->value();
    double YY=(XX*aspect1)/aspect2;//todo make aspect ratio changable
    miX=-XX;
    miY=-YY;
    maX= XX;
    maY= YY;
  }
  //printf("!%d == %d! %f < X < %f, %f < Y < %f \n",lin2d.size(),lines.size(),miX,maX,miY,maY);
  if (contEPSFile->text().isEmpty()){ 
   // qDebug()<<"contEPSFile->text() is Empty()";
      QString contourDescription1;
      for (int i=0; i<lin2d.size()/2; i++){
        //(wrt-min)/(max-min)
        if (chgl->contval.contains(i*2)) {
          if (chgl->contval.value(i*2)>0.00001f) {
            contourDescription1.append(QString("<font color=blue>blue contours at %1 e&Aring;<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
          }
          else if (chgl->contval.value(i*2)<-0.00001f){
            contourDescription1.append(QString("<font color=red>red contours at %1 e&Aring;<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
          }
          else {
            contourDescription1.append(QString("<font color=black>black contours at 0 e&Aring;<sup>-3</sup></font> <br>\n"));
          }
        }
      }
      emit bigmessage(QString("%1<br>").arg(contourDescription1));
    return;
  }

//      ||(!QFile::exists (contEPSFile->text()))) 
//qDebug()<<QDate::currentDate().toString(Qt::ISODate);
  FILE *f=fopen(contEPSFile->text().toLocal8Bit(),"wt");
  fprintf(f,"%s%d\n","%!PS-Adobe-3.0 EPSF-3.0\n%%BoundingBox: 0 0 1024 ",width);
  fprintf(f,"%s",
"%%Title: %s Contour Plot\n"
"%%Creator: %s by Christian B. Hubschle\n",PROGRAM_NAME,PROGRAM_NAME);
  fprintf(f,"%%%%CreationDate: %s\n",QDate::currentDate().toString(Qt::ISODate).toStdString().c_str());

  fprintf(f,"%s",
//"%%Pages: 1\n"
"%%EndComments\n"
"%%BeginProlog\n"
//"/l { newpath moveto lineto stroke } bind def\n"
//"/nm { newpath moveto } bind def\n"
"/cp {closepath} bind def\n"
"/ef {eofill} bind def\n"
"/gr {grestore} bind def\n"
"/gs {gsave} bind def\n"
"/sa {save} bind def\n"
"/rs {restore} bind def\n"
"/l {lineto} bind def\n"
"/m {moveto} bind def\n"
"/rm {rmoveto} bind def\n"
"/n {newpath} bind def\n"
"/s {stroke} bind def\n"
"/sh {show} bind def\n"
"/slc {setlinecap} bind def\n"
"/slj {setlinejoin} bind def\n"
"/slw {setlinewidth} bind def\n"
"/srgb {setrgbcolor} bind def\n"
"/rot {rotate} bind def\n"
"/sc {scale} bind def\n"
"/sd {setdash} bind def\n"
"/ff {findfont} bind def\n"
"/sf {setfont} bind def\n"
"/scf {scalefont} bind def\n"
"/sw {stringwidth} bind def\n"
//"%%EndSetup\n"
"%%EndProlog\n"
//"%%Page: 1 1\n"
//"gsave\n"
//"20 20 scale\n"
//"10 10 translate\n"
"save\n0 slc 0 slj\n/Helvetica ff 18 scf sf\n");
fprintf(f,"n 0 0 m 1024 0 l 1024 %d l 0 %d l cp clip s\n",width,width);
fprintf(f,"n 2 2 m 1023 2 l 1023 %d l 2 %d l cp s \n",width-2,width-2);
//"0 setgray\n"


  //Atoms
  for (int i=0; i<mole->showatoms.size(); i++){
    V3 ato=R2*(ROT*(mole->showatoms.at(i).pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
    if ((ato.x>maX)||(ato.x<miX)||(ato.y>maY)||(ato.y<miY)) continue;
    //Acol[xdinp[j].OrdZahl]
    if (fabs(ato.z)<0.2){
      double rad=(mole->showatoms.at(i).an>-1)?mole->arad[mole->showatoms.at(i).an]:0.15;
      QColor col=mole->AtomColor[mole->showatoms.at(i).an];
      fprintf(f,"n %G %G %G 0 360 arc cp gs %g %g %g srgb fill s gr%% %g %s\n",
          (ato.x-miX)/(maX-miX)*1024, (ato.y-miY)/(maY-miY)*width, rad*50,
          col.redF(), 
          col.greenF(), 
          col.blueF(),ato.z,
          mole->showatoms.at(i).Label.toStdString().c_str());
    }
    //    4 5 3 0 360 arc closepath
  }
  //Bonds
  for (int k=0;k<mole->showbonds.size();k++){
    V3 anf=R2*(ROT*(mole->showbonds[k].ato1->pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
    V3 end=R2*(ROT*(mole->showbonds[k].ato2->pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
    if ((anf.x>maX)||(anf.x<miX)||(anf.y>maY)||(anf.y<miY)) continue;
    if ((end.x>maX)||(end.x<miX)||(end.y>maY)||(end.y<miY)) continue;
    if ((fabs(anf.z)<0.2)&&(fabs(end.z)<0.2))
    fprintf(f,"n %G %G m %G %G l cp s\n",
        (anf.x-miX)/(maX-miX)*1024, (anf.y-miY)/(maY-miY)*width,
        (end.x-miX)/(maX-miX)*1024, (end.y-miY)/(maY-miY)*width);

  }

 // V3 v1,v2;
 //qDebug()<<contval;
  float red=0.0f,green=0.0f,blue=0.0f;
  QString contourDescription;
  for (int i=0; i<lin2d.size()/2; i++){
    //(wrt-min)/(max-min)
    if (chgl->contval.contains(i*2)) {
      if (chgl->contval.value(i*2)>0.00001f) {
        red=0.0f;green=0.0f;blue=1.0f;
        contourDescription.append(QString("% blue@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
//        printf("blau %f %d\n",contval.value(i*2),i);
      }
      else if (chgl->contval.value(i*2)<-0.00001f){
        red=1.0f;green=0.0f;blue=0.0f;
        contourDescription.append(QString("% red@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
//        printf("rot %f %d\n",contval.value(i),i);
      }
      else {
        red=0.0f;green=0.0f;blue=0.0f;
        contourDescription.append(QString("% black@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
//        printf("schwarz %f %d\n",chgl->contval.value(i*2),i);
      }
    }
    if ((lin2d.at(2*i  ).x>maX)||(lin2d.at(2*i  ).x<miX)||(lin2d.at(2*i  ).y>maY)||(lin2d.at(2*i  ).y<miY)) continue;
    if ((lin2d.at(2*i+1).x>maX)||(lin2d.at(2*i+1).x<miX)||(lin2d.at(2*i+1).y>maY)||(lin2d.at(2*i+1).y<miY)) continue;
    fprintf(f,"n %G %G m %G %G l cp gs %g %g %g srgb 0.5 slw s gr\n",
        (lin2d.at(2*i  ).x-miX)/(maX-miX)*1024, (lin2d.at(i*2  ).y-miY)/(maY-miY)*width, 
        (lin2d.at(2*i+1).x-miX)/(maX-miX)*1024, (lin2d.at(i*2+1).y-miY)/(maY-miY)*width,red,green,blue);
  }
  for (int i=0; i<mole->showatoms.size(); i++){
    V3 ato=R2*(ROT*(mole->showatoms.at(i).pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
    if ((ato.x>maX)||(ato.x<miX)||(ato.y>maY)||(ato.y<miY)) continue;
    if (fabs(ato.z)<0.2){
      double rad=(mole->showatoms.at(i).an>-1)?mole->arad[mole->showatoms.at(i).an]:0.15;
      fprintf(f,"gs n %G %G m (%s) true charpath 1 slw 1 setgray s gr\n",(ato.x-miX)/(maX-miX)*1024+rad*50, (ato.y-miY)/(maY-miY)*width,mole->showatoms.at(i).Label.toStdString().c_str());
      fprintf(f,"gs %G %G m (%s) sh gr\n",(ato.x-miX)/(maX-miX)*1024+rad*50, (ato.y-miY)/(maY-miY)*width,mole->showatoms.at(i).Label.toStdString().c_str());

    }
    //    4 5 3 0 360 arc closepath    
  }
  fprintf(f,"%s",
      "restore\n"
      "showpage\n"      
      "%%Trailer\n"
      "%%EOF\n"
      );
  fprintf(f,"%% plane: %s %s %s\n%% contours: %s\n",
      mole->showatoms.at(a1).Label.toStdString().c_str(),
      mole->showatoms.at(a2).Label.toStdString().c_str(),
      mole->showatoms.at(a3).Label.toStdString().c_str(),
      contourDescription.toStdString().c_str());
  //printf("->%d\n",fclose(f));
  //printf("ok3074\n");
  QString contourDescription1;
  for (int i=0; i<lin2d.size()/2; i++){
    //(wrt-min)/(max-min)
    if (chgl->contval.contains(i*2)) {
      if (chgl->contval.value(i*2)>0.00001f) {
        contourDescription1.append(QString("<font color=blue>blue contours at %1 e&Aring;<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
      }
      else if (chgl->contval.value(i*2)<-0.00001f){
        contourDescription1.append(QString("<font color=red>red contours at %1 e&Aring;<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
      }
      else {
        contourDescription1.append(QString("<font color=black>black contours at 0 e&Aring;<sup>-3</sup></font> <br>\n"));
      }
    }
  }
  emit bigmessage(QString("%1<br>").arg(contourDescription1));

}