eccm.c

Calculate the magnetic domain structures in exchange coupled magnetic media.

#include <stdio.h>
#include <math.h>
#define Pi 3.1415926
#define a  3.7e-10
#define N 10000
#define D 1e-4
 int main ()
{
FILE *fp;
fp=fopen("energy.data","w");
float k1,k2, M, A, H ;   
k1=1e7   ;
k2=4e5   ;
M=5e-1    ;
A=1e-11   ;
H=1e-2    ; 
float deltak=0;
int k, n;                                      
//printf ( "%s\n","Please Input number of steps" );
//scanf("%d",&N);
k=N-1;
printf ("%s\n","Please Input number of cells");
scanf("%d",&n);
deltak= (k2-k1)/n;
double theta[n], z[n],lambda,E_theta[n],K[n] ;
double  E[N],beta;
static float deltaE[N-1];
printf ("%s\n","Please Input lambda");
scanf("%f",&lambda);
printf("%s\n","Please Input beta:");
scanf("%f",&beta);

for(int i=0;i<=n;i++ )
 {
  E_theta[i]=0;
  K[i]=k1+deltak*i; 
  if (i==n)
  theta[i]=Pi/2;
  else
  theta[i]=Pi*i/(2*n);/*set initial values for n thetas and positions*/
 
    }
for(int j=0;j<N;j++)
     {  
   //printf("j now is %d",j);

       E[j]=0; 
       for(int i=0;i<n;i++)
          {
            if (i==0)
            {E_theta[i]=(K[i]*sin(2*theta[i])-A*sin(theta[i+1]-theta[i])/(a*a)+H*M*sin(theta[i]))*D;
            theta[i]=theta[i]-lambda*E_theta[i];
            }
            else
            {(E_theta[i]=K[i]*sin(2*theta[i])-A*(sin(theta[i+1]-theta[i])-sin(theta[i]-theta[i-1]))/(a*a)+H*M*sin(theta[i]))*D;
            theta[i]=theta[i]-lambda*E_theta[i];
            }
           }
        for(int i=0;i<n;i++)
            {
              E[j]=E[j]+(K[i]*sin(theta[i])*sin(theta[i])-A*cos(theta[i+1]-theta[i])-H*M*cos(theta[i]))*D;
         //  printf("j is %d",j);     
         // printf(" *** %e\n", E[j]);
    
               }
             if (j>=5)
              { for (int i=0;i<=4;i++)
                 { 
                   k=i; 
                   deltaE[k]=E[j-i]-E[j-1-i];
                        }
                    if(fabs(deltaE[0])<=beta&&fabs(deltaE[1])<=beta&&fabs(deltaE[2])<=beta&&fabs(deltaE[3])<=beta&&fabs(deltaE[4])<=beta)
                   { 
                    // break; 
                     printf("%e, %e,  %e, %e\n",E[j],E[j-1],E[j-2],E[j-3]);
                            }
                printf("The Energy is %e\n",E[j]);
          //      fprintf(fp, "The Energy is %e\n",E[j]);
    
                    //  else
                    //    continue;
                              }
                
    fprintf(fp, "The Energy is %e\n",E[j]);
}
                             
   }