SM_DipolarField.h

#ifndef SM_DipolarField_H
#define SM_DipolarField_H
 
//base classes
#include "SM_Object.h"
 
 
class SM_DipolarField : public virtual SM_Object  {
    
 
public:
private :
    typedef SM_DipolarField SelfClass;
    typedef SM_Object SuperClass;
    //multiple associations
    //---------------------
 
    
   
protected:
    // CONSTRUCTORS 
    SM_DipolarField(void) {
    }
    
    // DESTRUCTORS 
    virtual ~SM_DipolarField(void) {
    }
   
public:
    
    inline static CORE_UniquePointer<SelfClass> New() {
        return CORE_UniquePointer<SelfClass>(new SelfClass(),
                                             SelfClass::Delete());
    }
   
    inline static void ComputeUnscaledSpinDipolarField(const tReal *Xi,
                                                       const tReal* sXj,const tReal *eXj,
                                                       const tReal *sSj,
                                                       std::array<tReal,SM_Constants::DIM>& Hi,
                                                       std::array<tReal,SM_Constants::DIM>& U) {
        
        //iteratoron coordinate at particle i
        const tReal *iXi=null;
        
   
        //iterator on U which is the normalize direction U=Xj-Xi
        tReal *iU=U.data();
        const tReal *eU=iU;eU+=SM_Constants::DIM;
        
        tReal rij,rij3;//distance between particle i & j
 
        
        //scalar product 
        tReal sUSj;//<U,Sj>
 
        //iterators on coordinates of points j
        const tReal *iXj=sXj;
 
        //iterators on dierction of spin of points j
        const tReal *iSj=sSj;
        const tReal *iSjk=null;
 
        //iterator on Hi
        tReal *iHi=Hi.data();
 
        //Hi=0
        functions_array::reset(Hi);
        
        while (iXj!=eXj) {//loop on particle j
            
            if (iXj!=Xi) {//particles j must be different from particle i
                
                //|Xj-Xi|^3
                rij=0;
                
                //iterator on X at particle i
                iXi=Xi;
                //iterator on U=Xj-Xi
                iU=U.data();
                while (iU!=eU) {
                    
                    //U=Xj-Xi
                    (*iU)=(*iXj)-(*iXi);
                    
                    //rij=|Xj-Xi|^2
                    rij+=(*iU)*(*iU);
                    
                    //iterators at next coordinate
                    iU++;
                    iXj++;
                    iXi++;
                    
                }
                //|X_i-X_j|^2
                rij3=rij;
                
                //|X_i-X_j|
                rij=sqrt(rij);
 
               
                //|X_i-X_j|^3
                rij3*=rij;
                
                //Computes:
                //U=U/|U|
                
                //sUSj=<U,Sj>
                sUSj=0;
                
                
                //iterator on Pj-Pi
                iU=U.data(); 
                
                //iterator on Sj
                iSjk=iSj;
                while (iU!=eU) {//loop on coordinate
                    
                    //er=Xj-Xi/|Xj-Xi|
                    (*iU)/=rij;
                    
                    
                    //<U,Sj>
                    sUSj+=(*iSjk)*(*iU);
                    
                    //iterator on next coordinate
                    iU++;
                    iSjk++;
                    
                }
                
                //iterator on U=Pj-Pi/|Pj-Pi|
                iU=U.data();
                //iterator on H at point i
                iHi=Hi.data();
                
                while (iU!=eU) {//loop on coordinate
                    
                    //Hi[d]=\sum_j (3 <U,Sj>U-Sj)/r^3
                    (*iHi)+=(3*sUSj*(*iU)-(*iSj))/rij3;
                   
                    
                    //iterator at next coordinates
                    iHi++;
                    iSj++;
                    iU++;
                }
              
                
            }//end condition  Xi!=Xj
            else {
                //next particle j
                iXj+=SM_Constants::DIM;
                iSj+=SM_Constants::DIM;
            }
        }//end loop on particle j
        
 
    }
    inline static void ComputeUnscaledSpinDipolarFieldInSeparatedNetworkStorage(const tReal *Xi,
                                                                                const tReal* sYj,const tReal *eYj,
                                                                                const tReal *sSj,
                                                                                std::array<tReal,SM_Constants::DIM>& Hi,
                                                                                std::array<tReal,SM_Constants::DIM>& U) {
        //zero
        tReal eps=1.e-9;
        
        //iteratoron coordinate at particle i
        const tReal *iXi=null;
        
   
        //iterator on U which is the normalize direction U=Xj-Xi
        tReal *iU=U.data();
        const tReal *eU=iU;eU+=SM_Constants::DIM;
        
        tReal rij,rij3;//distance between particle i & j
 
        
        //scalar product 
        tReal sUSj;//<U,Sj>
 
        //iterators on coordinates of points j
        const tReal *iYj=sYj;
 
        //iterators on dierction of spin of points j
        const tReal *iSj=sSj;
        const tReal *iSjk=null;
 
        //iterator on Hi
        tReal *iHi=Hi.data();
 
        //Hi=0
        memset(iHi,0,sizeof(tReal)*SM_Constants::DIM);
        
        while (iYj!=eYj) {//loop on particle j
            
           
            
            //|Xj-Xi|^3
            rij=0;
            
            //iterator on X at particle i
            iXi=Xi;
            
            //iterator on U=Xj-Xi
            iU=U.data();
            while (iU!=eU) {
                
                //U=Yj-Xi
                (*iU)=(*iYj)-(*iXi);
                
                //rij=|Xj-Xi|^2
                rij+=(*iU)*(*iU);
                
                //iterators at next coordinate
                iU++;
                iYj++;
                iXi++;
                
            }
            //|Y_j-X_i|^2
            rij3=rij;
            
            //|Y_j-X_i|
            rij=sqrt(rij);
 
            //std::cout<<"Xi="<<functions_array::toString(3,Xi)<<" Yj="<<functions_array::toString(3,iYj-3)<<" rij="<<rij<<"\n";
            if (rij>eps) {
                
                //|Y_j-X_i|^3
                rij3*=rij;
                
                //Computes:
                //U=U/|U|
                
                //sUSj=<U,Sj>
                sUSj=0;
                
                
                //iterator on Yj-Xi
                iU=U.data(); 
                
                //iterator on Sj
                iSjk=iSj;
                while (iU!=eU) {//loop on coordinate
                    
                    //Yj-Xi/|Yj-Xi|
                    (*iU)/=rij;
                    
                    
                    //<U,Sj>
                    sUSj+=(*iSjk)*(*iU);
                    
                    //iterator on next coordinate
                    iU++;
                    iSjk++;
                    
                }
                
                //iterator on U=Yj-Xi/|Yj-Xi|
                iU=U.data();
                //iterator on H at point i
                iHi=Hi.data();
               
                while (iU!=eU) {//loop on coordinate
                    
                    //Hi[d]=\sum_j (3 <U,Sj>U-Sj)/r^3
                    (*iHi)+=(3*sUSj*(*iU)-(*iSj))/rij3;
                  
                    
                    //iterator at next coordinates
                    iHi++;
                    iSj++;
                    iU++;
                }
               
                
            }//end condition  rij>0
            else {
                //next particle j
                iSj+=SM_Constants::DIM;
            }
        }//end loop on particle j
        
 
    }
    
    
    inline static void ComputeSeparatedUnscaledSpinDipolarField(const tReal *Xi,
                                                                const tReal* sXj,const tReal *eXj,
                                                                const tReal *sSj,
                                                                std::array<tReal,SM_Constants::DIM>& Hi,
                                                                std::array<tReal,SM_Constants::DIM>& U) {
        
        //iteratoron coordinate at particle i
        const tReal *iXi=null;
        
        
        //iterator on U which is the normalize direction U=Xj-Xi
        tReal *iU=U.data();
        const tReal *eU=iU;eU+=SM_Constants::DIM;
        
        tReal rij,rij3;//distance between particle i & j
        
        
        //scalar product 
        tReal sUSj;//<U,Sj>
        
        //iterators on coordinates of points j
        const tReal *iXj=sXj;
        
        //iterators on dierction of spin of points j
        const tReal *iSj=sSj;
        const tReal *iSjk=null;
        
        //iterator on Hi
        tReal *iHi=Hi.data();
        
        //Hi=0
        memset(iHi,0,sizeof(tReal)*SM_Constants::DIM);
        
        while (iXj!=eXj) {//loop on particle j
            
            
            
            //|Xj-Xi|^3
            rij=0;
            
            //iterator on X at particle i
            iXi=Xi;
            //iterator on U=Xj-Xi
            iU=U.data();
            while (iU!=eU) {
                
                //U=Xj-Xi
                (*iU)=(*iXj)-(*iXi);
                
                //rij=|Xj-Xi|^2
                rij+=(*iU)*(*iU);
                
                //iterators at next coordinate
                iU++;
                iXj++;
                iXi++;
                
            }//end loop on U coordinates
            
            //|X_i-X_j|^2
            rij3=rij;
            
            //|X_i-X_j|
            rij=sqrt(rij);
            
            
            //|X_i-X_j|^3
            rij3*=rij;
            
            //Computes:
            //U=U/|U|
            
            //sUSj=<U,Sj>
            sUSj=0;
            
            
            //iterator on Pj-Pi
            iU=U.data(); 
            
            //iterator on Sj
            iSjk=iSj;
            while (iU!=eU) {//loop on U coordinate
                
                //Xj-Xi/|Xj-Xi|
                (*iU)/=rij;
                
                
                //<U,Sj>
                sUSj+=(*iSjk)*(*iU);
                
                //iterator on next coordinate
                iU++;
                iSjk++;
                
            }//end llo on U coordinate
            
            //iterator on U=Pj-Pi/|Pj-Pi|
            iU=U.data();
            //iterator on H at point i
            iHi=Hi.data();
            while (iU!=eU) {//loop on U coordinate
                
                //Hi[d]=\sum_j (3 <U,Sj>U-Sj)/r^3
                (*iHi)+=(3*sUSj*(*iU)-(*iSj))/rij3;
                
                
                //iterator at next coordinates
                iHi++;
                iSj++;
                iU++;
            }//end loop on U coordinate
            
            
        }//end loop on particle j
    }
        
    
    
