#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>


/* dünner Film mit Ausdehnung 2b in y und 2a in y, zentriert, b>=a
  printf("%e %e\n",1./0.,0./0.); ergibt: 1.#INF00e+00 -1.#IND00e+00
*/

double a,b,d,jc,acc;
int intervalls;
gsl_integration_workspace *w;

double f(double x,double y)
{
  double P,Q,sum1,sum2;
  P=sqrt((a-x)*(a-x)+(b-y)*(b-y));
  Q=sqrt(x*x+(b-a-y)*(b-a-y));

  sum1=sqrt(2.)*log((sqrt(2.)*P+a+b-x-y)/(sqrt(2.)*Q-a+b-x-y));
  sum2=log( fabs( ( (P+y-b)*(y-b+a)*(P+x-a)*x) / ( (y-b)*(Q+y-b+a)*(x-a)*(Q+x) ) ) );
  return sum1+sum2;
}

double Bz(double x,double y)
{
  return 1.e-7*jc*d*(f(x,y)+f(-x,y)+f(x,-y)+f(-x,-y));
}

/*
  overrides the default signal handler, no abort in case of errors
*/
void handler(const char * reason,const char * file,int line,int gsl_errno)
{
 int i;
 printf("%s in %s (%i): Error %i\n",reason,file,line,gsl_errno);
 getchar();
}

void print_field()
{
     int i,j;
     FILE *file;
     file=fopen("log.dat","w");

     for(i=0;i<=100;i++)
     {
             for(j=0;j<=100;j++)
             {
               fprintf(file,"%e\t",Bz(-a+2.*a/100.*j,b-2.*b/100*i));
             }
             fprintf(file,"\n");
     }
     fclose(file);
}

// y-> b-a+y
double int1(double y,void *params)
{
    double result,error;
    int status;
    gsl_function F;
    double func1(double x,void *params)
    {
     return fabs(Bz(x,b-a+y));
    }
    F.function=&func1;
    F.params=NULL;

    status=gsl_integration_qags(&F,0,y,0,acc,intervalls,w,&result,&error);
    return result;
}

double intint1()
{
    double result,error;
    int status;
    gsl_function F;
    F.function=&int1;
    F.params=NULL;

    status=gsl_integration_qags(&F,0,a,0,acc,intervalls,w,&result,&error);
    return result;
}

double int2(double x,void *params)
{
    double result,error;
    int status;
    gsl_function F;
    double func2(double y,void *params)
    {
     return fabs(Bz(x,y));
    }
    F.function=&func2;
    F.params=NULL;

    status=gsl_integration_qags(&F,0,b-a+x,0,acc,intervalls,w,&result,&error);
    return result;
}

double intint2()
{
    double result,error;
    int status;
    gsl_function F;
    F.function=&int2;
    F.params=NULL;

    status=gsl_integration_qags(&F,0,a,0,acc,intervalls,w,&result,&error);
    return result;
}

double g(double *k, size_t dim, void *params)
{
    return fabs(Bz(k[0],k[1]));
}

int main()
{
    gsl_error_handler_t *old_handler;
    int i,status;
    double dummy,result,error;
    FILE *file;
    a=0.0001;
    b=0.0002;
    d=0.000001;
    jc=1.e9;
    intervalls=1e3;
    acc=1.e-6;

    double xl[2] = {0,0};
    double xu[3] = {a,b};
    const gsl_rng_type *T;
    gsl_rng *r;
    gsl_monte_function G={&g,2,0};
    size_t calls = 500000;

    // override handler, allocate workspace memory
    old_handler=gsl_set_error_handler(&handler);
    w=gsl_integration_workspace_alloc(intervalls);

    file=fopen("file.log","w");

    printf("%e\n",(intint2()+intint1())/(a*b));


    gsl_rng_env_setup ();
    T=gsl_rng_default;
    // still alive
    printf("Allocating workspace for montecarlo integration\n");
    r=gsl_rng_alloc (T);
    printf("Dead\n");
    // dead
    {
     gsl_monte_miser_state *s= gsl_monte_miser_alloc (3);
     gsl_monte_miser_integrate(&G,xl,xu,3,calls,r,s,&result,&error);
     gsl_monte_miser_free (s);
    }
    printf("monte carlo: %e\n",result/(a*b));

    getchar();
}