#include "thinfilm_BLS.h"
#include <math.h>

void tf_buildmatrix(double **M, tf_Layer * l, tf_Global g)
{
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//Rezende and Cochran use positive j2 for biquadratic coupling !!
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

int i,j,n;
double cos_dth, cos_2dth;

n=g.n_layers;

//set everthing to zero
for(i=1;i<=2*n;i++)
  for(j=1;j<=2*n;j++)
      M[i][j] = 0;

//pointers to matrix for better readabillity:

#define Mxy(a,b) M[2*(a)-1][2*(b)]
#define Myx(a,b) M[2*(a)][2*(b)-1]
#define Mxx(a,b) M[2*(a)-1][2*(b)-1]
#define Myy(a,b) M[2*(a)][2*(b)]

//now build the matrix

//the gyromagnetic part:
for(i=1;i<=n;i++)
  {  Mxx(i,i) +=  -g.omega/l[i].gamma; 
     Myy(i,i) +=  -g.omega/l[i].gamma; }

//the zeeman terms:
for(i=1;i<=n;i++)
   Mxy(i,i) += g.B * cos(l[i].theta_M - g.theta_B), //Rezende: H1 
   Myx(i,i) += g.B * cos(l[i].theta_M - g.theta_B); //H3

//EB
for(i=1;i<=n;i++)
   Mxy(i,i) += l[i].X/l[i].M * cos(l[i].theta_M - l[i].theta_X), //H1
   Myx(i,i) += l[i].X/l[i].M * cos(l[i].theta_M - l[i].theta_X); //H3

//the interlayer exchange coupling 
//with the upper film
for (i=1;i<n;i++)
  cos_dth = cos(l[i].theta_M -l[i+1].theta_M), 
  cos_2dth = cos(2*(l[i].theta_M -l[i+1].theta_M)),
   Mxy(i,i) += 1.0 /l[i].M /l[i].d 				//H1
		*(l[i].j1*cos_dth + 2*l[i].j2*sqr(cos_dth)),
   Myx(i,i) += 1.0 /l[i].M /l[i].d 				//H3 
		*(l[i].j1*cos_dth +2*l[i].j2*cos_2dth),
   Mxy(i,i+1) += 1.0/l[i+1].M /l[i+1].d				//H2
	        *(-l[i].j1 - 2*l[i].j2*cos_dth),
   Myx(i,i+1) += 1.0/l[i+1].M /l[i+1].d				//H4
 		*(-l[i].j1*cos_dth -2*l[i].j2*cos_2dth);

//with the lower film:
for (i=2;i<=n;i++)
  cos_dth = cos(l[i].theta_M -l[i-1].theta_M),
  cos_2dth = cos(2*(l[i].theta_M -l[i-1].theta_M)),
   Mxy(i,i) += 1.0 /l[i].M /l[i].d                      	//H1
                *(l[i-1].j1*cos_dth + 2*l[i-1].j2*sqr(cos_dth)),
   Myx(i,i) += 1.0 /l[i].M /l[i].d                      	//H3
                *(l[i-1].j1*cos_dth + 2*l[i-1].j2*cos_2dth),
   Mxy(i,i-1) += 1.0/l[i-1].M /l[i-1].d                 	//G2
                *(-l[i-1].j1 - 2*l[i-1].j2*cos_dth),
   Myx(i,i-1) += 1.0/l[i-1].M /l[i-1].d                		//G4
                *(-l[i-1].j1*cos_dth -2*l[i-1].j2*cos_2dth);

//end interlayer exchange coupling

//the fourfold bulk anisotropy
for(i=1;i<=n;i++)
   Mxy(i,i) += 							//H1
	0.5*l[i].K/l[i].M*(3+cos(4*(l[i].theta_M-l[i].theta_K))),
   Myx(i,i) += 							//H3
	2*l[i].K/l[i].M*cos(4*(l[i].theta_M - l[i].theta_K));

//the uniaxial bulk anisotropy
for(i=1;i<=n;i++)
   Mxy(i,i) += 2 * l[i].U/l[i].M *
	sqr(cos (l[i].theta_M - l[i].theta_U)), 
   Myx(i,i) += 2 * l[i].U/l[i].M *
	cos (2*(l[i].theta_M - l[i].theta_U));

//the surface anisotropy
for(i=1;i<=n;i++)
        Mxy(i,i) += -2 * l[i].S / l[i].M / l[i].d;

//the thinfilm intralayer exchange:
for(i=1;i<=n;i++)
   Mxy(i,i) += l[i].D*sqr(g.q),
   Myx(i,i) += l[i].D*sqr(g.q);

if(perp_fmr !=1)
//dipolar interactions
//can be further optimized for the special cases:
//q=0 (FMR); q*d_i << 1; q*d_total <<1  ;rigid magnetization; ...
{
double **exp_sq, *sin_dMq;
exp_sq = matrix(1,n,1,n);
sin_dMq = vector(1,n);

for(i=1;i<=n;i++)
  {sin_dMq[i] = 0;	 
   for(j=1;j<=n;j++)
	     exp_sq[i][j] = 0;
 }

//calculate spacing:
for(i=1;i<n;i++)
  {
  exp_sq[i][i+1] = +l[i].s;
  for(j=i+2;j<=n;j++)
     exp_sq[i][j] = exp_sq[i][j-1] + l[j-1].s + l[j-1].d;
  }  

for(i=1;i<n;i++)
  for(j=i+1;j<=n;j++) 
{	exp_sq[i][j] = exp(-g.q*exp_sq[i][j]);
	exp_sq[j][i] = exp_sq[i][j];
}
		
for(i=1;i<=n;i++)
	sin_dMq[i] = sin(l[i].theta_M - g.theta_q);



//the single layer terms:
for(i=1;i<=n;i++)
   Mxy(i,i) += mu0*l[i].M*(1 - 0.5*g.q*l[i].d),
   Myx(i,i) += 0.5*mu0*l[i].M*g.q*l[i].d*sqr(sin_dMq[i]);


//the dipolar coupling:
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
{
	if(i != j)
	{ Mxy(i,j) += -0.5*mu0*l[i].M*g.q*l[i].d		//H2
			*(1 - 0.5*g.q*l[j].d)*exp_sq[i][j];	
          Myx(i,j) += 0.5*mu0*l[i].M*g.q*l[i].d			//H4
			*(1 - 0.5*g.q*l[j].d)*exp_sq[i][j]
			*sin_dMq[i]*sin_dMq[j];
       	 Mxx(i,j) +=   						//H5,G5
		Sign(j-i)*0.5*mu0*l[i].M*g.q*l[i].d
		*(1- 0.5*g.q*l[j].d)*exp_sq[i][j]
		*sin_dMq[j];
         Myy(i,j) +=  						//H6,G6
		-Sign(j-i)*0.5*mu0*l[i].M*g.q*l[i].d
                *(1- 0.5*g.q*l[j].d)*exp_sq[i][j]
                *sin_dMq[i];	

       }
}

free_vector(sin_dMq,1,n);
free_matrix(exp_sq,1,n,1,n);

} //end dipolar interaction

}

void tf_updateM(double **M, tf_Layer * l, tf_Global g)
{
int i,n;
n=g.n_layers;

//pointers to matrix for better readabillity:

#define Mxy(a,b) M[2*(a)-1][2*(b)]
#define Myx(a,b) M[2*(a)][2*(b)-1]
#define Mxx(a,b) M[2*(a)-1][2*(b)-1]
#define Myy(a,b) M[2*(a)][2*(b)]

//now build the matrix

//the gyromagnetic part:
for(i=1;i<=n;i++)
  {  Mxx(i,i) =  -g.omega/l[i].gamma;
     Myy(i,i) =  -g.omega/l[i].gamma; }
}
