/****************************************************************************
 **
 ** 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.
 **
 **
 ****************************************************************************/
#ifndef FOURXLE_H
#define FOURXLE_H 1
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <locale.h>
//#include <fftw3.h> //if you want you can try to uncoment this, but then you have to coment the following line
#include "kissfft/kiss_fftnd.h"
#include "molecule.h"
#include "chgl.h"
#include <QObject>
#include <QCheckBox>
#include <QMutex>
#include <QMutexLocker>

//! Rec is a reflection type of a fcf 6 file.
typedef struct {
  int  ih,//!< h
       ik,//!< k
       il;//!< l 
  float fo,//!< F observed
        so,//!< \f$\sigma(observed)\f$
        fc,//!< F calculated
        phi; //!< \f$\varphi\f$ 
} Rec;
#define LM 2000000
//! FNode is a struct for the 3 different edges of a cube 
struct FNode {
  V3 vertex;//!< vertex is the 3d position of the node
  V3 normal;//!< normal is the plane normal of the surface
  char flag;//!< this flag is set when the surface crosses this node
  inline operator char (){
    //! \returns the flag parameter. No Coala ;-)
    return flag;
  }
  inline operator V3(){
    //! \returns the vertex of the node.
    return vertex;
  }
  friend inline double Distance( const FNode& n1, const FNode& n2 ){
    /*! @param n1,n2 nodes vor which the squared distance should be calculated.
      \returns te positional squared distance between two nodes
      */
    return Norm(n1.vertex-n2.vertex);
  }
};

struct Polygn {
  V3 vertex;
  int n;
  int ii[13];
};

struct Ort {
  V3 vertex;
  V3 normal;
  GLfloat color;
  char direct;
};

struct Node {
  int index;
  unsigned char flag; 
  inline operator unsigned char (){return flag;}
};

//
//!  FourXle is an QObject that loads a fcf 6 file and creates iso-surfaces in meshed style \inherit QObject 
class FourXle:public QObject{
  Q_OBJECT
  public:
    int HKLMX;//!< Maximum of h,k and l values can be used to make a resulution cut. 
    float *datfo,//!<data pointer for the fobs map in real space
          *datfo_fc;//!<data pointer for the fobs-fcalc map in real space
    float f000;
    float overridef000;
    bool isBelo;
    FourXle(Molecule *mole_, ChGL *chgl_,QToolBar *toolView, double resol=2.5, double wght=1.1);
    /*!< @param mole_ pointer to the Molecule of ShelXle 
     *  @param chgl_ pointer to the ChGL widged of ShelXle
     *  @param toolView pointer to the View Toolbar of ShelXle
     *  @param resol resolution factor (the higher the more grid points)
     *  @param wght a parameter to downweight weak reflections
     */
    ~FourXle();
    QDialogButtonBox *buttonBoxMC;
    QPushButton *applyMC;          
    QPushButton *applyDF;          
    QPushButton *jnkButton;
    QLabel *wl, *pl, *dl, *ol, *rl, *sl, *lw, *ltz;
    QDialog *md;
    bool MYFCF;
    QDoubleSpinBox *maprad,*fomaps,*difmaps,*weak,*mapprec,*lineTrans,*lineWidth,*cutrange;
    QPushButton *foPlusButton,*foMinusButton,*diffPlusButton,*diffMinusButton, *defColorButton;
    QComboBox *mapSchnitt;
    QGridLayout *mdl;
    QVBoxLayout *amdl;
    QAction *keepiso;
    QAction *mapcontrol;
    QCheckBox *ownfcf6, *ownfcf6_;
    QString  fouName;
    void drawMap(int faps);
    int maptrunc;
    double map_radius; //!< radius of the map in Angstrom around the rotaion center (used only in a mode of FourXle).

    bool loadFouAndPerform(const char filename[],bool neu=true);
    bool loadMap(int NX, int NY, int NZ, float *dat);
    QDir resDir;
    QCheckBox *keepIso;//!<keeps the initial iso values good for idieal presentations
    float sigma[3];//!<sigma values 
    float iso[3];//!<iso values
    //  double mInimum[2],mAximum[2];
    V3 urs;//!< origin
    double rr,//!< resolution factor (the higher the more grid points) 
           rw;//!< a parameter to downweight weak reflections

