SM_ZeemanOperator.h

Go to the documentation of this file.

#ifndef SM_ZeemanOperator_H
#define SM_ZeemanOperator_H
 
#include "SM_TemplatedOperator.h"
 
#include "SM_Field.h"
 
class SM_ZeemanOperator : public SM_TemplatedOperator<SM_ZeemanOperator>  {
    
    //attributes
private :
    
  
   
    //association
 
    //nTimeSteps x nParticles x dim
    //Bext(t,i,k)=mBext[k+dim*(i+P*t)]
    //nTimeStep=1 for permanent Bext
    //nParticles=1 for constant field
 
    tIndex mZeemanTimeStepsNumber;
    tIndex mZeemanParticlesNumber;
    tDimension mDimension;
    std::valarray<tReal> mBext;
    
   
protected:
    // CONSTRUCTORS 
    SM_ZeemanOperator(void) {
        mZeemanTimeStepsNumber=1;
        mZeemanParticlesNumber=1;
        mDimension=3;
        setName("Zeeman");
    }
    
    // DESTRUCTORS 
    virtual ~SM_ZeemanOperator(void) {
    }
    
public : 
 
    //Instance building
    //=================
    inline static CORE_UniquePointer<SM_ZeemanOperator> New() {
        return CORE_UniquePointer<SM_ZeemanOperator>(new SM_ZeemanOperator(),
                                                     SM_ZeemanOperator::Delete());
    }
 
    virtual CORE_UniquePointer<SM_Operator> NewInstance() const override {
        return New();
    }
 
        //MEMORY
    
    virtual tMemSize getMemorySize() const {
        return sizeof(*this)+getContentsMemorySize();
    }
    virtual tMemSize getContentsMemorySize() const {
        tMemSize mem=SM_TemplatedOperator<SM_ZeemanOperator>::getContentsMemorySize();
        mem+=mBext.size()*sizeof(tReal);
        return mem;
    }
    
protected:
    inline void setZeemanParticlesNumber(const tIndex& p) {
        mZeemanParticlesNumber=p;
    }
    inline void setZeemanTimeStepsNumber(const tIndex& n) {
        mZeemanTimeStepsNumber=n;
    }
    
public:
    virtual void copy(const SM_Operator& op) {
        SM_Operator::copy(op);
        const SM_ZeemanOperator *zop=dynamic_cast<const SM_ZeemanOperator *>(&op);
        if (zop!=null) {
            mZeemanParticlesNumber=zop->getZeemanParticlesNumber();
            mZeemanTimeStepsNumber=zop->getZeemanTimeStepsNumber();
            mDimension=zop->getZeemanDimension();
            mBext=zop->getBextAsArray();
        }
    }
 
    inline const tIndex& getZeemanParticlesNumber() const {return mZeemanParticlesNumber;}
    
    inline const tIndex& getZeemanTimeStepsNumber() const {return mZeemanTimeStepsNumber;}
    inline const tDimension& getZeemanDimension() const {return mDimension;}
    
    virtual void discretize(const SM_System& system) override;
    
    inline const tReal* getBext(const tIndex& t) const {
        return &mBext[(t%mZeemanTimeStepsNumber)*mZeemanParticlesNumber*mDimension];
    }
    inline  tReal* getBext(const tIndex& t)  {
        return &mBext[(t%mZeemanTimeStepsNumber)*mZeemanParticlesNumber*mDimension];
    }
    inline const tReal* getBext() const {
        return &mBext[0];
    }
    inline tReal* getBext()  {
        return &mBext[0];
    }
protected:
    inline const std::valarray<tReal>& getBextAsArray()  const {
        return mBext;
    }
    inline std::valarray<tReal>& getBextAsArray()  {
        return mBext;
    }
public:
    inline void setBext(std::initializer_list<tReal>&& values) {
        mBext=values;
    }
   