    inline static void ComputeUnscaledSpinEnergy(const tReal *Xi,const tReal *Si,
                                                 const tReal* iXj,const tReal* eXj,const tReal* iSj,
                                                 tReal& Ei) {
 
      
        //iteratoron coordinate at particle i
        const tReal *iXi=null;
        
        //iterator on Si
        const tReal *iSi=null;
        
        
        
        //iterator on U which is the normalize direction U=Xj-Xi
        std::array<tReal,SM_Constants::DIM> U;
        tReal *iU=U.data();
        const tReal *eU=iU;eU+=SM_Constants::DIM;
        
        tReal rij,rij3;//distance between particle i & j
 
        
        //scalar product 
        tReal sUSj;//<U,Sj>
        tReal sUSi;//<U,Si>
        tReal sSiSj;//<U,Si>
        
        Ei=0;
        
        
        while (iXj!=eXj) {//loop on particle j
            
            if (iXj!=Xi) {//particles j must be different from particle i
                
                //|Pj-Pi|^3
                rij=0;
                
                //iterator on X at particle i
                iXi=Xi;
                //iterator on U=Xj-Xi
                iU=U.data();
                while (iU!=eU) {
                    
                    //U=Xj-Xi
                    (*iU)=(*iXj)-(*iXi);
                    
                    //rij=|Xj-Xi|^2
                    rij+=(*iU)*(*iU);
                    
                    //iterators at next coordinate
                    iU++;
                    iXj++;
                    iXi++;
                    
                }
                //|X_i-X_j|^2
                rij3=rij;
                
                //|X_i-X_j|
                rij=sqrt(rij);
                
                //|X_i-X_j|^3
                rij3*=rij;
                
                //Computes:
                //U=U/|U|
                
                //sUSj=<U,Sj>
                sUSj=0;
                
                //sUSi=<U,Si>
                sUSi=0;
                
                //sUSi=<Si,Sj>
                sSiSj=0;
                
                //iterator on Pj-Pi
                iU=U.data(); 
                
                //iterator on Si
                iSi=Si;
                while (iU!=eU) {//loop on coordinate
                    
                    //Xj-Xi/|Xj-Xi|
                    (*iU)/=rij;
                    
                    
                    //<U,Sj>
                    sUSj+=(*iSj)*(*iU);
                    
                    //<U,Si>
                    sUSi+=(*iSi)*(*iU);
                    
                    //<Si,Sj>
                    sSiSj+=(*iSi)*(*iSj);
                    
                    //iterator on next coordinate
                    iU++;
                    iSj++;
                    iSi++;
                    
                }
                
                  
                Ei+=(sSiSj-3*sUSi*sUSj)/rij3;
                
            }//end condition  Xi!=Xj
            else {
                //next particle j
                iXj+=SM_Constants::DIM;
                iSj+=SM_Constants::DIM;
            }
        }//end loop on particle j
        
 
    }
    inline static void ComputeSeparatedUnscaledSpinEnergy(const tReal *Xi,const tReal *Si,
                                                          const tReal* iXj,const tReal* eXj,const tReal* iSj,
                                                          tReal& Ei) {
 
      
        //iteratoron coordinate at particle i
        const tReal *iXi=null;
        
        //iterator on Si
        const tReal *iSi=null;
        
        
        
        //iterator on U which is the normalize direction U=Xj-Xi
        std::array<tReal,SM_Constants::DIM> U;
        tReal *iU=U.data();
        const tReal *eU=iU;eU+=SM_Constants::DIM;
        
        tReal rij,rij3;//distance between particle i & j
 
        
        //scalar product 
        tReal sUSj;//<U,Sj>
        tReal sUSi;//<U,Si>
        tReal sSiSj;//<U,Si>
        
        Ei=0;
        
        
        while (iXj!=eXj) {//loop on particle j
         
            
            //|Pj-Pi|^3
            rij=0;
            
            //iterator on X at particle i
            iXi=Xi;
            //iterator on U=Xj-Xi
            iU=U.data();
            while (iU!=eU) {
                
                //U=Xj-Xi
                (*iU)=(*iXj)-(*iXi);
                
                //rij=|Xj-Xi|^2
                rij+=(*iU)*(*iU);
                
                //iterators at next coordinate
                iU++;
                iXj++;
                iXi++;
                
            }
            //|X_i-X_j|^2
            rij3=rij;
            
            //|X_i-X_j|
            rij=sqrt(rij);
            
            //|X_i-X_j|^3
            rij3*=rij;
            
            //Computes:
            //U=U/|U|
            
            //sUSj=<U,Sj>
            sUSj=0;
            
            //sUSi=<U,Si>
            sUSi=0;
            
            //sUSi=<Si,Sj>
            sSiSj=0;
            
            //iterator on Pj-Pi
            iU=U.data(); 
            
            //iterator on Si
            iSi=Si;
            while (iU!=eU) {//loop on coordinate
                
                //Xj-Xi/|Xj-Xi|
                (*iU)/=rij;
                
                
                //<U,Sj>
                sUSj+=(*iSj)*(*iU);
                
                //<U,Si>
                sUSi+=(*iSi)*(*iU);
                
                //<Si,Sj>
                sSiSj+=(*iSi)*(*iSj);
                    
                //iterator on next coordinate
                iU++;
                iSj++;
                iSi++;
                    
            }
            
            
            Ei+=(sSiSj-3*sUSi*sUSj)/rij3;
            
       
        }//end loop on particle j
        
 
    }
  
