[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] (no subject)
From: |
Francesco Abbate |
Subject: |
Re: [Help-gsl] (no subject) |
Date: |
Fri, 2 Apr 2010 18:39:21 +0200 |
2010/4/1 <address@hidden>:
> Great. Don't want to waste anyone's time...
Don't worry, your question is ok in this mailing list and I've been
myself a student in physics... :-)
> So, I am trying to use the 'qag' numerical subroutine in a program I am
> writing and am having problems passing multiple parameters from my main
> program into the function that I am trying to integrate. I was able to do
> it with one parameter, as shown on the gsl webpage, and that works ok, but
> I can't get it to work with the three I need.
>
> I tried making a structure, as shown below, and it compiles but gives a
> bus error when I try to run the code. Here is an example of what I tried:
You example is basically correct. There are only two errors:
* you forgot to set the 'params' field in the gsl_function 'ENEDE'
* in the function ENEdE you are trying to calculate a power with a
negative base because you wrote
pow((alpha-beta)*(Epeak/100),(alpha-beta))
Here a modified examples of your code that seems to work. I've
replaced the base of pow with its absolute value to avoid an error but
this is quite arbitrary...
------------------------------------
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include <gsl/gsl_integration.h>
struct parameters_{
double Epeak;
double alphaband;
double betaband;
};
const double z=2.77;
const double omega_m=0.3;
const double omega_A=0.7;
const double H0=7.1;
static double
ENEdE(double En, void * params)
{
struct parameters_ p = *(struct parameters_ *)params;
double a = p.alphaband, b = p.betaband;
double d = a - b;
if (d * p.Epeak >= En)
return En*pow(En/100, a)*exp(-En/p.Epeak);
else
return En*pow(abs(d) * (p.Epeak/100), d) * exp(-d)*pow(En/100, b);
}
int main() {
const size_t wsize = 1000;
gsl_function ENEDE;
double band1=15.;
double band2=150.;
struct parameters_ parameters = {
.Epeak = 180,
.alphaband = 0.85,
.betaband = 2.3
};
gsl_integration_workspace * w = gsl_integration_workspace_alloc (wsize);
double k1,k1_err,k2,k2_err,k3,k3_err,lumdist;
ENEDE.function = &ENEdE;
ENEDE.params = ¶meters;
gsl_integration_qag (&ENEDE, band1, band2, 0, 1e-7, wsize, 5, w,
&k1, &k1_err);
gsl_integration_qag (&ENEDE, 1/(1+z), 10000/(1+z), 0, 1e-7, 1000, 5, w,
&k2, &k2_err);
printf("%g\n",k2/k1);
gsl_integration_workspace_free (w);
return 0;
}
-------------------------------------
Otherwise you could use "GSL shell" to explore your problem before
actually writing any code. Here how you could write your example in
GSL shell:
------------------ GSL shell example ---------------------------
p = {E= 180, alpha= 0.85, beta= 2.3}
function enede(En)
local d = p.alpha - p.beta
if d * p.E >= En then
return En*pow(En/100, p.alpha) * exp(-En/p.E)
else
return En*pow(abs(d) * (p.E/100), d) * exp(-d) * pow(En/100, d)
end
end
-- make a plot
fxplot(enede, 15, 150)
-- calculate the integral, it call gsl_integration_qag internally
k1, k1err = integ {f= enede, points = {15, 150}}
---------------------------- END of GSL shell example
------------------------------------
You should save the file calling it like 'enede.lua' or something like
that and load it in gsl shell by using:
> dofile('enede.lua')
You should use GSL shell 0.9.6-2.
With the plot I've seen that your function is very smooth between 15 and 150.
I hope that helps.
Best regards,
Francesco