Calculate demagnetization coefficient.

demag.c

/*
  demag.c
*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define LENG 100

static double PI = 3.141592653582;

double f(double t, double a, double b, double c);
double max_three(double a, double b, double c);
double DE_integral (double a, double b, double c, double err, double h, int J_max);

int main ()
{
  char input_a[LENG],input_b[LENG],input_c[LENG];
  int J_max, leng;  
  double a, b, c, max_abc;
  double N_a, N_b, N_c, h;
  double err;
  
  leng = LENG;
  a = 1.0;
  b = 1.0;
  c = 1.0;
  max_abc = 1.0;
  
  fprintf(stdout,"Input a, b, c axes length of a ellipsoid.\n");

  fprintf(stdout,"a = ");
  fgets(input_a, LENG, stdin);    
  sscanf(input_a, "%lf", &a);

  fprintf(stdout,"b = ");
  fgets(input_b, LENG, stdin);    
  sscanf(input_b, "%lf", &b);

  fprintf(stdout,"c = ");
  fgets(input_c, LENG, stdin);    
  sscanf(input_c, "%lf", &c);

  fprintf(stdout,"\na = %f b = %f, c = %f\n\n", a, b, c);

  if (a <= 0 || b<= 0 || c<=0) {
    fprintf(stdout, "Do not enter NEGATIVE numbers or ZERO!!\n");
    exit(1);
  }
  
  max_abc = max_three(a, b, c);
  a = a/max_abc;
  b = b/max_abc;
  c = c/max_abc;
  
  h = 0.001;
  J_max = 8000;
  err = 1.0e-05;
  
  N_a = DE_integral(a, b, c, h, err, J_max);
  N_b = DE_integral(b, c, a, h, err, J_max);
  N_c = DE_integral(c, a, b, h, err, J_max);
  
  fprintf(stdout,"\nN_a/4*pi = %.10f\n", N_a);
  fprintf(stdout,"N_b/4*pi = %.10f\n", N_b);
  fprintf(stdout,"N_c/4*pi = %.10f\n\n", N_c);
  
  fprintf(stdout,"N_a = %.10f\n", 4.0*PI*N_a);
  fprintf(stdout,"N_b = %.10f\n", 4.0*PI*N_b);
  fprintf(stdout,"N_c = %.10f\n", 4.0*PI*N_c);

  return (0);
}

double f(double t, double a, double b, double c) {
  double temp;
  
  temp = cosh(t)*exp(PI*sinh(t)/2.0)/((exp(PI*sinh(t)/2.0)+a*a)*sqrt((exp(PI*sinh(t)/2.0)+a*a)*(exp(PI*sinh(t)/2.0)+b*b)*(exp(PI*sinh(t)/2.0)+c*c)));
  
  return (temp);
}

double max_three(double a, double b, double c) {
  double temp1, temp2;
  if (a > b)
    temp1 = a;
  else
    temp1 = b;
  
  if (temp1 > c)
    temp2 = temp1;
  else
    temp2 = c;
  
  return (temp2);
}

double DE_integral (double a, double b, double c, double h, double err, int J_max){
  double sum,temp;
  int m, n, j, k;
  
  sum = 0.0;
  m = J_max+1;
  n = J_max+1;
  
  for (j = 1; j <= J_max; j++) {
    if (fabs(f(-j*h, a, b, c)) < err) {
      m = j;
      break;
    }
  }
  
  for (j = 1; j <= J_max; j++) {
    if (fabs(f(j*h, a, b, c)) < err) {
      n = j;
      break;
    }
  }
  
  for (k = -m ; k <= n; k++){
    temp= f(k*h, a, b, c);
    sum = sum + temp*PI*h*a*b*c/4.0;
  }

  if (m <= J_max && n<= J_max){
    fprintf(stdout, "The caliculation converged.\n");
  }else{
    fprintf(stdout, "The caliculation did not convergent.\n");
  }
  
  return (sum);
}