    inline static void ComputeSeparatedUnscaledSpinEnergy(const tReal* iXi,const tReal *eXi,const tReal *iSi,
                                                          const tReal* iXj,const tReal* eXj,const tReal* iSj,
                                                          tReal& E) {
 
      
        
        
     
        E=0;
        tReal Ei=0;
        while (iXi!=eXi) {//loop on particles i
 
            //compute the spin energy of the spin i
            Ei=0;
            ComputeSeparatedUnscaledSpinEnergy(iXi,iSi,iXj,eXj,iSj,Ei);
 
            E+=Ei;
            //next spin
            iXi+=SM_Constants::DIM;
            iSi+=SM_Constants::DIM;
        }
        
        
 
    }
 
    inline static void ComputeUnscaledSpinEnergy(const tReal* iXi,const tReal *eXi,const tReal *iSi,
                                                         const tReal* iXj,const tReal* eXj,const tReal* iSj,
                                                         tReal& E) {
 
      
        
        
     
        E=0;
        tReal Ei=0;
        while (iXi!=eXi) {//loop on particles i
 
            //compute the spin energy of the spin i
            Ei=0;
            ComputeUnscaledSpinEnergy(iXi,iSi,iXj,eXj,iSj,Ei);
 
            E+=Ei;
            //next spin
            iXi+=SM_Constants::DIM;
            iSi+=SM_Constants::DIM;
        }
        
        
 
    }
 
  
};
 
 
#endif