    QCheckBox *doMaps;//!< do or notthing todo is the question
    Molecule *mole;//!< pointer to the Molecule object
    ChGL *chgl;//!< pointer to the ChGL widget
    QString voxelstr;//!< a string to fix dimensions to desired values by comand line option '-voxel 92x92x92 '
    //void deleteLists();
    void killmaps();
    void jnk();
    float sigmaFo();
    float sigmaFoFc();
    void exportMaps(int na, const char filename[], const char atomlist[]);
    int n1,//!< dimension of the map in a diRection
        n2,//!< dimension of the map in b diRection
        n3,//!< dimension of the map in c diRection
        n4;//!< \f$ n4 = n1 \times n2\f$
    int n5;//!< \f$ n5 = n1 \times n2 \times n3\f$
    V3 dx,//!< vector in a direction for each map voxel
       dy,//!< vector in b direction for each map voxel
       dz;//!< vector in c direction for each map voxel
  QVector<Ort> Planorte;
  QVector<Polygn> Planpgns;
//  QList<double> data;
  inline int Intersect( GLfloat& vm, GLfloat& vp ){
    return vm*vp <= 0.0 && (vm<0.0 || vp<0.0);
  }
  Node *nodeX;
  Node *nodeY;
  Node *nodeZ;
  QCheckBox *extendUnit;
  void CalcPlaneVertex( int ix, int iy, int iz ,V3 n,V3 ap,V3 off);
  void makeFacesp(int nn, Node poly[] );
  void MakeElementp( int ix, int iy, int iz ,int s1, int s2);
  int IndexSelectedP( Node& node0, Node& node1, Node& node2, Node& node3 );
  void findContour(QList<V3> &lines, GLfloat value);
  void makePlane(QList<V3> &lines,int a1, int a2, int a3);
  double min,max;
  V3 orig,test3;
  V3 pnormal,aufpunkt;
  QLineEdit *contMapFile, *contourValueEdit, *contEPSFile;
  QComboBox *aspectRatios, *centerIsOn, *sourceMap;
  QDoubleSpinBox *cScopeBx, *heightScale;
  inline double winkel(V3 a,V3 b){
    return acos((a.x*b.x+a.y*b.y+a.z*b.z)/(sqrt(a.x*a.x+a.y*a.y+a.z*a.z)*sqrt(b.x*b.x+b.y*b.y+b.z*b.z)));
  }
    public slots:
      void bewegt(V3 v);
    void inimap();
    void change_iso(int numsteps,int diff);
    void chglDestroyed();
    void foactDestroyed();
    void fofcactDestroyed();
    void doMapsNow(bool b);//!<destroys maps or (re)calculates them freshly.
    void mapDefault();//!< restores the defaults of the Map Control dialog.
    void schnittart(int trunc);//!< the truncation type of the electron density maps changes to trunc @param trunc new truncation type.
    void controlMap();//!<Shows the Map-Control Dialog.
    void keepIsoValues(bool b);
    void killFourXleFirst();//!< attempt to run killmaps before chgl gets destryed
    void colorDLGFOP();//!<Color Dialog for Fobs plus. Positive values of Fo map are drawn with a mesh in this color.
    void colorDLGFOM();//!<Color Dialog for Fobs minus. Negative values of Fo map are drawn with a mesh in this color.
    void colorDLGDIP();//!<Color Dialog for Fobs-Fcalc plus. Positive values of Fo-Fc map are drawn with a mesh in this color.
    void colorDLGDIM();//!<Color Dialog for Fobs-Fcalc minus. Negative values of Fo-Fc map are drawn with a mesh in this color.
    void defaultColors();//!<Restores the default colors of the electron density maps (Fo+ blue Fo- yellow Fo-Fc+ green Fo-Fc- red).
    void openMapControl();//!< shows the Map Control dialog.
    void updateLines();
    void reportInfo();
signals:
    void bigmessage(const QString &);//!< to comunicate with the informaton window
    void reloadRes();
    void haltMW();//stop the map worker

  private:
    mutable QMutex mutex;
    double C[15],D[9],sy[12][192],wave;
    inline int Intersect( double& vm, double& vp ){ return vm*vp <= 0.0 && (vm<0.0 || vp<0.0); }
    inline int dex(int x,int y,int z);
    inline int dex3(int x,int y,int z);
    int oldatomsize;
    int acnt;
    QString info;
    V3 oc,dxc,dyc,dzc;
    V3  delDA[27];
#ifdef FFTW3_H
    fftwf_plan  fwd_plan;//!!!
    fftwf_complex *B;//!!!
#else
    kiss_fftnd_cfg fwd_plan;//!!!
    kiss_fft_cpx *B;//!!!
#endif
    FNode *nodex,*nodey,*nodez;
    int *noGoMap;
    Rec lr[LM];
    char cen,git;
    int nr,nc,ns;
    int mtyp,tri;
    QTime foti,suti;
    //          double minnor,maxnor;
    void gen_surface(bool neu,int imin=0, int imax=4);
    void CalcVertex( int ix, int iy, int iz);
    //    V3 CalcNormalX( int ix, int iy, int iz );
    //    V3 CalcNormalY( int ix, int iy, int iz );
    //    V3 CalcNormalZ( int ix, int iy, int iz );
    int IndexSelected( FNode& node0, FNode& node1, FNode& node2, FNode& node3);
    //    V3& VectorSelected( FNode& node0, FNode& node1, FNode& node2, FNode& node3);
    void MakeElement( int ix, int iy, int iz, int faps);
    void makeFaces(int n, int faps, FNode poly[] );
    char titl[80];/*fcmax=0,f000=0,resmax=99999.0,*/
    void trimm(char s[]);
    void deletes(char *s, int count);
    int readHeader(const char *filename);
    void sorthkl(int nr, Rec r[]);

    friend class GenMapWorker;
};



class GenMapWorker : public QThread{
  Q_OBJECT
  public:
    GenMapWorker(ChGL *chgl_,Molecule *mol_, FourXle *parent, bool neu_=true,int imin_=0, int imax_=4){
      chgl = chgl_;
      mole=mol_;
      fxle = parent;// maybe I need it later?
      neu=neu_;
      imin=imin_;
      imax=imax_;
    }
    ~GenMapWorker(){
      //printf("destruct GenMapWorker\n");
    }

    void run() {
      /* expensive or blocking operation  */
      gen_surface(neu, imin,imax);
    }
    void gen_surface(bool neu,int imin=0, int imax=4);
    ChGL *chgl;
    FourXle *fxle;//parent
    Molecule *mole;
    bool neu;
    int imin,imax;
};


#endif