VarData.h

/***********************************************************************
VarData: Header file for the VarData class.
 Copyright (C) 2008 by Giovanni Marelli <sabeiro@virgilio.it>


This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
***********************************************************************/
#ifndef VARDATA_H
#define VARDATA_H

#include "Matematica.h"
#include <vector> 
#include <list>

#ifdef USE_BOOST
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/bzip2.hpp>
using namespace  boost::iostreams;
#endif

;
#define CHAIN_IF_TYPE(c,t) ( ((c)&(t))==(t) )
#define VAR_IF_TYPE(c,t) ( ((c)&(t))==(t) )
#define VAR_ADD_TYPE(c,t) ( (c)|=(t) )
#define VAR_REM_TYPE(c,t) ( (c)&=(~t) )
enum ShiftType{
  SHIFT_NO      =0,
  SHIFT_CM      =1,
  SHIFT_NANO    =2,
  SHIFT_CM_NANO =3,
};
enum BackFoldType{
  BF_NO       =0,
  BF_PART     =1,
  BF_CHAIN    =2,
  BF_NANO     =3,
  BF_SKIP     =4,
  BF_TILT     =5,
};
enum BeadType{
  BEAD_PHOB   =0,
  BEAD_PHIL   =1,
  BEAD_OIL    =2,
  BEAD_NANO   ,
  BEAD_EVERY  ,
};
enum ChainType{
  CHAIN_EVERY   =0x00000,
  CHAIN_POLY    =0x00001,
  CHAIN_HOMO    =0x00002,
  CHAIN_ADDED   =0x00006,
  CHAIN_UP      =0x00005,
  CHAIN_DOWN    =0x00009,
  CHAIN_INNER   =0x00011,
  CHAIN_OUTER   =0x00021,
  CHAIN_STRETCH =0x00200,
  CHAIN_FLABBY  =0x00400,
  CHAIN_TILTED  =0x00800,
};
enum VarFileFormat{
  VAR_SYS_XVT        =  0x0001,
  VAR_SYS_TXVL       =  0x0002,
  VAR_SYS_XYZ        =  0x0004,
  VAR_SYS_XYZT       =  0x0008,
};
enum VarSysPropr{
  VAR_EDGE           =  0x0001,
  VAR_PART_ALLOCATED =  0x0002,
  VAR_CH_ALLOCATED   =  0x0004,
  VAR_MEMBRANE       =  0x0008,
  VAR_OPEN_TRUST     =  0x0010,
  VAR_LINK_ALLOCATED =  0x0020,
  VAR_CHAIN_DEF      =  0x0040,
  VAR_NBLOCK_ALL     =  0x0080,
};
enum VarCreate{
  VAR_ABSORBING  =  0x0001,
  VAR_REFLECTING =  0x0002,
  VAR_ADDED      =  0x0004,
  VAR_DISTRIBUTED=  0x0008,
  VAR_SOLVENT    =  0x0010,
  VAR_PLANAR     =  0x0020,
  VAR_PLANAR_PE  =  0x0040,
  VAR_TUBE       =  0x0080,
  VAR_VESICLE    =  0x0100,
  VAR_COATING    =  0x0200,
  VAR_TWOTAILS   =  0x0400,
  VAR_OBSTACLE   =  0x0800,
  VAR_OPPOSED    =  0x0f00,
};
enum BlockArch{
  ARCH_POINTS   = 0x0001,
  ARCH_LINES    = 0x0002,
  ARCH_TWOTAILS = 0x0004,
  ARCH_CLUSTER  = 0x0008,
};
enum NanoShape{
  SHAPE_NONE    = 0x00001,
  SHAPE_SPH     = 0x00002,
  SHAPE_HEI     = 0x00004,
  SHAPE_CYL     = 0x00008,
  SHAPE_TILT    = 0x00010,
  SHAPE_PILL    = 0x00020,
  SHAPE_WALL    = 0x00040,
  SHAPE_PORE    = 0x00080,
  SHAPE_SMOOTH  = 0x00100,
  SHAPE_CLUSTER = 0x00200,
  SHAPE_EXT     = 0x00400,
  SHAPE_JANUS   = 0x00800,
  SHAPE_CLINKS  = 0x01000,
  SHAPE_STALK   = 0x02000,
  SHAPE_TIP     = 0x04000,
  SHAPE_TORUS   = 0x08000,
  SHAPE_HARM    = 0x10000,
  SHAPE_UMBR    = 0x20000,
  SHAPE_BOUND   = 0x40000,
};
typedef struct{
  double Pos[3];
  double Bkf[3];
  double Vel[4];
  int Idx;
  int CId;
  int Typ;
}PART;
typedef struct{
  int *Link;
  int NLink;
}LINKS;
typedef struct{
  double Pos[4];
  double Dir[3];
  double Vel[3];
  //  Angle with the normal;
  double Angle;
  int Type;
  int NPCh;
  int InitBead;
  int EndBead;
}CHAIN;
typedef struct{
  char Name[60];
  int NPCh;
  int NPart;
  int NChain;
  int InitIdx;
  int EndIdx;
  int Arch;
  int nNano;
  int Asym;
} BLOCK;
typedef struct{
  char Name[60];
  double Pos[3];
  double Vel[3];
  double Size[3];
  int Topology;
  int NPart;
  int NChain;
  int NPCh;
  int InitIdx;
  int EndIdx;
} SOFT;
typedef struct{
  double Time;
  double Temp;
  double Beta;
  double Energy[3];
  double Edge[4];
  double InvEdge[4];
  double Cm[3];
  double Vel[4];
  double Pre[3];
  double vBB;
  double SurfTens;
  double chiN;
  double kappaN;
  double rho;
  double kappaBend;
  double kappaSpring;
  double SpringRest;
  double ReOverCutOff;
  double Deltat;
  double WFuncStraight2;
  double WFuncStraight3;
  unsigned long Step;
  int NPart;
  int NChain;
  int NPCh;
  int NAllocP;
  int NAllocC;
  int NType;
  int NLink;
  int NBlock;
  //  Number of nanoparticles
  int NNano;
}GENERAL;
class Properties{
 public:
  Properties(){RePhob=0.;RePhil=0.;VolPhob=0.;VolPhil=0.;FactPhob=0.;FactPhil=0.;GyrPhob=0.;GyrPhil=0.;ChDiff=0.;};
  void Print(){printf("RePhob %lf RePhil %lf\n",RePhob,RePhil);
  printf("RadGyrPhob %lf RadGyrPhil %lf\n",GyrPhob,GyrPhil);
  printf("StructFactPhob %lf StructFactPhil %lf\n",FactPhob,FactPhil);
  printf("VolPhob %lf VolPhil %lf\n",VolPhob,VolPhil);
  printf("ChainDiff %lf\n",ChDiff);};
  double RePhob;//End to end distance
  double RePhil;
  double VolPhob;//volume fraction
  double VolPhil;
  double FactPhob;//structure factor
  double FactPhil;
  double GyrPhob;//radius of gyration
  double GyrPhil;
  double ChDiff;//lateral diffusion of the chains
  Properties operator+(const Properties&) const;
  Properties operator*(const double&) const;
};
class MatInt{
 private:
  double *IntMatr2;
  double *IntMatr3;
  int NType;
  int NOrd;
 public:
  MatInt(int NType,int NOrd);
  void FillEntries(double *Matr,int Ord);
  int IntType(int t1,int t2);
  int IntType(int t1,int t2,int t3);
  double Coeff(int t1,int t2);
  double Coeff(int t1,int t2,int t3);
  void SetCoeff(double Co,int t1,int t2);
  void SetCoeff(double Co,int t1,int t2,int t3);
  void Print();
  int pNType(){return NType;};
  void Rescale(double SFactor,int Order);
};
typedef struct{
  char ArchFile[60];
  double Pos[3];
  double Bkf[3];
  double Vel[3];
  double Force[3];
  double Axis[3];
  double AMom[3];
  double AMomTemp[3];
  double AVel[3];
  double Mass;
  double Rad;
  double Hamaker;
  double Height;
  double Gamma;
  double Zeta;
  double Viscosity;
  double OffSet;
  double CutOff;
  double Coating;
  double BaseLine;
  double DistThr;
  double ForThr;
  double PotThr;
  double Area;
  int Shape;
  int NLink;
  int NHeight;
  int NCircle;
  int nBlock;
}NANO;
typedef struct {
  double x[3];
} XYZ;
typedef struct {
  XYZ p[8];
  XYZ n[8];
  double val[8];
} GRIDCELL;
typedef struct {
  XYZ p[3];
  XYZ c;
  int v[3];
  XYZ n[3];
} VAR_TRIANGLE;
typedef struct {
  XYZ p[2];
  XYZ c;
  int v[2];
  XYZ n[2];
} VAR_LINE;
//-----------------Class-VarData----------------------------------
class VarData{
 private:
  GENERAL *Gen;
 public:
  Matematica *Mat;
  MatInt *MInt;
  VarData();
  ~VarData();
  //-------------------------VarData.cpp-------------------
  void VarMessage(const char * s, ...);
  bool Open(char *InFile,int BF);
  bool OpenRisk(char *InFile,int BF);
  bool OpenTrust(char *InFile,int BF);
  void AllocPart();
  void AllocChain();
  Properties SysProperties();
  void SysInfo(char *cSystem);
  void SysDef(char *cSystem);
  char *SysState();
  void SetCoeff();
  void SetCoeff(double *v2,double *v3);
  double TwoPartDist(int p1,int p2,double *RelDist);
  double TwoPartDist(double *Pos,int p2,double *RelDist);
  double TwoPartDist2(int p1,int p2,double *RelDist);
  double TwoPartDist2(double *Pos,int p2,double *RelDist);
  int TwoPartDist(int p1,int p2,double *RelDist,double CutOff);
  int TwoPartDist(double *Pos,int p2,double *RelDist,double CutOff);
  //----------------------VarDataWrite.cpp-----------------
  bool Write(char *OutFile);
  bool WriteTxvl(char *OutFile);
  bool WriteXvt(char *OutFile);
  bool WriteXyz(char *OutFile);                  
  void HeaderInteraction(FILE *FileToWrite);
  void StringNano(char *NString,int n);
  int HeaderNano(FILE *FileToWrite);
  int HeaderSoft(char *Line);
  void WriteLinkedSurf(FILE *FWrite,double *Plot,int NSample,int NType,double *Bound,int *PId);
  void WriteSurf(FILE *F2Write,double **Plot,int NSample,int OffSet);
  void ShapeId(int iShape,char *Shape);
#ifdef USE_BOOST
  void OpenComprWrite(char *FName,filtering_ostream ZipFile);
  void OpenComprRead(char *FName,filtering_istream ZipFile);
  void CloseCompr(filtering_istream ZipFile);
  void CloseCompr(filtering_ostream ZipFile);
#endif
  //------------------------VarDataBackFold.cpp------------
  bool BackFold(int How);
  int BfDefChain();
  int BfEdge();
  void DistFromNp();
  void ShiftRef(int BackFold);
  int StalkPos(double *OldPos);
  void BfPep();
  void BackBone(double *Line,int NBin);
  void StalkLineProf(double *Line,int NBin);
  void StalkPos2(double *OldPos,double *CmStalk);
  void StalkPos3(double *OldPos,double *CmStalk);
  int StalkPos4(double *OldPos,double *CmStalk);
  double NormalWeight(VAR_TRIANGLE *Triang,double *Weight,int NGrid,int NTri);
  void ConnectLineChain(VAR_LINE *Triang,int NGrid,int NTri);
  void ConnectLineChain2(VAR_LINE *Triang,int NGrid,int NTri);
  void ConnectLineChain3(VAR_LINE *Triang,int NGrid,int NTri);
  double PorePos();
  //-----------------VarDataString.cpp------------------
  int Fetch(char *str,char *mask,char *fmt, ... );
  int BraketPos(char *str,char *mask,int *sPos,int *sLen);
  int Fetch(char *str,char *mask,int NArg,double *Val);
  bool ReadString(const char *String,char *cLine,double *Value);
  bool ReadString(const char *String,double *Value,char *line);
  bool ReadString(const char *String,char *cLine,int *Value);
  int  ReadVal(char *pLine,double *Value);
  //-----------------VarDataRead.cpp-------------------
  int ReadLineXvt(char *cLine,double *Pos,int *Type);
  bool ReadConf(char *InFile);
  void ReadHeader(FILE *FileToRead);
  void ReadHeaderTxvl(FILE *FileToRead);
  void ReadHeaderXvt(FILE *FileToRead);
  int ReadPart(FILE *FileToRead);
  int ReadPartTxvl(FILE *FileToRead);
  int ReadPartXvt(FILE *FileToRead);
  int ReadPartXyz(FILE *FileToRead);
  int ReadPartXyzt(FILE *FileToRead);
  int ReadPassThru(FILE *FileToRead);
  int ReadSoft(FILE *ConfFile);
  void ReadNano(FILE *ConfFile,int NCircle,int NHeight);
  int NanoString(char *cLine,int n);
  void SubNanoHeader(char *cFile);
  int ShapeId(char *Shape);
  //--------------------VarDataCreate.cpp------------------------
  int DefSoft(char *nome2,char *ConfF);
  int TrialSys();
  bool CreateSoft(int *arch,double Thickness,int s);
  void CreateTube(int *arch,double Thickness,int s);
  void CreatePlanar(int *arch,double Thickness,int s);
  void CreateVesicle(int *arch,double Thickness,int s);
  void CreateCoating(int *arch,double Thickness,int s);
  void CreateObstacle(int *arch,double Thickness,int s);
  int CheckNano(double *Pos,int s);
  void AddProtein(int NCircle,int NHeight,int nNano,char *filename);
  void CreateProtein(int nNano,int nStart);
  void AddStuffing(char *filename,int nStuffing,int nNano);
  void AddSolvent(char *filename,int nWater);
  void AddChains(char *filename,double Thickness);
  void AddCholesterol(char *filename,double Thickness,int s);
  void DefBlock(int *NChStep,int How);
  void DefRest(int *arch,int s);
  int PutPart(int j,int p,int HalfLim,double sigma);
  void FindNeighbours(char *FileName);
  //---------------------VarDataEl.cpp----------------------------
  void SwapChain(int c1,int c2,int b);
  void SwapChain(int c1,int c2);
  void SwapPart(int p1,int p2);
  void ChangeNChain(int NChain,int b);
  bool ShiftSys(int How);
  void SampleSurface(double *Plot,int NSample,int Type);
  MOMENTI SampleSurfacePart(double *Plot,int NSample,int Type);
  MOMENTI SampleSurface(Matrice *Plot,int NSample,int Type);
  MOMENTI SampleSurfaceMem(int NSample);
  void LoadDensFile(double **Plot,int NBin);
  int SpatialDerivative(Matrice *Surface,Matrice *Resp,SPLINE Weight,int NSample);
  void ShiftBlock(Vettore *Shift,int b);
  void RotateBlock(Vettore *Axis,Vettore *Origin,int b);
  void MirrorBlock(Vettore *Px1,Vettore *Px2,Vettore *Px3,int b);
  void Transform(int block);
  void Point2Shape(int iShape);
  typedef double(VarData::*NANO_DIST)(double *Pos,int n);
  NANO_DIST Nano_Dist;
  double NanoDist2(double *Pos,int n){return (*this.*Nano_Dist)(Pos,n);}
  double NanoDist2(double x,double y,double z,int n);
  double FieldNo(double *Pos,int n);
  double FieldSphere(double *Pos,int n);
  double FieldElips(double *Pos,int n);
  double FieldParab(double *Pos,int n);
  double FieldCyl(double *Pos,int n);
  double FieldTransMem(double *Pos,int n);
  double FieldJanus(double *Pos,int n);
  double FieldTorus(double *Pos,int n);
  double FieldTilt(double *Pos,int n);
  double FieldBound(double *Pos,int n);
  double FieldTiltWall(double *Pos,int n);
  //---------------------VarDataExp.cpp----------------------------
  int PairCorrelation(double *Point,int NSample,int How,int Type);
  int PairCorrelationRound(double **Point,int NSample,int Type);
  int PairCorrelationSquare(double **Point,int NSample,int Type);
  int PairCorrelationPep(double **Point,int NSample,int Type);
  int Scattering2d(double **Point,int NSample,int Type);
  int Scattering2D(double **Point,int NSample,int Type);
  void Spettro2d(double *Points,int NSample,int Type);
  void Spettro2d(double *Plot,int NSample);
  //---------------------VarDataContour.cpp------------------------
  int DensityProfile(int coord,int NSample,int NType,double *dDensity);
  int Core(double ***Plot,int NSample,double Border[3][2]);
  int RadDistr(int NSample,double *Plot,double Border[2],int How);
  int Worm(int Partition,int NSample,double *Border,double *dPoint);
  void VolumeCircSlab(double *VolContr,int NSample);
  void Stalk(int NSample,int NLevel,double **Plot,double Threshold);
  //---------------------VarDataPos.cpp----------------------------
  int Arrange(int **Triangle,int Vertex);
  int Folding();
  int OrderPos();
  int CalcnPos(double *Pos);
  int Neighbour(double *Pos);
  int NChainPSquare(double *Plot);
  int LateralFluctuation(double *Plot,int LatValue);
  int Voronoi();
  int PosVectInt(double *Pos);
  //------------------VarInterp.cpp----------------------------
  int InterParab(PART *PmIn,PART *PmOut,int NIn,int nOut);
  int InterParab2(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterCubica(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterForth(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterSpline3(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterSpline4(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterBSpline(PART *PmIn,PART *PmOut,int NIn,int NOut);
  int InterBSpline2D(double **PlIn,double **PmOut,int NIn,int NOut);
  int InterBSpline2D(double *PlIn,double *PmOut,int NIn,int NOut);
  int InterBSpline1D(double *PlIn,double *PmOut,int NIn,int NOut);
  int InterPoly(PART *PmIn,PART *PmOut,int NIn,int nOut);
  int InterDerMatrix(PART *Pm, int NMass,SPLINE Weight,double Offset);
  void SmoothGrid(int NSample,char *FWrite);
  void SmoothGrid(int NSample);
  void ConvoluteMatrix(double *Plot,int NGrid,Matrice *Mask,int NDim);
  void ConvoluteMatrix1(double *Plot,int NGrid,Matrice *Mask);
  void ConvoluteMatrix2(double *Plot,int NGrid,Matrice *Mask);
  void ConvoluteMatrix3(double *Plot,int NGrid,Matrice *Mask);
  //-------------------VarDataComm.cpp----------------------
  int SetNPart(int NewNPart);
  int SetNChain(int NewNCh);
  int SetNLink(int NewNCh);
  void SetNPCh(int NewNCh);
  void SetNType(int NewNType);
  int AllocLinks(int NewNCh);
  int SetNBlock(int Val);
  int SetNNano(int Val);
  void Copy(PART *P1,PART *P2,int NPartOld);
  void Copy(CHAIN *C1,CHAIN *C2,int NChainOld);
  //-------------------VarDataMarchCubes.cpp----------------
  VAR_TRIANGLE *MarchingCubes(double *Plot,int NSample,double IsoLevel,int *NTri);
  VAR_LINE *MarchingSquares(double *Plot,int NSample,double IsoLevel,int *NTri);
  //-------------------VarDataCGAL.cpp----------------------
  //#ifdef __CGAL_h__
  void AreaDistr(double *Distr,double *RadDistr,int NSample);
  //#endif //CGAL
  double pTime(){return Gen->Time;};
  double pDeltat(){return Gen->Deltat;};
  double pTemp(){return Gen->Temp;};
  double pBeta(){return Gen->Beta;};
  double pEnergy(int d){return Gen->Energy[d];};
  double pEdge(int d){return Gen->Edge[d];};
  double pInvEdge(int d){return Gen->InvEdge[d];};
  double pVol(){return Gen->Edge[0]*Gen->Edge[1]*Gen->Edge[2];};
  double pCm(int d){return Gen->Cm[d];};
  double pVelMax(int d){return Gen->Vel[d];};
  double pchiN(){return Gen->chiN;};
  double pkappaN(){return Gen->kappaN;};
  double pkBen(){return Gen->kappaBend;};
  double pkSpr(){return Gen->kappaSpring;};
  void SetkBen(double Val){Gen->kappaBend = Val;};
  void SetkSpr(double Val){Gen->kappaSpring = Val;};
  void SetSprRest(double Val){Gen->SpringRest = Val;};
  double pSprRest(){return Gen->SpringRest;};
  double prho(){return Gen->rho;};
  double pReOverCutOff(){return Gen->ReOverCutOff;};
  double pWei2Par(){return Gen->WFuncStraight2;};
  double pWei3Par(){return Gen->WFuncStraight3;};
  int pStep();
  int pNPart();
  int pNChain();
  int pNChain(int b);
  int pNPCh();
  int pNPCh(int c);
  int pNType();
  int pNLink();
  int pNNano();
  int pNBlock();
  int pNAllocP();
  int pNAllocC();
  void SetEdge(double Val,int d){Gen->Edge[d] = Val;Gen->InvEdge[d] = 1./Val;};
  void SetCNorm(int d){
    if(d<0 || d > 2) return;
    CNorm = d;
    CLat1 = (d+1)%3;
    CLat2 = (d+2)%3;
  };
  void SetScaleF(double *Scale){
    for(int d=0;d<3;d++){
      ScaleF[d] = Scale[d];
    }
  }
  void SetShiftPos(double *RefPos){
    for(int d=0;d<3;d++){
      ShiftPos[d] = RefPos[d];
    }
  }
  void SetIfNormalize(int If){
    IfNormalize = If;
  }
  void SetDeltat(double Val){Gen->Deltat = Val;};
  void SetStep(int Val){Gen->Step = Val;};
  void SetTemp(double Val){Gen->Temp = Val;Gen->Beta = 1./Val;};
  void SetTime(double Val){Gen->Time = Val;};
  void IncrStep(){Gen->Step++;};
  double pPos(int p,int d);
  double pChPos(int p,int d);
  void pPos(int p,double *Pos);
  double *pPos(int p);
  double pPosNoBkf(int p,int d);
  double pVel(int p,int d);
  void SetPos(int p, double *Pos);
  void SetPos(int p,int d,double Pos);
  void SetVel(int p, double *Vel);
  void SetType(int p,int t);
  int pType(int p);
  int pChain(int p);
  double pNanoPos(int n,int d);
  void SetBkf(int p);
  void SetNanoBkf(int n);
  void pPos(double *Pos);
  char cWhat2Draw[STRSIZE];
  NANO *Nano;
  PART *Pm;
  LINKS *Ln;
  CHAIN *Ch;
  SOFT *Soft;
  BLOCK *Block;
  double *PlotMem;
  double ShiftPos[3];
  double ScaleF[3];
  int NSoft;
  int NPartNearSphere;
  int NAddChain;
  int NAddChol;
  int NSolvent;
  int NStuffing;
  int CNorm;
  int CLat1;
  int CLat2;
  int NChType;
  int NPType;
  int NEdge;
  int SysType;
  int SysFormat;
  int SysCreate;
  int IfNormalize;
  int IfPlotMem;
};
#endif //VARIABILI_H