This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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/


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]