This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: about QAWF integration -- Fourier integrals
- From: xie qiong <xieq at yahoo dot com>
- To: gsl-discuss at sources dot redhat dot com
- Date: Sat, 9 Mar 2002 15:36:50 -0800 (PST)
- Subject: Re: about QAWF integration -- Fourier integrals
Hello,
My entire code is as follows:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_integration.h>
#define PI 3.14159265358979
#define TOL 1e-6
#define OMEGAMAX 100
#define OMEGASPACE 0.1
#define FACTOR 10
#define ARRAYSIZE OMEGAMAX*FACTOR
double myfunc(double x, void *params)
{
double a = *(double *)params;
return exp(- log(x) * log(x)/(2.0 * a * a)) /
(sqrt(2.0 * PI) * a * x);
}
int main ()
{
int i, m;
double w0, error1 = 0.00, error2 = 0.00;
double imag[ARRAYSIZE] = {0.0};
double real[ARRAYSIZE] = {0.0};
double a = 0.23026 * 1.0;
FILE *fp;
gsl_integration_workspace * wp;
gsl_integration_workspace * c_wp;
gsl_integration_qawo_table * wf =
gsl_integration_qawo_table_alloc(0.0,
10000.0,
GSL_INTEG_COSINE,
10000);
gsl_function F;
F.function = &myfunc;
F.params = &a;
for ( w0 = 0.0, m = 0; w0 < OMEGAMAX; w0 +=
OMEGASPACE, m++ ) {
wp = gsl_integration_workspace_alloc(20000);
c_wp = gsl_integration_workspace_alloc(20000);
gsl_integration_qawo_table_set(wf,
w0,
10000.0,
GSL_INTEG_COSINE);
gsl_integration_qawf(&F,
1e-30,
TOL,
20000,
wp,
c_wp,
wf,
real+m,
&error1);
printf("\t\t\t\t%e\n", error1);
printf("%e", real[m]);
gsl_integration_workspace_free(wp);
gsl_integration_workspace_free(c_wp);
wp = gsl_integration_workspace_alloc(20000);
c_wp = gsl_integration_workspace_alloc(20000);
gsl_integration_qawo_table_set(wf, w0,
10000.0, GSL_INTEG_SINE);
gsl_integration_qawf(&F, 1e-30, TOL, 20000,
wp, c_wp, wf, imag+m, &error2);
printf(" %e\n", imag[m]);
gsl_integration_workspace_free(wp);
gsl_integration_workspace_free(c_wp);
}
gsl_integration_qawo_table_free(wf);
fp = fopen("myfttesti.txt", "w");
if ( fp == NULL ) {
printf("Cann't open trapzi10.txt file\n");
return 1;
}
for (i = 0; i < ARRAYSIZE; i++)
fprintf(fp, "%e\n", imag[i]);
fclose(fp);
fp = fopen("myfttestr.txt", "w");
if ( fp == NULL ) {
printf("Cann't open trapzr10.txt file\n");
return 1;
}
for ( i = 0; i < ARRAYSIZE; i++ )
fprintf(fp, "%e\n", real[i]);
fclose(fp);
}
--Qiong
gsl-discuss@sources.redhat.com
--- Brian Gough <bjg@network-theory.co.uk> wrote:
> xie qiong writes:
> > Hello,
> >
> > Thank you for your reply.
> >
> > The following is my code. Even though I changed
> the
> > tolerance to 0.1, core dump still happens.
> >
> >
>
> The program does not compile (missing definitions of
> OMEGAMAX, etc).
__________________________________________________
Do You Yahoo!?
Try FREE Yahoo! Mail - the world's greatest free email!
http://mail.yahoo.com/