SM_ZeemanOperator.h

#ifndef SM_ZeemanOperator_H
#define SM_ZeemanOperator_H
 
#include "SM_Operator.h"
 
#include "SM_Field.h"
 
class SM_ZeemanOperator : public SM_SliceOperator  {
    
    //attributes
private :
    
    //super class type
    typedef SM_SliceOperator SuperClass;
    
    //association
 
    //nTimeSteps x nParticles x dim
    //Hext(t,i,k)=mHext[k+dim*(i+P*t)]
    
    //nTimeSteps=1 for permanent Hext
    //nParticles=1 for constant field
    tIndex mTimeStepsNumber;
    tIndex mParticlesNumber;
 
    //external field
    std::valarray<tReal> mHext;
    
   
protected:
    // CONSTRUCTORS 
    SM_ZeemanOperator(void) {
        this->setName("Zeeman");
        mTimeStepsNumber=1;
        mParticlesNumber=1;
    }
    
    // DESTRUCTORS 
    virtual ~SM_ZeemanOperator(void) {
    }
    
public : 
 
 
    //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SuperClass::getContentsMemorySize();
        mem+=mHext.size()*sizeof(tReal);
        return mem;
    }
     //Instance building
    //=================
 
   
    inline static CORE_UniquePointer<SM_ZeemanOperator> New() {
        return CORE_UniquePointer<SM_ZeemanOperator>(new SM_ZeemanOperator(),
                                                     SM_ZeemanOperator::Delete());
    }
 
    
    virtual CORE_UniquePointer<SM_SliceOperator> NewInstance() const override {
        return New();
    }
    
protected:
    inline void setParticlesNumber(const tIndex& p) {
        mParticlesNumber=p;
    }
    inline void setTimeStepsNumber(const tIndex& n) {
        mTimeStepsNumber=n;
    }
    
public:
    virtual void copy(const SM_SliceOperator& op) {
        SuperClass::copy(op);
        const SM_ZeemanOperator *zop=dynamic_cast<const SM_ZeemanOperator*>(&op);
        if (zop!=null) {
            mParticlesNumber=zop->getParticlesNumber();
            mTimeStepsNumber=zop->getTimeStepsNumber();
            mHext=zop->getHext();
        }
    }
 
    inline const tIndex& getParticlesNumber() const {return mParticlesNumber;}
    
    inline const tIndex& getTimeStepsNumber() const {return mTimeStepsNumber;}
   
    
    inline const tReal* getHextAtTimeIndex(const tIndex& timeIndex) const {
        return &mHext[(timeIndex%mTimeStepsNumber)*mParticlesNumber*SM_Constants::DIM];
    }
    inline  tReal* getHextAtTimeIndex(const tIndex& timeIndex)  {
        return &mHext[(timeIndex%mTimeStepsNumber)*mParticlesNumber*SM_Constants::DIM];
    }
    inline const tReal* getHextValues() const {
        return &mHext[0];
    }
    inline tReal* getHexValues()  {
        return &mHext[0];
    }
protected:
    inline const std::valarray<tReal>& getHext()  const {
        return mHext;
    }
    inline std::valarray<tReal>& getHext()  {
        return mHext;
    }
public:
    inline void setHext(std::initializer_list<tReal>&& values) {
        mHext=values;
    }
   
    inline void setHext(const std::vector<tReal>& values) {
        mHext.resize(values.size());
        tReal *Vi=&mHext[0];
        for_each(values.begin(),values.end(),
                 [&](const auto &v) {
                     (*Vi)=v;
                     Vi++;
                 });
    }
    
    inline  void setHext(const std::array<tReal,SM_Constants::DIM>& values) {
      
        std::valarray<tReal>& Hext=getHext();
        Hext.resize(SM_Constants::DIM);
        memcpy(&Hext[0],values.data(),SM_Constants::DIM*sizeof(tReal));
    }
 
    virtual void adimensionize(const tReal& cH) override {
        //begin iterator on Hext
        tReal *iHext=&mHext[0];
 
        //end iterator on Hext
        const tReal *eHext=iHext;
        eHext+=mHext.size();
 
        //dim factor
        tReal invCH=1.;
        invCH/=cH;
        
        
        while (iHext!=eHext) {
            (*iHext)*=invCH;
            iHext++;
        }
    }
    virtual void discretize(const SM_Network& network, const SM_Material& material) override {
        setParticlesNumber(network.getParticlesNumber());
        setTimeStepsNumber(1);
    }
 
 
    //magnetics fields
    //================
   
    
   
    
    virtual void computeMagneticFieldSlice(const tIndex& timeIndex,
                                           const SM_Network& network,
                                           const SM_Material& material,
                                           const tIndex& startIndex,
                                           const tIndex& endIndex,
                                           const tReal *S,
                                           const tBoolean& alpha,
                                           const tIndex& nH,
                                           tReal *H) const;
    //spin energy
    //===========
 
    virtual tReal computeSpinEnergy(const tIndex& i,
                                    const tIndex& timeIndex,
                                    const SM_Network& network,
                                    const SM_Material& material,
                                    const SM_RealField& S) const override {
        //begin iterator on direction of magnetic moment of spin i
        const tReal *Si=&S[0];Si+=i*SM_Constants::DIM;
        //end iterator on direction of magnetic moment of spin i
        const tReal *eSi=Si;eSi+=SM_Constants::DIM;
        //energi of the spin i
        tReal Ei;
        //compute spin energy
        ComputeSpinEnergy(material.getAtomicSpinMoment()/material.getAdimensionizedDerivativeEnergyFactor(),//Mu_s=\tilde mu_s. Mu_B or only \tilde mu_s
                          &getHextAtTimeIndex(timeIndex)[i*SM_Constants::DIM],//Hext
                          Si,eSi,//S at i
                          Ei);
        return Ei;
    }
 
    inline static void ComputeSpinEnergy(const tReal& mu_s,
                                         const tReal *Hext,
                                         const tReal *iS,const tReal *eSi,
                                         tReal& E) {
        ComputeUnscaledSpinEnergy(Hext,iS,eSi,E);
        E*=mu_s;
        E*=-1;
    }
  
protected:
    inline static void ComputeUnscaledSpinEnergy(const tReal *Hext,
                                                 const tReal *iS,const tReal *eS,
                                                 tReal& E) {
 
        E=0;
        while (iS!=eS) {//loop on particles in [start,end[
            
            //E= \sum_{i=0}^{i=P-1}  < H_{ext}(P_i),S_i> 
            E+=(*Hext)*(*iS);
            
            //iterator at next values
            Hext++;
            iS++;
        }//end loop on particles
    }
    
    
    //total energy
    //============
public:
  
    virtual  tReal computeEnergySlice(const tIndex& timeIndex,
                                      const SM_Network& network,
                                      const SM_Material& material,
                                      const tIndex& startIndex,
                                      const tIndex& endIndex,
                                      const tReal *S) const;
 
 
 
 
    //string representation
    //=====================
    
    virtual tString toString() const override {
        tString ret=SuperClass::toString();
        ret+=" Hext:"+functions_array::toString(getHext())+"\n";
        return ret;
    }
  
};
 
 
#endif