Hi list!
The attached program crashes after allocating memory (marked "still alive"
int the file). Ihave no idea why, allocating memory for the classical
integrations
works fine. I'm using dev-cpp 4.9 and gsl 1.6.
Any ideas?
LG Flo
#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();
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