    inline void setBext(const std::vector<tReal>& values) {
        mBext.resize(values.size());
        tReal *Vi=&mBext[0];
        for_each(values.begin(),values.end(),
                 [&](const auto &v) {
                     (*Vi)=v;
                     Vi++;
                 });
    }
    
      
    virtual void computeMagneticField(const tIndex& t,
                                      const SM_Network& network,
                                      const tReal *mu,
                                      tReal *B) const override {
        return computeTemplatedMagneticField(t,network,mu,B);
    }
    virtual void computeMagneticFieldAndEnergy(const tIndex& t,
                                               const SM_Network& network,
                                               const tReal *mu,
                                               tReal *B,
                                               tReal& E) const override{
        return computeTemplatedMagneticFieldAndEnergy(t,network,mu,B,E);
    }
     
    inline void computeTemplatedMagneticField(const tIndex& t,
                                              const SM_Network& network,
                                              const tReal *mu,
                                              tReal *B) const {
        tReal E=0;
        computeTemplatedMagneticFieldAndEnergy(t,network,mu,B,E);
        
    }
    inline void computeTemplatedMagneticFieldAndEnergy(const tIndex& t,
                                                       const SM_Network& network,
                                                       const tReal *mu,
                                                       tReal *B,
                                                       tReal& E) const {
        //get the number of particles of the network
        const tIndex& nParticles=network.getParticlesNumber();
 
        //get the dimension of the netwok
        const tDimension& dim=network.getDimension();
 
        //iterator on dimension
        tDimension d;
 
        
        
        //Bext values at time t=0
        const tReal *Bext_i,*Bext_id;
        tBoolean incBext=(mZeemanParticlesNumber>1);
        
        //iterator on B
        tReal *Bi_d=B;
 
        //iterator on Mu
        const tReal *mu_id=mu;
        //energy of the operator
        E=0;
        //Bext at particule i  
        Bext_i=getBext(t);
        
        for (tIndex i=0;i<nParticles;i++) {
 
            //Bext at coordinate d
            Bext_id=Bext_i;
            
            for (d=0;d<dim;d++) {
                
                //zeeman magnetic field contribution to total field
                (*Bi_d)+=(*Bext_id);
                
                //energy of the operator
                E-=(*Bext_id)*(*mu_id);
 
                //next coordinate
                Bext_id++;//Bext value at next coordinate
                Bi_d++;//B value at next coordinate
                mu_id++;//mu value at next coordinate
            }
 
            //Bext at next perticles
            Bext_i+=incBext;
            
        }
 
        
    }
    
     
    virtual  tReal computeEnergy(const SM_Network& network,
                                 const tReal *mu,
                                 const tReal *H) const override {
        //get the size of the network
        tIndex i,N=network.getParticlesNumber()*network.getDimension();
        tReal E=0;
        const tReal *mu_i=mu;
        const tReal *H_i=H;
        for (i=0;i<N;i++) {
            E-=(*mu_i)*(*H_i);
            mu_i++;
            H_i++;
        }
        return E;
    }
    
    virtual  tReal computeEnergy(const tIndex& t,
                                 const SM_Network& network,
                                 const tReal *mu) const override{
        return computeTemplatedEnergy(t,network,mu);
    }
    inline tReal computeTemplatedEnergy(const tIndex& t,
                                        const SM_Network& network,
                                        const tReal *mu) const {
 
        //get the number of particles of the network
        const tIndex& nParticles=network.getParticlesNumber();
 
        //get the dimension of the netwok
        const tDimension& dim=network.getDimension();
 
        //iterator on dimension
        tDimension d;
        
        //Bext values uniform on all points
        const tReal *Bext_id,*Bext_i=getBext(t);
        tBoolean incBext=(mZeemanParticlesNumber>1);
        
        //magnetic moment at point i
        const tReal *Mui_d=mu;
        //energy E=-\sum_i <Bext,Mu_i>
        tReal E=0;
        for (tIndex i=0;i<nParticles;i++) { //loop on particles
            //external field
            Bext_id=Bext_i;
            for (d=0;d<dim;d++) {//loopon dimension
                E-=(*Bext_id)*(*Mui_d);
                
                Bext_id++;
                Mui_d++;
            }
 
            //BNext at next partiule
            Bext_i+=incBext;
        }
 
        return E;
    }//end method
  
};
 
 
#endif