EMM_SLElementaryDemagnetizedMatrix.h

Go to the documentation of this file.

#ifndef EMM_SLElementaryDemagnetizedMatrix_H
#define EMM_SLElementaryDemagnetizedMatrix_H

#include "CORE_Object.h"

DEFINE_SPTR(EMM_SLElementaryDemagnetizedMatrix);

class EMM_SLElementaryDemagnetizedMatrix : public virtual CORE_Object {
  
  
  SP_OBJECT(EMM_SLElementaryDemagnetizedMatrix);
  
  // ATTRIBUTES
  
  
  
 private:
    static const tReal m4GaussianWeights[];
    static const tReal m4GaussianPoints[];
    
 protected: 
  // METHODS
  
  // CONSTRUCTORS 
    
  EMM_SLElementaryDemagnetizedMatrix(void);
   
    
    
  // DESTRUCTORS 
    
    
  virtual ~EMM_SLElementaryDemagnetizedMatrix(void);   
    
   
 private:

  // CONSTRUCTORS
  
 public:
  inline static SP::EMM_SLElementaryDemagnetizedMatrix New() {
    SP::EMM_SLElementaryDemagnetizedMatrix ret(new EMM_SLElementaryDemagnetizedMatrix(),
                                       EMM_SLElementaryDemagnetizedMatrix::Delete()); 
    ret->setThis(ret);
    return ret;
  }; 
 
  
  
 public:
    
 
  // SET methods
    
  // GET methods 

  tReal KaaG(const int& x,
             const int& y,
             const int& z,
             const int& dx,
             const int& dy,
             const int& dz,
             tBoolean& succeeds) const;

  tReal KabG(const int& x,
             const int& y,
             const int& z,
             const int& dx,
             const int& dy,
             const int& dz,
             tBoolean& succeeds) const;
  
  tReal KabDG(const int& dx,
              const int& dy,
              const int& dz,
              tBoolean& succeeds) const;

  tReal KaaDG(const int& dx,
              const int& dy,
              const int& dz,
              tBoolean& suceeds) const;
    

public:
    static inline tReal integrateKxx(const int& index,const int& dx,const int& dy,const int& dz,
                                     const int& r,const int& s,const int& t,
                                     tReal work[4]) {
        tReal int01Kxx=0;
        int i,j,k;
        for (i=0;i<4;i++) {
            for (j=0;j<4;j++) {
                for (k=0;k<4;k++) {
                    work[0]=(m4GaussianPoints[i]-1)/2;
                    work[1]=(m4GaussianPoints[j]-1)/2;
                    work[2]=(m4GaussianPoints[k]-1)/2;
                    int01Kxx+=(m4GaussianWeights[i]*m4GaussianWeights[j]*m4GaussianWeights[k]/8)*
                        Kxx(index,dx,dy,dz,work,r,s,t);
                   
                }
            }
        }
        return int01Kxx;
    }
    static inline tReal integrateKxy(const int& index,const int& dx,const int& dy,const int& dz,
                                     const int& r,const int& s,const int& t,
                                     tReal work[4]) {
        tReal int01Kxy=0;
        int i,j,k;
        for (i=0;i<4;i++) {
            for (j=0;j<4;j++) {
                for (k=0;k<4;k++) {
                    work[0]=(m4GaussianPoints[i]-1)/2;
                    work[1]=(m4GaussianPoints[j]-1)/2;
                    work[2]=(m4GaussianPoints[k]-1)/2;
                    int01Kxy+=(m4GaussianWeights[i]*m4GaussianWeights[j]*m4GaussianWeights[k]/8)*
                        Kxy(index,dx,dy,dz,work,r,s,t);
                }
            }
        }
        return int01Kxy;
    }

    static inline tReal Kxx(const int& k,
                            const int &dx,
                            const int &dy,
                            const int &dz,
                            tReal x[4],
                            const int& r,
                            const int& s,
                            const int& t) {
        x[0]=dx+x[0]+r;
        x[1]=dy+x[1]+s;
        x[2]=dz+x[2]+t;
      
        x[3]=x[k]*sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
        
        tReal kxx=1.;
        if (x[k]!=0) {
            kxx=-atan(x[(k+1)%3]*x[(k+2)%3]/x[3]);
        }else {
            kxx=-M_PI/2;
        }

        if ((r+s+t)%2==1) {//(-1)^{r+s+t}=-1
            kxx*=-1;
        }
        return kxx;
    }
    static inline tReal Kxy(const int& k,
                            const int &dx,
                            const int &dy,
                            const int &dz,
                            tReal x[4],
                            const int& r,
                            const int& s,
                            const int& t) {
        x[0]=dx+x[0]+r;
        x[1]=dy+x[1]+s;
        x[2]=dz+x[2]+t;

        //x[3]=x[0]*x[0]+x[1]*x[1]
        x[3]=x[(k+1)%3];
        x[3]*=x[3];
        x[3]+=x[k]*x[k];
       
        
        tReal kxy=x[(k+2)%3];//kxy=x[2]
        
        if (x[3]>0) {
            //kyx=x[2]+sqrt(x[2]*x[2]+x[0]*x[0]+x[1]*x[1])
            kxy+=sqrt(kxy*kxy+x[3]);
            //kxy=log((x[2]+sqrt(x[2]*x[2]+x[0]*x[0]+x[1]*x[1]))/sqrt(x[0]*x[0]+x[1]*x[1]))
            kxy=log(kxy/sqrt(x[3]));
        } else {
            if (kxy>0) { 
                kxy=log(2*kxy);
            } else {
                kxy=-CORE_Object::getRealInfinity();
            }
        }
       
        if ((r+s+t)%2==1) {//(-1)^{r+s+t}=-1
            kxy*=-1;
        }
       
        return kxy;
    }
    
 public:
  // OTHERS methods
 

 
    

};

#endif