help-gsl
[Top][All Lists]
Advanced

[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   = &parameters;

  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




reply via email to

[Prev in Thread] Current Thread [Next in Thread]