#include<iostream>
#include<math.h>
#include<stdlib.h>
using namespace std;

/* Calcul du rotationnel et de la divergence d'un champ de vecteurs */

const int N = 3;

void un_champ(double* x, double* v){
  v[0] = x[0] + x[1] + x[2];
  v[1] = x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
  v[2] = sqrt(x[0]) + sqrt(x[1]) + sqrt(x[2]);
}

void rot_div(void (*champ)(double*, double*), double* x, double h,
	     double* rot, double* div){
  double *v0, *v1;
  double** der;
  int i, j;
  
  v0 = (double*)malloc(N * sizeof(double));
  v1 = (double*)malloc(N * sizeof(double));
  der = (double**)malloc(N * sizeof(double*));
  for (i = 0; i < N; i++){
    der[i] = (double*)malloc(N * sizeof(double));
    x[i] = x[i] - h;
    (*champ)(x,v0);
    x[i] = x[i] + 2*h;
    (*champ)(x,v1);
    x[i] = x[i] - h;
    for (j = 0; j < N; j++)
      der[i][j] = (v1[j] - v0[j]) / (2*h);
  }
  *div = 0.;
  for (i = 0; i < N; i++){
    *div += der[i][i];
    rot[i] = der[(i+1)%N][(i+2)%N] - der[(i+2)%N][(i+1)%N];
  }
  for (i = 0; i < N; i++)
    free(der[i]);
  free(der);
  free(v0);
  free(v1);
}

int main(){
  double h = 1.e-10;
  double *x, *rot;
  double div;
  int i;

  x = (double*)malloc(N * sizeof(double));
  x[0] = 1.; x[1] = 2.; x[2] = 3.;
  rot = (double*)malloc(N * sizeof(double));
  rot_div(un_champ, x, h, rot, &div);
  cout << "Divergence : " << div << endl;
  cout << "Rotationnel : ";
  for (i = 0; i < N; i++)
    cout << rot[i] << " ";
  cout << endl;
  free(x);
  free(rot);
  return 0;
}