ScatteringData.h

/*  ObjCryst++ Object-Oriented Crystallographic Library
    (c) 2000-2002 Vincent Favre-Nicolin vincefn@users.sourceforge.net
        2000-2001 University of Geneva (Switzerland)

    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 _OBJCRYST_SCATTERINGDATA_H_
#define _OBJCRYST_SCATTERINGDATA_H_

#include "ObjCryst/CrystVector/CrystVector.h"

#include "ObjCryst/ObjCryst/General.h"

#include "ObjCryst/ObjCryst/SpaceGroup.h"
#include "ObjCryst/ObjCryst/ScatteringPower.h"
#include "ObjCryst/ObjCryst/Scatterer.h"
#include "ObjCryst/ObjCryst/Crystal.h"

//#include <stdlib.h>
#include <string>
//#include <iomanip>
//#include <cmath>
//#include <typeinfo>
//#include <fstream>
//#include <ctime>

namespace ObjCryst
{
#ifndef HAVE_SSE_MATHFUN
//initialize tabulated values of cosine
void InitLibCrystTabulCosine();
void DeleteLibCrystTabulCosine();
void InitLibCrystTabulExp();
void DeleteLibCrystTabulExp();
#endif
extern const RefParType *gpRefParTypeScattData;
extern const RefParType *gpRefParTypeScattDataScale;
extern const RefParType *gpRefParTypeScattDataProfile;
extern const RefParType *gpRefParTypeScattDataProfileType;
extern const RefParType *gpRefParTypeScattDataProfileWidth;
extern const RefParType *gpRefParTypeScattDataProfileAsym;
extern const RefParType *gpRefParTypeScattDataCorr;
extern const RefParType *gpRefParTypeScattDataCorrInt;
extern const RefParType *gpRefParTypeScattDataCorrIntPO_Direction;
extern const RefParType *gpRefParTypeScattDataCorrIntPO_Fraction;
extern const RefParType *gpRefParTypeScattDataCorrIntPO_Amplitude;
extern const RefParType *gpRefParTypeScattDataCorrInt_Ellipsoid;
extern const RefParType *gpRefParTypeScattDataCorrIntAbsorp;
extern const RefParType *gpRefParTypeScattDataCorrIntPolar;
extern const RefParType *gpRefParTypeScattDataCorrIntExtinc;
extern const RefParType *gpRefParTypeScattDataCorrPos;
extern const RefParType *gpRefParTypeScattDataBackground;

extern const RefParType *gpRefParTypeRadiation;
extern const RefParType *gpRefParTypeRadiationWavelength;

class NiftyStaticGlobalObjectsInitializer_ScatteringData
{
   public:
      NiftyStaticGlobalObjectsInitializer_ScatteringData()
      {
         if (mCount++ == 0)
         {
            #ifndef HAVE_SSE_MATHFUN
            InitLibCrystTabulCosine();
            InitLibCrystTabulExp();
            #endif
            gpRefParTypeScattData= new RefParType(gpRefParTypeObjCryst,"Scattering Data");
            gpRefParTypeScattDataScale= new RefParType(gpRefParTypeObjCryst,"Scale Factor");
            gpRefParTypeScattDataProfile= new RefParType(gpRefParTypeScattData,"Profile");
            gpRefParTypeScattDataProfileType= new RefParType(gpRefParTypeScattDataProfile,"Type");
            gpRefParTypeScattDataProfileWidth= new RefParType(gpRefParTypeScattDataProfile,"Width");
            gpRefParTypeScattDataProfileAsym= new RefParType(gpRefParTypeScattDataProfile,"Asymmetry");
            gpRefParTypeScattDataCorr= new RefParType(gpRefParTypeScattData,"Correction");
            gpRefParTypeScattDataCorrInt= new RefParType(gpRefParTypeScattDataCorr,"Intensities");
            gpRefParTypeScattDataCorrIntPO_Direction= new RefParType(gpRefParTypeScattDataCorrIntPO_Direction,"Preferred orientation direction");
            gpRefParTypeScattDataCorrIntPO_Fraction= new RefParType(gpRefParTypeScattDataCorrIntPO_Fraction,"Preferred orientation fraction");
            gpRefParTypeScattDataCorrIntPO_Amplitude= new RefParType(gpRefParTypeScattDataCorrIntPO_Amplitude,"Preferred orientation amplitude");
            gpRefParTypeScattDataCorrInt_Ellipsoid= new RefParType(gpRefParTypeScattDataCorrInt_Ellipsoid,"Preferred orientation ellipsoid");
            gpRefParTypeScattDataCorrIntAbsorp= new RefParType(gpRefParTypeScattDataCorrInt,"Absorption");
            gpRefParTypeScattDataCorrIntPolar= new RefParType(gpRefParTypeScattDataCorrInt,"Polarization");
            gpRefParTypeScattDataCorrIntExtinc= new RefParType(gpRefParTypeScattDataCorrInt,"Extinction");
            gpRefParTypeScattDataCorrPos= new RefParType(gpRefParTypeScattDataCorr,"Reflections Positions");
            gpRefParTypeScattDataBackground= new RefParType(gpRefParTypeScattData,"Background");
            gpRefParTypeRadiation= new RefParType(gpRefParTypeObjCryst,"Radiation");
            gpRefParTypeRadiationWavelength= new RefParType(gpRefParTypeRadiation,"Wavelength");
         }
      }
      ~NiftyStaticGlobalObjectsInitializer_ScatteringData()
      {
         if (--mCount == 0)
         {
            #ifndef HAVE_SSE_MATHFUN
            DeleteLibCrystTabulCosine();
            DeleteLibCrystTabulExp();
            #endif
            delete gpRefParTypeScattData;
            delete gpRefParTypeScattDataScale;
            delete gpRefParTypeScattDataProfile;
            delete gpRefParTypeScattDataProfileType;
            delete gpRefParTypeScattDataProfileWidth;
            delete gpRefParTypeScattDataProfileAsym;
            delete gpRefParTypeScattDataCorr;
            delete gpRefParTypeScattDataCorrInt;
            delete gpRefParTypeScattDataCorrIntAbsorp;
            delete gpRefParTypeScattDataCorrIntPolar;
            delete gpRefParTypeScattDataCorrIntExtinc;
            delete gpRefParTypeScattDataCorrPos;
            delete gpRefParTypeScattDataBackground;
            delete gpRefParTypeRadiation;
            delete gpRefParTypeRadiationWavelength;
            gpRefParTypeScattData=0;
            gpRefParTypeScattDataScale=0;
            gpRefParTypeScattDataProfile=0;
            gpRefParTypeScattDataProfileType=0;
            gpRefParTypeScattDataProfileWidth=0;
            gpRefParTypeScattDataProfileAsym=0;
            gpRefParTypeScattDataCorr=0;
            gpRefParTypeScattDataCorrInt=0;
            gpRefParTypeScattDataCorrIntAbsorp=0;
            gpRefParTypeScattDataCorrIntPolar=0;
            gpRefParTypeScattDataCorrIntExtinc=0;
            gpRefParTypeScattDataCorrPos=0;
            gpRefParTypeScattDataBackground=0;
            gpRefParTypeRadiation=0;
            gpRefParTypeRadiationWavelength=0;
         }
      }
   private:
      static long mCount;
};
static NiftyStaticGlobalObjectsInitializer_ScatteringData NiftyStaticGlobalObjectsInitializer_ScatteringData_counter;
//######################################################################
//######################################################################
class Radiation: public RefinableObj
{
   public:
      Radiation();
      Radiation(const RadiationType rad,const REAL wavelength);
      Radiation(const string &XRayTubeElementName,const REAL alpha2Alpha2ratio=0.5);
      Radiation(const Radiation&);
      ~Radiation();
      virtual const string& GetClassName() const;

      void operator=(const Radiation&);

      RadiationType GetRadiationType()const;
      void SetRadiationType(const RadiationType);
      void SetWavelengthType(const WavelengthType &type);
      WavelengthType GetWavelengthType()const;
      const CrystVector_REAL& GetWavelength()const;
      void SetWavelength(const REAL );
      void SetWavelength(const string &XRayTubeElementName,const REAL alpha2Alpha2ratio=0.5);

      REAL GetXRayTubeDeltaLambda()const;
      REAL GetXRayTubeAlpha2Alpha1Ratio()const;

      const RefinableObjClock& GetClockWavelength()const ;
      const RefinableObjClock& GetClockRadiation()const ;
      virtual void XMLOutput(ostream &os,int indent=0)const;
      virtual void XMLInput(istream &is,const XMLCrystTag &tag);
      //virtual void XMLInputOld(istream &is,const IOCrystTag &tag);
      void Print()const;
      REAL GetLinearPolarRate()const;
      void SetLinearPolarRate(const REAL f);
   private:
      void InitOptions();
      RefObjOpt mRadiationType;
      RefObjOpt mWavelengthType;
      CrystVector_REAL mWavelength;
      string mXRayTubeName;
      REAL mXRayTubeDeltaLambda;
      REAL mXRayTubeAlpha2Alpha1Ratio;
      //Clocks
         RefinableObjClock mClockWavelength;
         RefinableObjClock mClockRadiation;
      REAL mLinearPolarRate;
   #ifdef __WX__CRYST__
   public:
      virtual WXCrystObjBasic* WXCreate(wxWindow*);
      friend class WXRadiation;
   #endif
};

//######################################################################
//######################################################################
class ScatteringData: virtual public RefinableObj
{
   public:
      ScatteringData();
      ScatteringData(const ScatteringData &old);
      ~ScatteringData();
      virtual ScatteringData* CreateCopy()const=0;

      virtual void SetHKL( const CrystVector_REAL &h,
                           const CrystVector_REAL &k,
                           const CrystVector_REAL &l);
      virtual void GenHKLFullSpace2(const REAL maxsithsl,
                                   const bool unique=false);
      virtual void GenHKLFullSpace(const REAL maxTheta,
                                   const bool unique=false);

      RadiationType GetRadiationType()const;
      virtual const Radiation& GetRadiation()const=0;

      virtual void SetCrystal(Crystal &crystal);
      const Crystal& GetCrystal()const ;
      Crystal& GetCrystal() ;
      bool HasCrystal()const;

      long GetNbRefl() const;
      const CrystVector_REAL& GetH() const;
      const CrystVector_REAL& GetK() const;
      const CrystVector_REAL& GetL() const;
      const CrystVector_REAL& GetH2Pi() const;
      const CrystVector_REAL& GetK2Pi() const;
      const CrystVector_REAL& GetL2Pi() const;
      const CrystVector_REAL& GetReflX() const;
      const CrystVector_REAL& GetReflY() const;
      const CrystVector_REAL& GetReflZ() const;

      const CrystVector_REAL& GetSinThetaOverLambda()const;
      const CrystVector_REAL& GetTheta()const;
      const RefinableObjClock& GetClockTheta()const;

      const CrystVector_REAL& GetFhklCalcSq() const;
      std::map<RefinablePar*, CrystVector_REAL> & GetFhklCalcSq_FullDeriv(std::set<RefinablePar *> &vPar);
      const CrystVector_REAL& GetFhklCalcReal() const;
      const CrystVector_REAL& GetFhklCalcImag() const;

      const CrystVector_REAL& GetFhklObsSq() const;

      const map<const ScatteringPower*,CrystVector_REAL> &GetScatteringFactor() const;

      CrystVector_REAL GetWavelength()const;

      void SetIsIgnoringImagScattFact(const bool b);
      bool IsIgnoringImagScattFact() const;
      // Set an option so that only low-amgle reflections (theta < angle)
      // are used. See DiffractionData::mUseOnlyLowAngleData
      //virtual void SetUseOnlyLowAngleData(const bool useOnlyLowAngle,const REAL angle)=0;
      virtual void PrintFhklCalc(ostream &os=cout)const;
      virtual void PrintFhklCalcDetail(ostream &os=cout)const;

      virtual void BeginOptimization(const bool allowApproximations=false,
                                     const bool enableRestraints=false);
      virtual void EndOptimization();
      virtual void SetApproximationFlag(const bool allow);
      virtual void SetMaxSinThetaOvLambda(const REAL max);
      REAL GetMaxSinThetaOvLambda()const;
      virtual long GetNbReflBelowMaxSinThetaOvLambda()const;
      const RefinableObjClock& GetClockNbReflBelowMaxSinThetaOvLambda()const;
   protected:
      virtual void PrepareHKLarrays() ;
      virtual CrystVector_long SortReflectionBySinThetaOverLambda(const REAL maxSTOL=-1.);
      CrystVector_long EliminateExtinctReflections();

      //The following functions are used during the calculation of structure factors,
         virtual void CalcSinThetaLambda()const;
         REAL CalcSinThetaLambda(REAL h, REAL k, REAL l)const;
         virtual const CrystMatrix_REAL& GetBMatrix()const;
         void CalcScattFactor()const;
         void CalcTemperatureFactor()const;
         virtual void CalcResonantScattFactor()const;
         void CalcGlobalTemperatureFactor() const;

      void CalcStructFactor() const;
      void CalcStructFactor_FullDeriv(std::set<RefinablePar *> &vPar);
      void CalcGeomStructFactor() const;
      void CalcGeomStructFactor_FullDeriv(std::set<RefinablePar*> &vPar);
      void CalcLuzzatiFactor()const;
      void CalcStructFactVariance()const;

      long mNbRefl;
      CrystVector_REAL mH, mK, mL ;
      mutable CrystVector_long mIntH, mIntK, mIntL ;
      mutable CrystVector_REAL mH2Pi, mK2Pi, mL2Pi ;
      mutable CrystVector_REAL mX, mY, mZ ;

      CrystVector_int mMultiplicity ;

      CrystVector_int mExpectedIntensityFactor;

      mutable CrystVector_REAL mFhklCalcReal, mFhklCalcImag ;
      mutable std::map<RefinablePar*, CrystVector_REAL> mFhklCalcReal_FullDeriv, mFhklCalcImag_FullDeriv ;
      mutable CrystVector_REAL mFhklCalcSq ;
      mutable std::map<RefinablePar*, CrystVector_REAL> mFhklCalcSq_FullDeriv;

      Crystal *mpCrystal;

      REAL mGlobalBiso;
      mutable CrystVector_REAL  mGlobalTemperatureFactor;

      bool mUseFastLessPreciseFunc;

      //The Following members are only kept to avoid useless re-computation
      //during global refinements. They are used \b only by CalcStructFactor()

         mutable CrystVector_REAL mSinThetaLambda;

         mutable CrystVector_REAL mTheta;

         mutable map<const ScatteringPower*,REAL> mvFprime,mvFsecond;

         mutable map<const ScatteringPower*,CrystVector_REAL> mvTemperatureFactor;

         mutable map<const ScatteringPower*,CrystVector_REAL> mvScatteringFactor;

         mutable map<const ScatteringPower*,CrystVector_REAL> mvRealGeomSF,mvImagGeomSF;
         mutable map<RefinablePar*,map<const ScatteringPower*,CrystVector_REAL> > mvRealGeomSF_FullDeriv,mvImagGeomSF_FullDeriv;

      //Public Clocks
         RefinableObjClock mClockHKL;
         mutable RefinableObjClock mClockStructFactor;
         mutable RefinableObjClock mClockStructFactorSq;
      //Internal Clocks
         mutable RefinableObjClock mClockTheta;
         mutable RefinableObjClock mClockScattFactor;
         mutable RefinableObjClock mClockScattFactorResonant;
         mutable RefinableObjClock mClockGeomStructFact;
         mutable RefinableObjClock mClockThermicFact;

         RefinableObjClock mClockGlobalBiso;
         mutable RefinableObjClock mClockGlobalTemperatureFact;

      bool mIgnoreImagScattFact;

      // Maximum sin(theta)/lambda
         REAL mMaxSinThetaOvLambda;
         mutable long mNbReflUsed;
         mutable RefinableObjClock mClockNbReflUsed;

      // Maximum Likelihood
         mutable map<const ScatteringPower*,CrystVector_REAL> mvLuzzatiFactor;
         mutable CrystVector_REAL mFhklCalcVariance;
         mutable RefinableObjClock mClockLuzzatiFactor;
         mutable RefinableObjClock mClockFhklCalcVariance;
      CrystVector_REAL mFhklObsSq;
      RefinableObjClock mClockFhklObsSq;
   #ifdef __WX__CRYST__
      //to access mMaxSinThetaOvLambda
      friend class WXDiffractionSingleCrystal;
      friend class WXPowderPattern;
   #endif
};

}//namespace ObjCryst
#endif // _OBJCRYST_SCATTERINGDATA_H_